飞行器气动加热特性的计算与分析是高超声速飞行器设计中的关键问题之一，随着高超声速技术的迅速发展，高马赫中低空多体分离方案越来越受到重视，此设计方案的关键问题之一是分离过程中的气动热特性的快速准确预测。
目前在高超声速气动加热计算方法研究方面，针对单体外形布局或飞行器某一部件的工程方法[1, 2, 3, 4, 5]、数值方法[6, 7, 8]和工程与数值相结合的方法[9, 10]的发展已相对较成熟。随着高超声速飞行器技术的发展，多体分离问题涉及到的高度越来越低、分离时的飞行速度越来越快，多体分离过程中的气动加热问题逐渐变得重要起来，目前国内外对这一问题的研究多为实验方面[11]，因此，发展一种多体分离动态气动加热计算技术就变得非常有必要。
本文将单体飞行器布局气动加热计算技术与多体分离过程的气动热计算相结合，提出了一种多体分离动态气动加热计算技术，对多体分离过程中的气动加热问题进行计算与分析。单体飞行器布局气动加热计算技术是基于普朗特边界层理论，将整个流场分为边界层外的无粘流动和边界层内的粘性流动两部分，使用数值方法计算无粘解，气动热工程方法计算边界层传热，并耦合结构传热[12, 13]，发展了一种可用于全机外形复杂流场气动加热的计算方法；针对多体分离问题，提出了不依赖于计算网格类型和分离体个数的分离体划分策略与热防护参数设置方法，结合本文发展的气动热计算技术，形成了针对高超声速多体分离过程中各分离体表面热流密度分布以及防热层的结构温度分布的计算方法。这种方法结合了数值模拟和纯工程方法各自的优点，提高了数值模拟的效率，扩充了纯工程方法的应用范围，能方便的应用于高超声速多体分离的气动力、气动热的快速分析。
采用本文提出的数值计算技术，针对假定高超声速多体分离布局，进行了多体分离过程中的结构温度分布特性与热流密度分布特性的数值模拟，给出了动态分离过程中的计算与分析结果。 1 复杂外形气动加热计算技术
根据普朗特边界层理论，在高超声速条件下，流体的粘性和导热可以看成集中作用在物体表面的薄层即边界层内。因此，可将整个流场分为边界层外的无粘流动和边界层内的粘性流动两部分，而在使用工程方法对边界层内传热问题进行分析时，只需要边界层外缘参数的分布，由此，对气动热的求解也可以化为高超声速无粘流场的求解和边界层内粘性主导区域的气动热的求解。 1.1 无粘外流场的求解
在表面热流计算中，所关心的仅仅是边界层外缘气动特性的分布，并且在高超声速下的边界层很薄，因此本文在计算边界层外缘参数时，采用了直接求解Euler方程的方法，将Euler方程数值解的物面参数分布作为粘性流动边界层的外缘参数。
本文采用课题组已有重叠网格并行计算程序完成分离过程无粘外流解计算，输出物面网格点上的流动参数，如速度、压强、密度等，再根据完全气体状态方程，计算出其他所需流动参数，如温度、马赫数等。 1.2 表面热流密度计算
表面热流密度的计算分为驻点区和非驻点区两部分分别计算。
驻点区热流密度的计算，已有多种适用于不同条件的成熟计算方法，如：SCALA方法、LEES方法、Fay-Riddell方法[14]、Kemp-Riddl方法等。本文采用Fay-Riddell方法。Fay-Riddell公式主要用于轴对称的平衡边界层驻点热流密度的计算，公式如下：
或 其中Pr=0.71，Le=1.0。
非驻点区的表面热流计算，在工程上也有局部相似解法(LEES钝头体层流热流密度分布公式)和参考焓方法等方法，本文采用平板传热模型计算传热系数及热流密度，对较复杂外形绕流的计算进行简化处理。 1.3 结构传热计算
为了达到快速预测的目的，防热层结构温度分布根据一维热传导方程构建差分方程求解。根据防热层毕奥数的大小，防热层分为热薄壁和热厚壁分别使用不同的传热计算公式[15]。热薄壁传热公式如下：
式中ρδ为材料密度；cδ为材料比热容；δ为材料厚度；α为传热系数；ε为材料辐射系数。 2 气动加热与结构传热耦合计算
高超声速流动中，飞行器受到的气动加热过程是一个持续非稳态过程，壁面热边界条件在达到平衡前处于时刻变化中。边界层外的气流向飞行器表面传入的热流密度qw除了与飞行器气动参数有关，还与物面的温度Tw有关，而Tw的计算又与热流密度qw相关，因此，要准确模拟飞行器受热情况，气动加热和结构传热必须进行耦合计算[16]。
其中Taw为绝热壁温。
气动加热和结构传热耦合计算步骤如图 1所示。
 |
| 图 1 气动加热和结构传热耦合计算步骤Fig. 1 Calculation steps of aerodynamic heating coupled with structure heating |
在求解时间步tn时，首先以上一时间步tn－1的壁面温度计算结果Tn－1w作为tn时的物面热边界条件，计算气动热环境，得到tn时间步的热流密度qwn ；然后以热流密度qwn为热边界条件，计算结构传热，得到tn时间步的壁面温度Twn，Twn又作为下一时间步tn+1的物面热边界条件。气动热环境计算和结构传热计算交替进行，又相互影响，实现气动加热耦合结构传热。
按照上述方法，只要给定初始时刻壁温，按时间步长推进求解，就可进行高超声速全机外形复杂流动下的气动力、气动热及结构温度分布特性的快速计算。在此方面，课题组已经发展了相应的软件系统“高速飞行器气动加热计算与分析”V4.0。 3 多体分离过程气动加热计算技术
随着高超声速飞行器技术的发展，可以预测高马赫数中低空多体分离方案必将成为发展方向之一，必然涉及多体分离过程中气动加热的快速预测技术，目前国内外对此问题的数值模拟方法研究还较少。
本文在国内首次对多体分离过程中的气动加热计算技术开展了研究，发展了分离体划分策略、防热材料设置等方法，并首次完成了多体分离过程气动加热计算软件。
3.1 分离体划分策略本文提出的方法是：将飞行器的每一个分离的整体作为一个分离体，首先确定分离体的个数，然后确定每一个分离体所包含的部件序号，这些参数可在程序输入卡中进行设置，对分离体个数、计算网格类型等没有限制。
计算时，运行程序中的物面提取模块，可提取出无粘外流解中的物面部件，根据物面网格信息，完成分离体的划分。
3.2 分离体热防护参数设置方法完成分离体的划分后，整个飞行器按部件分为了不同的分离体，对每一分离体单独进行热防护参数设置，设置方法与单体热防护设置方法相同，用选择计算部件加上坐标范围划分热防护区域，进行热防护参数设置，如图 2。
|
| 图 1 热防护系统设置Fig. 1 Thermal protection system setting |
飞行器外形采用RAMC-II试验模型，总体外形为球锥体，头部半径为0.1524m，半锥角为9°，总长度为1.295m。分别进行巡航状态瞬时传热、长时加热计算。本文图中热流密度单位W/m2，温度单位K。
4.1.1 RAMC-II巡航状态瞬时传热巡航状态瞬时传热，设定巡航状态为：迎角0°、飞行高度60km，Ma=4.0、5.0、6.0，初始壁温247K。表 1给出了本文驻点热流密度计算结果与工程算法Kemp-Riddell驻点加热公式计算结果的对比，误差在10%以内。
| Ma | Kemp-Riddle/(kW·m-2) | 本文计算结果/(kW·m-2) |
| 4 | 14.23 | 14.50 |
| 5 | 27.98 | 31.41 |
| 6 | 48.98 | 53.77 |
巡航状态长时加热，设定巡航状态为：迎角6°、飞行高度60km，Ma=6.0，初始壁温247K，巡航时间1000s，飞行器防热层为2mm的Ti，表面发射率为0.8。巡航第1000s时驻点温度为945K，外表面温度分布如图 3(a)，子午面温度分布如图 3(b)，其中T_RAD为使用同样外流解直接计算辐射平衡温度结果，可以看到，子午面上温度分布与辐射平衡温度符合很好。
 |
| 图 3 巡航第1000s气动热参数分布Fig. 3 Aerodynamic heating variables distribution(t=1000s) |
由于缺乏类似算例对比，本文假定了一个机体/挂弹组合体方案。飞行方案为：组合体巡航飞行150s时释放挂弹；计算要求为：从0s开始至安全分离时刻整个过程中的机体(JT)、挂弹(DD)气动热特性。设计计算状态与初始参数设定如下：
(a) 组合体以马赫数5.0、迎角3°、飞行高度40km巡航飞行150s；计算从飞行时间0s~150s的全机气动加热特性；初始时刻组合体表面温度为250K。
(b) 第150s挂弹分离，在气动力与重力作用下按六自由度刚体运动。分离时刻参数为：弹抛出速度0m/s、弹初始姿态角3°、弹初始姿态角速度为0°/s。分离时间段为150s~160s。计算整个分离过程中机体与挂弹的气动加热特性。分离初始时刻分离体结构温度分布及气动热环境参数为为巡航段尾端第150s时刻的值。
(c) 分离体划分策略：机体和挂弹分别作为一个分离体，根据网格信息，每个分离体包含一个网格部件。
(d) 机体与挂弹采用相同的热防护方案，防热材料均为厚度为2mm的Ti，表面发射率为0.8。
 |
| 图 4 飞行器外形Fig. 4 Aircraft configuration |
图 5为飞行器巡航第150s时全表面温度分布和热流密度分布。由计算结果可知机体在头部和机翼前缘气动加热严重，温度较高，最高头部温度为1004K，挂弹在弹头处达到最高温度1088K，巡航第150s时机体头部和挂弹头部温度基本达到平衡，热流密度在头部均接近为0W/m2。
 |
| 图 5 巡航第150s全表面气动热参数分布Fig. 5 Distribution of surface aerodynamic heating variables(t=150s) |
图 6为飞行器从150s开始分离，直至安全分离整个过程中飞行器外表面温度分布和热流密度分布，分离过程中机体在头部达到最高温度1004K，挂弹在弹头达到最高温度1112K。
 |
| 图 6 分离过程全表面气动热参数分布Fig. 6 Distribution of surface aerodynamic heating variables during separation |
图 7(a)和图 7(b)分别给出了机体头部一点和挂弹头部一点的表面温度和热流密度随时间变化曲线(JT表示机体、DD表示外挂导弹)。
 |
| 图 7 机体头部和挂弹头部特征点气动热参数变化曲线Fig. 7 Variation of aerodynamic heating variables of characteristic point on JT and DD |
整个气动加热计算过程在普通台式计算机上CPU耗时约25min(不记无粘外流解计算时间)，可知本文发展的多体分离动态气动加热计算技术，可对复杂外形高速飞行器分离过程气动加热进行快速计算与分析。
4.2.2 分离过程干扰区内气动加热分离过程中，各分离体对其他分离体的流场会产生干扰，对其气动加热产生影响。
图 8为分离过程中机体与挂弹相接附近不同时刻的热流密度分布。图 9(a)为机体受挂弹干扰区内两个点上的热流密度随时间变化，Point1在左，Point2在右；图 9(b)为挂弹弹身上两个点上的热流密度随时间变化，Point3在左，Point4在右。
 |
| 图 8 分离过程机体与挂弹热流密度分布(单位：W/m2)Fig. 8 Heat flux distribution during separation |
 |
| 图 9 机体和挂弹点上热流密度随时间变化(单位：W/m2)Fig. 9 Variation of heat flux of points on JT and DD |
从图 8和图 9(a)中可以看出，机身受干扰区域内，分离开始，高热流区向后移动，Point1最大热流出现时间为151.89s，最大值7095W/m2，Point2最大热流出现时间为153.69s，最大值7668 W/m2。从图 9(b)中可以看到，弹身上两点的在分离过程中热流密度值出现跃增，峰值热流密度均超过初始时刻热流密度，其中Point4达到37505 W/m2，远远大于初始时刻热流密度22771 W/m2。
可见，分离过程中各分离体对其他分离体在气动热方面会产生不容忽视的影响。 4.3 不同热防护参数下分离过程气动加热
采用与算例4.2相同的计算外形与外流解，对机体和挂弹分别使用不同的热防护方案，对分离过程气动加热进行分析。
4.3.1 分离参数组合体巡航飞行1000s后，释放挂弹，巡航参数和分离时刻参数与算例4.2相同。
4.3.2 分离体划分及热防护区域划分
机体和挂弹分别作为一个分离体，每个分离体按坐标划分为两个热防护区域，如图 10所示，机体将头部和机身分为两个热防护区域(Tps1和Tps2)，挂弹将弹头和弹身分为两个热防护区域(Tps3和Tps4)。
 |
| 图 10 热防护区域划分Fig. 10 Thermal protection area division |
对实际热防护系统进行模拟，不同的热防护区域使用不同的热防护方案，如表 2，初始物面温度250K。
热防护区域 | 热防护类型 | 防热材料 | |
| Tps1Tps3 | 热薄壁+隔热层 | 热薄壁：Ti厚度：0.002m层数：1发射率：0.8 | 隔热层：SiO2厚度：0.02m层数：20 |
| Tps2Tps4 | 热薄壁 | 热薄壁：Ti厚度：0.002m层数：1发射率：0.8 | |
图 11(a)、图 11(b)和图 11(c)分别为分离初始时刻1000s时飞行器防热层外、内表面温度分布和热流密度分布。图 12为分离过程热流密度分布。
 |
| 图 11 分离初始时刻气动热参数分布Fig. 11 Aerodynamic heating variables distribution (t=1000s) |
 |
| 图 12 分离过程热流密度分布(单位：W/m2)Fig. 12 Heat flux distribution during separation |
本文耦合课题组发展的高超声速长时气动加热快速计算技术与多体分离数值模拟结果，完成了高超声速多体分离过程飞行器全机、部件等的气动加热特性计算技术。论文结论如下：
(1) 结合无粘外流解和气动热工程算法，耦合结构传热计算，发展了一种针对全机复杂外形的长时气动加热计算技术，完成了“高速飞行器气动加热计算与分析”软件，并在相关航空航天型号设计单位得到了成功应用。
(2) 提出了不依赖于计算网格类型和分离体个数的多体分离动态气动热加热计算方法，提出了分离体划分策略和按分离体进行热防护参数设置的方法，将多体分离和气动加热计算技术相结合，实现了分离过程中气动加热的数值模拟。
| [1] | 杨恺, 高效伟. 高超声速气动热环境工程算法. 导弹与航天运载技术, 2010, (4): 19-23. |
| [2] | HAMILTON H H, WEILMUENSTER K J, DEJARNETTE F R. Approximate method for computing laminar and turbulent convective heating on hypersonic vehicles using unstructured grids. In: 41st Annual AIAA Thermophysics Conference. AIAA 2009-4310, 2009. |
| [3] | 薛鹏飞, 龚春林, 谷良贤. 基于轴对称比拟的气动热计算方法研究. 计算机仿真, 2012, 29(6): 97-100. |
| [4] | ZHAO J S, GU L X, MA H Z. A rapid approach to convective aeroheating prediction of hypersonic vehicles. Science China Technological Sciences, 2013: 1-15. |
| [5] | 中国人民解放军总装备部军事训练教材编辑工作委员会. 高超声速气动热和热防护. 北京: 国防工业出版社, 2003: 46-136. |
| [6] | SINHA K K, REDDY D S. Effect of chemical reaction rates on aeroheating predictions of re-entry flows. Journal of Thermophysics and Heat Transfer, 2011, 25(1): 21-33. |
| [7] | 巩伟杰, 唐硕, 李世珍. 高超声速飞行器气动加热三维数值分析方法研究. 飞行力学, 2011, 29(2): 78-81. |
| [8] | PRAKASH A, ZHONG X. Numerical simulation of planetary reentry aeroheating over blunt bodies with non-equilibrium reacting flow and surface reactions. In: Proceeding of the 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition. Orlando, Florida. AIAA 2009-1542. |
| [9] | HAMILTON H H, GREENE F A, DEJARNETTE F R. Approximate method for calculating heating rates on three-dimensional vehicles. Journal of Spacecraft and Rockets, 1994, 31(3): 345-354. |
| [10] | 王杰. 高超声速飞行器气动加热计算技术. . 南京: 南京航空航天大学, 2011. |
| [11] | LIECHTY D S. Aeroheating characteristics for a two-stage-to-orbit concept during separation at Mach 6. AIAA Paper 2005-5139, 2005. |
| [12] | BOVA S W, HOWARD M A. Coupling strategies for high-speed aeroheating problems. Sandia National Laboratories, 2011. |
| [13] | 任青梅, 杨志斌, 成竹, 等. 气动加热与热响应耦合分析及试验研究. 航天器环境工程ISTIC, 2012, 29(3): 235-242. |
| [14] | FAY J A. Theory of stagnation point heat transfer in dissociated air. Journal of the Aerospace Sciences, 1958, 25(2): 73-85. |
| [15] | QJ1734-89. 战术导弹气动加热工程计算方法. 中华人民共和国航天工业部, 1989. |
| [16] | 吴洁, 阎超. 气动热与热响应的耦合研究. 导弹与航天运载技术, 2009, (4): 35-39. |