摘要
空空导弹热力环境的精确预示需要进行流、热、固多场耦合分析，并研究耦合效应对流场与结构-温度场的影响。采用分区求解方法，建立空空导弹流-热-固多场耦合模型，对导弹结构变形与温度、压力相互之间的耦合关系以及耦合效应对温度、压力仿真结果的影响进行了分析。研究表明：气动力与气动热引起弹体结构弯曲变形，且结构变形会造成导弹温度场和压力场发生变化，其中气动力作用产生的结构变形为引起温度和压力变化的主要原因。耦合效应对导弹气动热力学环境预示精度产生影响，当导弹长细比较小，或飞行速度较低，或飞行攻角较小时，耦合效应影响小；当导弹为大长细比并以高速、大攻角飞行时，考虑耦合效应预示精度高。
Abstract
With regard to accurately predicting thermodynamic environment for air-to-air missile, fluid-thermal-solid coupling analysis should be adopted, and coupling effects on flow field and structure-temperature field should be investigated. A fluid-thermal-solid multi-field coupling simulation model for air-to-air missile was established via using partition algorithm, and the coupling relationships among missile structural deformation, temperature and pressure were analyzed and the coupling effects on temperature and pressure simulation results were examined as well. Results show that, missile bending deformation is caused by aerodynamic and aerothermal, and the structure deformation could lead to the variation of missile temperature field and pressure field. The main reason for the variation is induced by the structure deformation from aerodynamic effect. The predicting accuracy of missile thermodynamic environment can be affected by the coupling effects. If the missile slenderness ratio or the angle of attack is small, or the flight speed is low, the coupling effects have a little influence on thermodynamic environment for supersonic missile. If the missile slenderness ratio and the angle of attack are large, and the flight speed is high, high predicting accuracy can be obtained by considering coupling effects.
空空导弹高速飞行时，一方面在气动热、气动力作用下弹体结构产生变形，一方面结构变形又反过来对导弹温度场和压力场产生影响，流、热、固多场之间相互作用，耦合效应明显[1-3]。精确预测超音速导弹气动热、气动力，提供结构热应力等计算初始载荷，对导弹结构设计具有指导意义，但是因为耦合效应的影响，传统仿真方法无法保证精确预示，需要采用流-热-固多场耦合分析方法[4]。
针对飞行器气动热力学环境预测，传统方法将飞行器壁面假设为绝热壁面，采用流体软件先得到壁面温度与流场分布，再将流体计算结果传递给固体进行结构热应力计算，此方法较为简单，易于实现，但未考虑流、固表面之间的热流交换，温度仿真结果需要工程系数修正[5-8]。基于此不足，在飞行器温度计算方面，夏刚等将二阶精度的AUSM+格式与结构传热的Galerk方法相结合，实现了飞行器的热-结构耦合计算[9]。屈程等将直接模拟蒙特卡罗方法与结构传热计算方法相结合，实现了飞行器流、固之间的热流交换[10]。近年来，高阶计算格式的应用也越来越广泛，如四阶 MDADF-HY 格式、五阶 WENO 格式等，可以提高飞行器气动加热仿真计算的准确性[11-12]。以上方法虽然实现了流、固之间的热流交换，但并未考虑结构变形对气动热环境的影响，当飞行器耦合效应明显时，此类方法将不再适用。在飞行器气动力环境计算方面，肖军等通过气动方程采用LU-SGS格式和隐式的线性多步法实现飞行器的气动弹性计算[13]。通过建立基于计算流体动力学（computational fluid dynamics，CFD）/计算结构动力学（computational structural dynamics，CSD）的松耦合方法，也能够实现飞行器的气动弹性载荷计算[14-15]。上述方法在进行气动力环境的仿真分析时，并未考虑气动热对飞行器的影响。晏筱璇等[16]结合自回归滑动平均（auto regressive moving average，AEMA）模型与切比雪夫方法本征正交分解（proper orthogonal decomposition-Chebyshev，POD-Chebyshev）气动热降阶模型，搭建了气动热弹性分析框架。赵建宁等[17]建立了分区耦合的间断/连续伽辽金有限元方法及其算法框架实现流-热-固耦合计算。Chen等[18]发展了一种高超声速流场热结构耦合分析平台，同时引入自适应时间步长计算方法提高计算效率。Miller等[19]开发了一种分区求解方法，重点通过研发时间推进程序实现流体热结构分析。此类多场耦合仿真分析方法能够在耦合效应明显的情况下较为精确地得到飞行器的气动热、气动力载荷，但主要在一些简单的二维模型上展开。
目前仍然缺乏对于空空导弹流-热-固相互间的耦合关系，以及耦合效应对气动热力学环境预示精度的影响的研究。而成功完成导弹设计，需要对耦合关系以及耦合效应进行深入研究，对导弹热防护以及结构强度设计具有十分重要的意义，否则会因为设计不当造成导弹飞行品质降低、结构破坏等不利后果。
本文采用流-热-固多场耦合的分区求解方法，分析了空空导弹的温度场、压力场和结构变形，研究其相互之间的耦合关系，并分析了耦合效应对气动热力学环境的影响。
1 数值计算方法
1.1 控制方程
流体计算通过求解三维Navier-Stokes方程获得：
(1)
式中，Q为守恒变量，G为无黏通量，GV为黏性通量，t为物理时间，为某一固定区域Ω的边界，dS为面积微元，n为控制边界法向单位矢量。
温度场通过数值求解三维热传导方程获得，基于能量守恒定律和Fourier定律，得到结构瞬态热传导方程：
(2)
式中，ρ为密度，c为比热容，T为温度，λ为导热系数。
应力应变通过求解热弹性力学方程，基于静气动弹性求解如下矩阵方程：
(3)
式中，K为结构刚度矩阵，KT为热传导矩阵，KuT为热弹性刚度矩阵，us为位移向量，T为温度向量，F为力载荷向量，QT为热载荷向量。
1.2 耦合计算方法
基于分区耦合方法，流-热-固耦合分析模型如图1所示。流场分析得到热流和压力，结构场分析得到壁面温度和结构位移。耦合界面进行实时参数交换，热流和压力联合影响结构热力行为，壁面温度和结构位移影响流场气动热力学行为。一个时间的基本迭代过程为：①将结构壁面温度和结构位移通过耦合界面传递给流体域，进行流体分析得到热流和壁面压力；②再将热流和壁面压力通过耦合界面传递给固体域进行热-结构分析，得到结构壁面温度和结构位移；③循环①~②步，直到流体域、固体域的计算物理量收敛。
图1多场耦合分析模型
Fig.1Multi-field coupling analysis model
1.3 仿真校验
以经典圆管绕流试验作为算例[20]。实验所用圆管内、外半径分别为25.4 mm和38.1 mm，材料参数如表1所示。圆管内壁设为等温壁，温度值与初始环境温度一致，为21.4℃。来流参数为：温度为-31.5℃，压力为638 Pa，马赫数为6.47。图2为二维计算网格，流场与圆管网格交界面为耦合界面，流固通过耦合界面进行数据传输实现多场耦合。图3和图4分别为计算初始时圆管表面压力和热流分布与文献[20]的实验结果对比，图中p0、q0分别为驻点压强与热流，θ为物面到圆心的连线与x轴的夹角。可以看出，压强p和热流q分布与文献[20]符合得很好。图5为驻点在2 s时的温度值与文献的对比，对比结果较好。由此验证了此耦合方法的有效性。
表1材料参数
Tab.1 Material parameters
图2二维计算网格
Fig.2Two-dimensional computational grid
图3壁面归一化压力分布与实验结果对比
Fig.3Comparison between normalized surface pressure distributions and experimental results
图4壁面归一化热流分布与实验结果对比
Fig.4Comparison between normalized surface heating rate distributions and experimental results
图5驻点温度随时间变化
Fig.5Variation of stagnation point temperature with time
2 仿真计算模型
以长细比（导弹弹体长度与弹径之比，即L/D）为10的无翼空空导弹为例，其几何模型如图6所示，弹体结构采用钛合金材料，材料参数见表2。
图6导弹几何模型
Fig.6Missile model
表2钛合金材料参数
Tab.2 Titanium alloy material parameters
多场耦合仿真边界条件设置如图7所示，导弹外部环境为外流场、弹体内部空间环境为内流场。假设无内部热源，不考虑辐射，对称面采用对称边界条件，其他面采用外部流边界条件。在外部流边界上进行压力、温度和速度设置，其中压力、温度基于大气高度25 km参数进行设置，即2 549 Pa、-51℃，速度基于导弹实际飞行速度设置，假设无侧偏角，仅考虑攻角，则速度在Y方向分量为0。弹体尾部采用固定边界条件，弹体结构内、外壁与流场交界面分别设置内、外场耦合边界条件。初始环境温度为20℃，初始环境压力与25 km高度大气压力相同。采用非结构化网格，图8为网格示意图，其中流场网格数约20万，导弹结构场网格数约1万。进行瞬态计算，耦合时间步长取0.000 5 s，迭代总步数为80 000，导弹飞行时间为40 s。
图7边界条件
Fig.7Boundary condition
在进行空空导弹仿真分析时，对数值仿真模型进行网格无关性验证和时间独立性检验。以长细比为20、攻角为20°、马赫数为6的最严酷飞行条件为例，弹体前端迎风侧位置作为监测点，其他仿真条件与上文一致。对A1、B1、C1、D1四种模型进行网格无关性验证，对A2、B2、C2三种模型进行时间独立性检验，得出结果分别如表3、表4所示。由表可知，当总网格数在20万以上，耦合时间步长取0.000 5 s时，仿真结果变化差异不大，所以本文计算模型网格数量满足网格无关性要求，且选取时间步长合适。
图8计算网格
Fig.8Computation grid
表3网格无关性验证
Tab.3 Grid independence verification
表4时间独立性验证
Tab.4 Time independence verification
3 仿真结果分析
以导弹长细比分别为5、10、15、20（通过改变弹身长度实现），飞行攻角分别为0°、5°、10°、15°、20°，飞行马赫数分别为3、4、5、6，飞行高度为25 km展开研究，仿真具体条件基于第2节计算模型。计算条件选取准则为：目前战斗机飞行高度在0~25 km，且在高空中空空导弹更容易实现高马赫。由此选取飞行高度为25 km，飞行马赫数最高为6。在终止迭代步，各载荷仿真值达到最大，以此时刻进行分析。
3.1 热力环境分析
因导弹在不同长细比下热力规律相似，以长细比为10的导弹弹体前端位置为例，进行研究。
图9为导弹在Ma=6时，不同攻角下的温度云图。攻角为0°时，温度沿径向均匀分布，轴向头部温度高；随着攻角的增大，沿径向温度梯度逐渐增大，且迎风面温度升高。图10为导弹在攻角为20°时，不同飞行马赫数下的温度云图。温度分布规律相似，迎风面温度高；随着导弹飞行速度的提高，径向温度梯度增大，热环境变严酷。
图9不同攻角下的温度分布
Fig.9Temperature distribution at different attack angle
图10不同马赫数下的温度分布
Fig.10Temperature distribution at different Mach number
图11为导弹在Ma=6时，不同攻角的压力云图；图12为导弹在攻角为20°时，不同飞行马赫数下的压力云图；导弹压力分布规律与温度相似。
图11不同攻角下的压力分布
Fig.11Pressure distribution at different attack angle
图12不同马赫数下的压力分布
Fig.12Pressure distribution at different attack Mach number
图13为导弹在Ma=6时，不同攻角下放大多倍的变形图。攻角为0°时，结构沿轴向拉伸；随着攻角的增大，结构沿径向弯曲且弯曲量逐渐增大。图14为导弹在攻角为20°时，不同飞行马赫数下放大多倍的变形图。导弹飞行马赫数越高，径向弯曲变形越大。
图13不同攻角下的结构变形
Fig.13Structural deformation at different angle of attack
图14不同马赫数下的结构变形
Fig.14Structural deformation at different Mach number
3.2 耦合关系分析
导弹高速飞行时，气动加热和气动力共同作用使弹体结构产生变形，尤其是弹体径向弯曲变形，引起导弹飞行攻角发生改变，反过来引起压力和温度分布变化。需要研究导弹长细比、飞行马赫数以及攻角三种因素下，温度与气动力引起的结构变形对温度、压力仿真结果的影响，有助于为导弹热防护与结构设计提供指导。将计算情况分为：热/力变形共同参与耦合计算、仅热变形参与耦合计算和仅力变形参与耦合计算。以图6弹身前端迎风面位置，即气动热力学环境严酷区域为例展开分析。基于控制替代法进行贡献率计算，定量化分析热因素（仅热变形）、力因素（仅力变形）对耦合计算的影响；贡献率越高，说明该因素引起的结构变形对仿真结果的影响越大。
贡献率计算公式为：
(4)
其中，Ca表示热因素的贡献率，Cb表示力因素的贡献率，Hab表示耦合计算结果，Ha表示热因素计算结果，Hb表示力因素计算结果。分母为0时表示热因素和力因素未引起导弹结构发生弯曲变形或弯曲变形量可以忽略，耦合效应弱，这说明各因素贡献率均为0。
温度影响分析如图15所示。以导弹马赫数为6、攻角为20°为例，研究不同长细比对温度的影响，如图15（a）所示，贡献率分析见表5。当导弹长细比为5时，热变形与力变形对温度仿真结果影响小，热、力因素不产生贡献，耦合作用弱；随着导弹长细比增加，热变形与力变形均对温度结果产生影响，且力变形影响更大（即力因素贡献大）；当导弹长细比达到20时，力变形结果与热/力变形结果温度相差不到5℃，而热变形结果与热/力变形结果相差达到21℃。
以导弹长细比为20、攻角为20°为例，研究不同飞行马赫数对温度的影响，如图15（b）所示，贡献率分析见表6。随着飞行马赫数的增加，热变形结果与热/力变形结果的温度值差异不断增大，而力变形引起的差异相对较小且增长缓慢，且力因素占主要贡献。
图15温度影响分析
Fig.15Temperature influence analysis
表5不同长细比对温度的贡献率
Tab.5 Contribution rate of different slenderness ratio to temperature
以导弹长细比为20、马赫数为6为例，研究攻角对温度的影响，如图15（c）所示，贡献率分析见表7。当导弹攻角为0°时，热、力因素不产生贡献；随着攻角增大，热变形结果与热/力变形结果的温度值差异不断增大，在攻角为20°时达到20℃；力变形结果与热/力变形结果的温度值差异也逐渐增大，但差异相对较小且增长缓慢，在攻角为20°时仅为约5℃；其中力因素占主要贡献。
表6不同飞行马赫数对温度的贡献率
Tab.6 Contribution rate of different flight Mach number to temperature
表7不同攻角对温度的贡献率
Tab.7 Contribution rate of different angle of attack to temperature
压力影响分析如图16所示。图16（a）为导弹在马赫数为6、攻角为20°下，长细比对压力的影响，贡献率分析见表8。当导弹长细比为5时，热、力因素不产生贡献，耦合作用弱；随着导弹长细比的增加，热变形与热/力变形的压力仿真结果相差达到2 kPa以上，而力变形结果与热/力变形结果压力相差仅约200 Pa，力因素贡献率大。
图16（b）为导弹在长细比为20、攻角为20°下，马赫数对压力的影响，贡献率分析见表9。随着飞行马赫数的增加，热变形结果与热/力变形结果的压力值差异不断增大，而力变形引起的差异较小，力因素占主要贡献。
图16（c）为导弹在长细比为20、马赫数为6下，攻角对压力的影响，贡献率分析见表10。当导弹攻角为0°时，热、力因素不产生贡献，耦合作用弱；随着攻角增大，热变形结果与热/力变形结果压力值差异不断增大；而力变形引起的差异相对较小且增长缓慢；力因素占主要贡献。
图16压力影响分析
Fig.16Pressure influence analysis
综上，各因素对温度和压力的影响相似。当导弹长细比小于10，或飞行马赫数低于4，或攻角小于5°时，气动加热和气动力造成弹体结构发生变形所引起的导弹的温度场和压力场改变有限；而当导弹长细比大于15且飞行马赫数达到5、攻角大于5°时，气动加热和气动力造成弹体结构发生变形所引起的导弹温度场和压力场改变较大，且此影响主要由气动力引起。
表8不同长细比对压力的贡献率
Tab.8 Contribution rate of different slenderness ratio to pressure
表9不同飞行马赫数对压力的贡献率
Tab.9 Contribution rate of different flight Mach number to pressure
表10不同攻角对压力的贡献率
Tab.10 Contribution rate of different angle of attack to pressure
3.3 耦合效应分析
考虑耦合效应能够得到导弹气动热、气动力载荷的精确解，但相比于不考虑耦合效应，计算效率降低。在实际工程应用中，允许仿真结果在许用范围内存在一定误差，当不考虑耦合效应进行仿真分析时，满足仿真结果在误差许用范围内，可忽略耦合效应进行导弹热力环境高效精准预示。研究不同长细比、飞行速度与攻角情况下耦合效应对导弹热力环境的影响，为在工程应用中选取合适的仿真计算方法提供参考，有助于高效得到导弹热防护与结构设计的初始载荷。同样以图6所示弹身前端迎风侧位置为例展开分析。
定义温度差为：
(5)
式中，Tc为考虑耦合效应的温度计算结果，Tnc为不考虑耦合效应的温度计算结果。
定义压力差为：
(6)
式中，Pc为考虑耦合效应的压力计算结果，Pnc为不考虑耦合效应的压力计算结果。
耦合效应对温度差的影响如图17所示。图17（a）为导弹在马赫数为6、攻角为20°时，不同长细比下耦合效应对温度差的影响。长细比为5时，温度差在0.1%以内；随着长细比的增加，耦合效应对温度差的影响逐渐增大，当长细比增加到20时，温度差达到约6%（25℃）。
图17（b）为导弹在长细比为20、攻角为20°时，不同飞行马赫数下耦合效应对温度差的影响，当导弹以Ma=3飞行时，温度差约2.7%（4℃）；随着飞行马赫数的增大，温度差也不断增大。
图17（c）导弹在长细比为20、马赫数为6时，不同攻角下耦合效应对温度差的影响。当导弹攻角为0°时，温度差在0.1%以内；随着攻角增加，温度差随之增大。
耦合效应对压力差的影响如图18所示。图18（a）为导弹在马赫数为6、攻角为20°时，不同长细比下耦合效应对压力差的影响。当导弹长细比为5时，压力差小于0.1%；随着长细比的增加，耦合效应对压力差的影响逐渐增大，当导弹长细比为20时压力差达到约9%（3 000 Pa）。
图17耦合效应对温度的影响
Fig.17Coupling effects on temperature
图18（b）为导弹在长细比为20、攻角为20°时，不用速度下耦合效应对压力差的影响。当导弹以Ma=3飞行时，压力差仅为1%（70 Pa）；随着飞行马赫数增大，压力差也不断增大。
图18（c）为导弹在长细比为20、马赫数为6时，不同攻角下耦合效应对压力差的影响。当导弹攻角为0°时，压力差在0.1%以内；随着攻角增加，压力差也随之增大。
综上，导弹在不同马赫数、攻角和弹体长细比的情况下耦合效应对温度差和压力差的影响相似。当导弹以低速（Ma＜4）、小攻角（α=0°，5°）飞行或者弹体长细比较小（L/D＜10）时，耦合效应对温度差和压力差的影响小，在工程分析中可考虑忽略耦合效应。当导弹高速（Ma＞5）、大攻角（α=10 °，15°，20°）飞行以及弹体长细比较大（L/D＞15）时，在工程分析中需要考虑此效应。忽略耦合效应将导致导弹的气动热力环境预示误差较大，若采用非耦合方法计算时，仿真结果需要进行工程系数修正。
图18耦合效应对压力差的影响
Fig.18Coupling effects on pressure difference
4 结论
1）空空导弹高速飞行时，头部热力学环境较为严酷，弹体结构发生轴向拉伸和径向弯曲，且随着导弹飞行速度与攻角的增大，导弹迎风面热力学环境更加严酷，弯曲变形量也随之增大。
2）空空导弹气动热力学环境的改变由气动热变形与气动力变形共同引起，且主要由气动力因素引起，且随着导弹飞行速度增高、攻角增大、弹体结构长细比加大，气动力引起的变形对导弹热力学环境影响增大。
3）导弹低速（Ma＜4）飞行，或飞行攻角较小（α=0°，5°），或弹体长细比较小（L/D＜10）时，忽略耦合效应带来的温度差与压力差均在2.7%以内，因此耦合效应对导弹的气动热力学环境预示精度影响较小，工程实际中可忽略耦合效应，以提高计算效率。
4）当导弹长细比较大（L/D＞15）并以高速（Ma＞5）、大攻角（α=10°，15°，20°）飞行时，忽略耦合效应带来的温度差与压力差分别最高达到6%与9%，耦合效应明显，在导弹热防护与结构设计中，对于初始热、力载荷的仿真计算需要考虑此效应。
5）本文目前的研究对象为无翼导弹，没有考虑舵翼对研究产生的影响。下一步，将继续深化相关研究，考虑舵翼等因素对导弹多场耦合研究产生的影响，更好地为导弹热力环境预示提供指导。
