近空间飞行器拥有快速响应、突防能力强等特点，可在较短时间内对目标实施精准打击，是各个军事大国研究的重点，对国家安全、空间技术发展乃至科学技术的进步具有极大的战略意义[1-2]。在近空间飞行器高速飞行过程中，周围空气受到激波的强烈压缩，对飞行器表面产生剧烈的气动加热作用。舵翼结构作为飞行器关键部件，产生很高的温度与温度梯度，使其结构产生变形或破坏，直接影响飞行稳定性，严重时造成偏离飞行轨道甚至掉落，影响精确打击效果[3-4]。因此，亟需开展舵翼结构在气动热力结构耦合作用下的响应研究。
M.P.Galanin等利用OpenFOAM平台的RhoCentral求解器结合有限单元法构建了气动热力耦合传热迭代算法，并通过算例验证了算法的正确
等建立了导弹舵翼结构的热流固耦合分析模型，研究了在马赫数为2.3和3.7时气动力和气动热对舵翼结构变形的影响[6]。董维中等利用气动物理流场计算软件AEROPH_Flow建立了温度分布与气动热的耦合计算，实现了真实飞行条件下的高温气体非平衡效应和气动热环境的精确模拟[7]。陈梅洁等采用顺序耦合方法，基于Transition k-k1-w湍流模型计算了转捩位置、热流数值，得到了温度场和应力场的分布[8]。周印佳等采用分区求解方法，通过耦合界面的实时数据传递，实现了超高温陶瓷材料的流-热-固耦合计算[9]。李佳伟等将流场与结构温度场采用统一的控制方程进行求解，避开了传统耦合方法在迭代计算时的大量数据交换[10]。徐世南等人采用分区求解方法，开展了大长细比导弹气动力、气动热与结构变形间的耦合效应定量化研究[11]。
在实验方面，夏吝时[12]对金属尖化前缘、吴大方[13]对飞行器耐高温复合材料翼面进行了热环境下的热力性能实验。欧朝等对锥-柱-裙构件在气动环境下的结构响应开展了实验，并分析了表面粗糙度对激波的影响[14]。蔡礼港等对舵翼结构在电弧风洞中的烧蚀速率开展了实验[15]。
大多数研究者在模拟方面仅针对较低马赫数工况开展，实验方面基于单纯热环境进行强度分析，或对气动环境进行流场验证，而多物理场耦合模拟结果难以得到实验验证。本文根据集束射流风洞装置，基于计算流体力学、结构力学、空气动力学及传热学基本理论，利用ANSYS Workbench建立一体化的多物理场耦合数值模拟框架，得到舵翼结构在近空间飞行环境中的气动加热和结构响应过程，并在集束射流风洞装置中开展实验，验证数值模拟结果。
1 实验装置与数值模拟
1.1 实验装置
本文实验采用北京科技大学自主研制的集束射流气动热环境模拟实验装置如图1所示，其工作原理如图2所示。装置由气体及航空煤油供应系统、燃烧器、射流喷管、实验舱、引射管道消音塔组成。利用引射系统在实验舱内形成负压环境，空气、氧气、航空煤油在燃烧器内反应产生的高温高压气体，通过调节燃烧组分控制来流温度，经过Laval喷管的膨胀和加速作用，在喷管出口形成速度4～6马赫、总温≤2 500 K的气流，来模拟构件在近空间飞行状态下多物理场耦合作用的结构响应。
图1 集束射流气动热环境模拟实验装置图
图2 集束射流气动热环境模拟实验原理示意图
1.2 数值模拟
1) 流场控制方程
流场按照质量、动量和能量守恒方程来计算。
质量守恒方程的一般形式是由欧拉根据达朗贝尔的研究结果得到的，也被称为连续性方程：
(1)
式(1)中， ρf为流体密度(kg/m3)； t为时间(s)； ▽代表拉普拉斯算子， V为速度矢量(u、v、w)，则
(2)
式(2)中，u、v、w分别指流体微元在x、y、z三个方向的速度分量(m/s)。
流体的动量守恒方程是指动量对时间的变化率与作用在该流体微元上的各种外力之和是相等的，其一般形式为：
(3)
式(3)中， p指作用在所划分的流体微元上的压强(Pa)； fx、 fy、 fz分别指的是流体微元体积力3个方向的分量(N)； τxx、τxy、τxz等九个量为流体微元粘性应力张量的各个分量是由分子间的粘性作用引起的(Pa)。对于牛顿流体，N-S方程中的粘性应力张量的表达式为：
τxx=λ▽·V+2μux
τyy=λ▽·V+2μvy
τzz=λ▽·V+2μwz
τxy=τyx=μ(uy+vx)
τxz=τzx=μ(uz+wx)
τyz=τzy=μ(vz+wy)
(4)
式(4)中， μ为粘性系数(kg/m·s); λ为体积粘性系数(kg/m·s)，根据斯托克斯假设，通常取
针对高速流动，在空气与舵翼结构的相互作用过程中，伴随着各种形式能量的相互转换，但整个系统的总能量是保持不变的，这一规律即为能量守恒定律。能量守恒方程可表示为：
(5)
式(5)中， E为总能量(J)； e为内能，e=cvT(J)；cv为气体定容比热容(J/kg·K)； fw=fxu+fyv+fzw为单位质量的流体体积力所做的功(J)。 Qx、Qy、Qz为3个方向的压力和粘性力做的功(J)； 
为单位质量的流体的辐射热量(J)。 qx、qy、qz为3个方向的热流量，其中k为热传导系数(W/m·K)； T为温度(K)。 各部分表达式如下。
(6)
对于可压缩粘性流动，将方程(1)、(3)和(5)合称为纳维-斯托克斯方程(N-S方程)。N-S方程中有5个方程式，而有ρf、P、V(u、v、w)、T共6个未知量，未知量个数大于方程个数，则方程组无法进行求解，需添加状态方程，以使方程组封闭。
本文使用完全气体状态方程
p=Z(R*/M)ρfT
(7)
式(7)中， M表示气体摩尔质量(kg/mol)； R*为摩尔气体常数，R*=8.314 5 J/(mol·K)； Z为压缩因子，若不计气体分子之间的作用力及分子的体积，则压缩因子Z=1。
2) 结构传热方程
根据傅里叶定律和能量守恒定律，可建立热传导问题的方程，即导热微分方程：
(8)
式(8)中， T为温度场(K)； t为时间(s)； ρ为材料密度(kg/m3)； λ为材料的导热系数(W/(m·k))； C为比热容(J/(kg·K))； Q为物体内部的热源强度(W/kg)。
3) 结构热响应方程
结构受到气动力和气动加热的同时作用，气动热使结构表面温度升高，通过热量的传导使结构出现温度梯度并发生热膨胀，结构受到约束产生热应力，与气动力共同作用，使结构产生变形。结构的响应方程为：
[σ]=[D][ε]
(9)
式(9)中，[σ]为应力矩阵，[D]为结构材料弹性矩阵，[ε]为结构的总应变矩阵。
2 舵翼结构气动热力耦合模拟
飞行器在近空间飞行环境中的气动加热过程是一个持续的、非稳态的过程，并且是流场、温度场和结构场相互耦合作用的过程，本文采用流-热-固耦合解法进行气动热力耦合问题的求解。求解模拟过程如图3所示，同时建立流体域模型和固体域模型，利用基于有限体积法的Fluent求解器，通过求解连续方程、动量守恒方程和能量守恒方程，计算流体域的温度、压强、速度及耦合面上的温度分布。然后通过System Coupling模块将流体域网格节点上的热流和压力数据插值映射到固体域表面网格上，并作为结构场求解的边界条件。利用Transient Structural求解器计算结构上的温度、应力、应变与位移分布。然后再此通过System Coupling模块将结构的温度和位移场数据插值映射到流体域网格上，并以此作为边界条件进行流体域的求解，直到达到所需的耦合计算时间。
图3 流热固耦合数值求解模拟过程框图
2.1 喷管流场模拟
本文采用数值模拟和风洞实验相结合的方法证明研究结果的准确性，因此首先根据装置建立喷管流场模型，模拟得到喷管出口处的流场参数，再将该流场参数作为舵翼结构多物理场耦合计算的初始边界条件，模拟舵翼周围气动环境和结构响应。
按照实验装置实际尺寸在SolidWorks软件中建立喷管模型，使用ICEM软件进行网格划分。为保证计算结果准确高效，对模型进行结构化网格划分，由于喷管为对称结构，所以本文只建立喷管的1/2模型，网格数量为31 826，网格质量均大于0.98。喷管入口气流为氧气、航空煤油、空气燃烧后的混合气体，成分体积分数为O2 21%、CO2 11%、H2O 9%、N2 59%，喷管入口边界条件为压力入口，参数如表1所示，出口边界条件为压力出口，壁面为绝热壁面边界条件，对称面为对称边界条件，采用剪切压力传输k-w湍流模型，喷管网格如图4所示。
图4 喷管流体域网格示意图
喷管流场速度分布云图如图5所示，在喷管出口处形成了速度、温度、压强稳定的菱形区域，作为实验区域，并选择此区域流场参数作为舵翼结构多物理场耦合计算的初始边界条件，具体参数值如表1所示。
图5 喷管流场马赫数分布云图
表1 喷管流场参数
参数数值入口温度/K2 400入口压强/MPa4.07出口马赫数/Ma5.128出口温度/K513.42出口压强/Pa3 494.92出口密度/(kg·m-30.022比热容/(J·(kg·K)-1)1.33
2.2 舵翼结构建模
舵翼结构三维模型在SolidWorks软件中建立，尺寸如图6所示，取前缘处为原点建立坐标轴，下文相关位置信息以此图作为参考。将三维模型导入ICEM软件进行结构化网格划分，得到其流体域和固体域网格，如图7所示，由于在舵翼结构的外表面会存在边界层，边界层处的流场参数变化梯度较大，为能准确模拟舵翼结构与流场之间的相互作用，对边界层网格进行加密处理。网格总数为831 246，其中舵翼结构网格数为94 826，流体域网格数为736 420。采用ICEM CFD中的Determinant(2×2×2)值判断网格质量，得到网格质量均大于0.85，远大于建议值0.1，因此网格质量满足计算要求。
图6 舵翼结构示意图
图7 舵翼结构网格示意图
模型边界条件如图8所示，流体域壁面边界选择无滑移边界条件，初始壁温设置为300 K；流体域外边界设置为压力远场边界条件，选择理想气体为自由来流，来流参数选用喷管出口流场参数。为考虑舵翼结构内部传热及变形，将舵翼结构与外流场交界面设置为流-固耦合面。舵翼结构的对称面添加对称约束，尾部施加固定约束。舵翼结构材料为D6AC高强度合金结构钢，其材料参数值如表2所示。
图8 舵翼结构的边界条件示意图
表2 D6AC钢材料参数
温度/K3736739731 273比热容/(J·(kg·K)-1)523601644678导热系数/(W·(m·K)-1)43.735.519.028.5密度/(kg·m-3)7 900杨氏模量/Pa1.2×1013热膨胀系数/K-11.68×10-5泊松比0.3
2.3 模拟结果
1) 流场模拟结果
在舵翼结构多物理场耦合数值模拟中，设置时间推进步长Δt=0.000 1 s，耦合计算时间为40 s。舵翼结构20 s时外流场温度和压力分布云图如图9，由于舵翼结构头部圆弧半径很小，对来流的阻滞作用较弱，在舵翼结构头部轴线上x=-0.65 mm处形成了附体激波，压力由3 495 Pa升至109 681 Pa，气流速度急剧降低使动能转换为热能，使温度由513 K升至2397 K，激波前后温度与压力数据如图10所示。
图9 舵翼外流场温度分布和压力分布云图(20 s)
图10 舵翼结构前缘流场激波温度与压力分布曲线
如图11所示，由于舵翼前缘激波作用明显，流场升温剧烈，在1.8 s时间内由494 K升至2 300 K，最终稳定在2 397 K。舵翼结构尾部下端与夹持部分过渡不平滑，发生激波与膨胀波的掺混、引射，以及附面层的分离，贴近壁面区域的气体低速流动，外部高速流动，形成涡流，在温度监测区域形成一个局部高温区。尾部温度升温相对缓慢，在10 s时间内由300 K升至1 542 K，最终稳定在1 723 K。
图11 舵翼外流场温度变化曲线
2) 舵翼结构模拟结果
如图12(a)、图9(a)温度分布云图所示，受流场气动热作用，舵翼结构前缘驻点温度最高，沿尾部方向逐渐降低。随着时间的推进，舵翼结构温度逐渐升高，且高温区域不断向尾部移动，温度梯度逐渐减小。当t=10 s时，前缘驻点温度为1 467 K，当t=40 s时，升高到1 723 K。舵翼结构沿x方向温度分布曲线如图13(a)所示。
由图12 (c)、图9(b)所示，舵翼结构应变分布与温度分布趋势相同，前缘驻点应变最大，沿尾部方向逐渐降低。随着时间的推进，应变随温度升高而逐渐升高，且大应变区域不断向尾部移动，应变梯度逐渐减小。当t=10 s时，前缘驻点应变为17 864 με，当t=40 s时，升高到25 983 με。舵翼结构沿x方向应变分布曲线如图13(b)所示。
图12 舵翼结构温度分布和应变分布云图
图13 舵翼结构温度与应变沿x方向分布曲线
为了分析气动热和气动力对结构响应的影响大小，将气动热和气动力分别作用在舵翼结构上，变形分布如图14所示。取模拟过程中t=5 s和t=10 s两个时间点，在气动热作用下，舵翼前缘变形由0.418 69 mm增加至0.655 93 mm；在气动力作用下，舵翼前缘变形由3.095 1×10-5 mm 增加至3.149 6×10-5 mm。可见，气动热作用下结构变形值与增大幅度均远远大于气动力的作用，因此，气动热是造成舵翼结构破坏的主要因素。
3 舵翼结构风洞实验
对集束射流气动热环境模拟实验装置喷管出口处的流场参数进行测量，将气体压力测试排架安装在喷管出口，测量气流压力参数，如图15所示。
图14 气动力和气动热引起的舵翼结构变形分布云图
图15 流场压力测试实验装置
测量探头使流场产生激波，测量得到的压力为激波后的流场总压，流场实际静压值为：
(13)
式(13)中，P为测点静压值(Pa)； PT为探头测量的激波后总压值(Pa)； M为流场马赫数[16]。
实验测得激波后的流场总压为116 282.94 Pa，计算得到测点静压为3 388.46 Pa，相对模拟结果3 494.92 Pa误差为3.5%。
舵翼结构高温应变与温度测量实验如图16、图17所示，将高温应变片布置在舵翼结构尾部背风面上部、将S型热电偶布置在背风面下部，将舵翼试件安装在支架上，通过铠装高温电缆分别接入高温应变信号采集系统与温度信号采集系统。
图16 舵翼结构高温应变与温度测量实验装置示意图
图17 舵翼结构高温应变与温度测量实验装置
舵翼结构应变、温度实验测量结果与模拟曲线如图18所示，数值模拟与实验测量得到的温度数据变化曲线趋势是一致的，在时间历程内最大误差为6.41%。由于实验条件限制，在高温应变测量实验中只得到了前20 s的数据，从曲线中可以看出，数值模拟得到的舵翼结构尾部应变变化曲线与实验测得的结果是相匹配的，在20 s的时间历程内最大误差为9.73%。
图18 应变、温度实验测量与模拟曲线
4 结论
1) 流场在舵翼结构头部轴线上x=-0.65 mm处形成了附体激波，压力由3 495 Pa升至109 681 Pa；前缘在1.8 s时间内由494 K升至2 300 K，最终稳定在2 397 K；尾部形成高温涡流，在10 s时间内由300 K升至1 542 K，最终稳定在1 723 K。
2) 舵翼结构温度场与应变场规律相似，前缘驻点处最大，沿尾部方向逐渐降低，随着时间推移整体温度和应变逐渐升高，在40 s时分别达到1 723 K和25 983 με。
3) 在t=5 s和t=10 s时，气动热单独作用下舵翼前缘变形为0.418 69 mm和0.655 93 mm；气动力单独作用下，则为3.095 1×10-5mm 和3.149 6×10-5 mm。可见，气动热对舵翼结构变形的影响远大于气动力的作用。
4) 流场静压测量结果为3 388.46 Pa，较模拟结果相差3.5%；舵翼结构温度与高温应变测量结果与模拟基本一致，在时间历程内最大误差分别为6.41%和9.73%。证明本文使用的多物理场耦合数值模拟方法基本正确，建立的模型基本可行，得到的结果符合科学实验规律，对未来采用模拟方法对近空间飞行器材料构件的结构响应分析有着重要的意义。
[1] 陈雄昕，刘卫华，罗智胜，等.高超音速飞行器气动热研究进展[J].航空兵器，2014，(6)： 8-13.
[2] 彭治雨，石义雷，龚红明，等.高超声速气动热预测技术及发展趋势[J].航空学报，2015，36(1)： 325-345.
[3] 屈程，王江峰.稀薄流高超声速飞行器气动加热耦合计算[J].国防科技大学学报，2016，38(5)： 112-120.
[4] 叶友达，张涵信，蒋勤学，等.近空间高超声速飞行器气动特性研究的若干关键问题 1[J].力学学报，2018，50(6)： 1292-1310.
[5] GALANIN M P，ZHUKOV V T，KLYUSHNEV N V，et al.Implementation of an iterative algorithm for the coupled heat transfer in case of high-speed flow around a body[J].Computers & Fluids，2018，172(3)：483-491.
[6] OGNJANOVI? O V，MAKSIMOVI? S M，VIDANOVI? N D，et al.Numerical aerodynamic-thermal-structural analyses of missile fin configuration during supersonic flight conditions[J].Thermal Science，2017，21(6)： 3037-3049.
[7] 董维中，高铁锁，丁明松，等.高超声速飞行器表面温度分布与气动热耦合数值研究[J].航空学报，2015，36(1)： 311-324.
[8] 陈梅洁，程珙，田枫林，等.高超声速飞行器流-热-固耦合数值模拟[J].化工学报，2015，66(S1)： 89-94.
[9] 周印佳，孟松鹤，解维华，等.高超声速飞行器热环境与结构传热的多场耦合数值研究[J].航空学报，2016，37(9)： 2739-2748.
[10] 李佳伟，王江峰，杨天鹏，等.高超声速飞行器前缘流-热-固一体化计算[J].国防科技大学学报，2018，40(6)： 9-16.
[11] 徐世南，吴催生.高超声速导弹多场耦合仿真[J].宇航学报，2019，40(7)： 768-775.
[12] 夏吝时，杨凯威，孔维宣，等.金属尖化前缘模型主动式热疏导性能试验研究[J].装备环境工程，2017，14(10)： 82-86.
[13] 吴大方，王岳武，蒲颖，等.高超声速飞行器复合材料翼面结构 1100 ℃ 高温环境下的热模态试验[J].复合材料学报，2015，32(2)： 323-331.
[14] 欧朝，吉洪亮，肖涵山，等.MF-1 模型飞行试验结构与热防护关键问题研究[J].空气动力学学报，2017，35(5)： 742-749.
[15] 蔡礼港，常秋英，杨超，等.高超声速火箭弹舵翼气动热烧蚀预测[J].兵器装备工程学报，2018 (08)： 53-57.
[16] 傅职忠.超音速时修正空速与马赫数的换算及加速度因子的计算[J].中国民航学院学报，1996，14(3)： 10-20.