热防护系统耦合分析方法与应用

桂业伟，刘磊，杜雁霞

(中国空气动力研究与发展中心 空气动力学国家重点实验室，四川 绵阳 621000)

0 引言

热防护技术是保证航天飞行器在上升/再入段受气动加热环境下不至发生结构过热甚至烧毁的一项关键技术，也是保证导弹在气动加热环境下正常工作的一项关键技术。随着导弹技术的发展，特别是高超声速巡航弹的出现，长时间飞行导致气动力、热与结构间的多物理场耦合效应更加显著，对导弹的热防护系统进行耦合分析已势在必行。在长时间高马赫数飞行条件下，弹体结构可能在气动力载荷与结构热/力作用的综合影响下产生结构变形。变形又会影响气动加热及结构的热/力响应特性，这就形成了复杂的多场耦合关系，如图1所示。开展一体化热防护系统耦合分析方法研究，对深刻把握热防护系统多场耦合规律及其效应，准确评估热防护结构安全性，提高防热结构设计水平，从而提高导弹性能有着极为重要的意义。

图1 气动力/热与结构间的相互影响关系Fig.1 Relationship between aerodynamic force/ thermal and structure

针对防热结构内部力、热与结构高度耦合的现象，国外较早开展了热防护一体化分析方法的研究。Thornton E. A.和Dechaumphai P.在1988年就提出了多场耦合的问题[1]，并发表了关于多场耦合分析研究方面的成果。Chen[2]将耦合计算方法应用到航天领域，对再入火星大气层的探测器进行了流场和结构传热的松耦合计算，并对松耦合迭代计算中迭代次数和计算精度的关系问题进行了研究。德国科学家[3]通过未来空间运载系统技术(technology for future space transportation systems,TETRA)和未来空间运载系统应用技术(selected system and technologies for future space transportation systems applications,ASTRA)的支持，较为深入地开展了气动力/热与结构耦合问题研究，并采用如图2所示的缝隙模型验证了相关的耦合方法。

图2 耦合试验模型Fig.2 Coupling test model

国内学者在20世纪90年代初曾进行过弹头沿弹道烧蚀耦合计算方面的研究工作，以边界层理论或者其他近似计算方法沿弹道给出冷壁条件下的气动热环境，经热壁修正后用于结构烧蚀及热响应计算。这类方法属于工程方式的耦合计算，有一定的实际应用价值。20世纪90年代末至本世纪初，国内学者[4-5]率先在国内开展了气动热和热响应的耦合研究。随后研究人员[6-7]又针对飞行器防热结构，开展了气动热/热响应/热应力多场耦合一体化预测方法研究，逐步发展形成了初具雏形的热防护综合分析方法及软件系统。本文较为系统的归纳了多种热防护系统耦合分析策略与方法，并列举了耦合分析方法的实际工程应用。

1 多场耦合分析思路与策略

在计算中使用合理的耦合策略与分析方案，是多场解决问题的有效途经。本文根据气动力/热与结构多物理场间的耦合关系，提出了如下4种具有较高工程实用价值的耦合策略。通过这几种耦合分析方案的实际运用，能较好地解决目前工程中实际面临的大部分热防护系统耦合分析问题。

1.1 气动热与热响应的耦合(双向面耦合)

气动热与结构热响应之间通过飞行器表面的热量传递紧密耦合在一起，属于典型的双向面耦合。研究表明，气动热与结构热响应之间的双向面耦合可通过3种解耦方式实现：

(1) 热壁修正法

根据沿弹道冷壁热流，实时修正结构热响应计算时进入防热结构的热流值。

(2) 锚点迭代法

根据当前锚点壁温求解热壁热流，在通过热壁热流计算下一锚点的结构壁温。如此迭代推进。

(3) 整体耦合迭代法

沿弹道计算飞行器表面的冷壁热环境；以热环境为输入条件计算沿弹道的防热结构热响应情况；再根据获得的表面温度沿弹道变化结果，计算得到壁温变化情况下的热流分布情况。重复以上过程直到热环境和热响应计算结果收敛。

1.2 热响应与热应力/热变形的耦合(双向体耦合)

结构温度的非均匀上升、温度梯度的形成会使结构产生热应力，同时还会给结构带来热变形。应力/应变会对一些材料特性(接触热阻、导热率等)造成一定影响，如在气动力等共同外部因素下结构发生大变形，则会进一步影响到热响应规律。对于长时间加热的新型高超声速飞行器来说，结构热响应与热应力/热变形的耦合通常是无法忽略的。

对于结构热响应和热应力/热应变的耦合问题，计算的结构模型一致，但结构应力弹性波的传播和热量传递的时间特征尺度相差甚远。在结构发生温度改变后，热应力产生变化的时间对传热来说可忽略不计。因此，在计算分析时，热响应计算得到的温度信息与热应力计算需要的温度信息直接点点对接即可。变形场和温度场也可采用同一套计算网格，以减少三维差值误差。计算时间上，通过分析材料特性和温度梯度的大小，在典型时刻点(或等间距的时刻点)计算应力/应变场信息，并将计算结果(网格变形量、材料变化量)反馈给热响应计算模块进行迭代耦合求解。

1.3 气动力/热与结构的多场耦合(双向面体耦合)

结构变形量的大小受热壁热流、飞行器外形、结构温度与材料参数的影响，同时，飞行器气动力/热的变化量也直接受到结构变形量的大小影响。结构变形造成的飞行器外形变化对气动力/热特性造成变化属于面耦合，气动热特性的变化对热响应规律造成变化也属于面耦合，热响应规律的变化，造成结构热应力和热变形属于体耦合，因此考虑结构形变的气动力/热/结构多场耦合最为复杂。

由于涉及的物理场较多，计算量极大，耦合策略有较大难度。变形场/应力场的计算频率、气动力/气动热的更新频率将直接影响耦合效率和计算精度。研究表明，考虑飞行弹道和计算过程中结构内温度梯度的特点，计算典型状态下(或等间距的时刻点)结构变形场，并根据变形信息及时修正气动力/热信息是较好的耦合计算策略。

1.4 舱内外耦合热分析(双向面体耦合)

由于高能密度仪器设备的采用，舱内仪器设备发热对舱段结构传热的影响已不容忽视，当考虑舱内热环境对结构热响应的影响时，需建立结构热响应与舱内热环境的耦合分析方法。对于舱内多部件复杂系统，通常基于集总参数法，采用热网络节点法建立相应的温度场预测方法。结构热响应部分一般采用基于有限体积或有限元的数值计算方法。由于结构传热与舱内温度场的时间尺度通常在同一量级，因此，耦合的重点在于空间物理场的数据传递。采用高效、高精度的数据插值方法，实现热网络节点与网格点之间的数据交互，进而实现结构热响应与舱内热环境耦合分析的有效途径之一。

2 各物理场数值计算方法

2.1 气动热计算方法

(1) 气动热工程计算方法

对于有攻角飞行器气动加热，工程上通常采用跟踪流线法来计算。跟踪流线法主要思想是利用轴对称比拟概念，在弹体上布满流线，利用小横向流假设沿每根流线按等价的轴对称体处理，将三维边界层问题简化为流线坐标下的轴对称问题。层流热流用修正的Lees公式或工程方法计算，湍流热流用经过形状因子和压缩因子修正的平板湍流参考焓方法，用工程方法近似计算激波形状。舵面迎风面热流按钝前缘平板计算，前缘驻点线的热流用修正的无限后掠圆柱理论计算，干扰区的热环境通常用试验数据拟合关联方法计算。

(2) 气动热数值计算方法

三维非定常可压缩Navier-Stokes方程的直角坐标无量纲守恒形式可写为

式中：Q为守恒状态变量；E，F，G为无粘通量向量；Ev，Fv，Gv为粘性通量向量；Re为方程无量纲过程产生的无量纲系数，即通常所称的Reynolds数。

采用全层流计算，粘性系数由Sutherland公式给出，Prandtl数取Pr=0.72，对于空气的比热比计算中取常数γ=1.4。无粘通量向量采用Hanel修正的Van Leer[8-9]通量向量分裂方法计算，重构采用带有Van-Albada限制器的MUSCL[10-12]方法，定常计算时间上采用LU-SGS方法迭代进行。

2.2 舱内热环境计算方法

基于节点热网络法，舱段内部考虑导热、对流及辐射等过程的传热控制过程可描述为

2.3 热响应计算方法

直角坐标系的结构热响应控制方程可表示为

对上式在控制单元内进行积分可以得到

式中：V为单元体积；n为单元边界的外法向；N为单元的面总数(k=1,2,…,N)；Sk为单元k面的面积；qk为单元k面的热流。

时间方向采用二阶TVD-Runge-Kutta方法进行离散可以得到如下形式：

式中：i为单元序号；n为n时刻；N为单元边界面的个数；V为单元的体积；Sk为单元边界面的面积；nk为单元边界面的外法向单位向量。

2.4 热变形/热应力计算方法

对物体受热产生的应力问题，物体由于热膨胀只产生线应变，而剪切应变为0。这种由于热变形产生的应变可以看作初应变ε0，表达式为

ε0=α(φ-φ0).

在这种情况下应力应变关系就可表示为

σ=D(ε-ε0)，

式中：α为材料的热膨胀系数(1/℃)，对各向异性复合材料来说，各方向的膨胀系数α通常是不相同的；φ0为结构的初始温度场；φ为结构的稳态或瞬态温度场。φ可由温度场分析得到的单元结点温度φi插值求得。对求解域进行有限元离散，根据弹性力学变分原理就得到有限元求解方程为Ka=P。

3 应用实例与分析

3.1 高超飞行器前体考核试验状态的确定

新一代高超声速飞行器由于升阻比要求，前体要求通常采用尖锐前缘设计。这样的设计会造成更为严峻的结构温升和应力/应变问题，给考核试验提出了新要求。为更为接近真实飞行情况，考核试验通常需要选取适合的试验状态参数。本文提出的耦合分析方法，可以为考核试验提供合理的状态参数，有较大的实际应用价值。图3为飞行条件下40 s时刻前体温度分布情况。

图3 飞行条件下40 s时刻前体温度分布Fig.3 Temperature nephogram of leading edge at 40 s

计算试验状态的温度场时，不同时段的热流、焓值等参数是以阶跃的形式加载的，实际飞行中热流、焓值的变化是实时的。这里采用工程方法计算了飞行器沿轨道的表面热环境，采用三维数值方法计算局部部件沿轨道时间的温度场分布和应力场分布，得到飞行器局部部件沿轨道的温度和应力分布规律。根据该规律，并考虑地面风洞设备的实际能力，通过调整焓值和驻点压力，计算了不同时间、不同焓值下的飞行器局部试验件的温度和应力分布，并与飞行状态相对比得到了合适的试验状态参数。由图4可知，通过优化后，试验状态3的温升曲线与实际飞行的温升曲线相差较小。同时，试验状态参数3得到的热应力结果与真实弹道结果也较为接近。

3.2 高超舵翼结构的多场耦合分析

高超声速巡航飞行器通常存在舵翼结构，而舵翼结构在发生偏转或存在迎角情况时会承受较高的气动力/气动热载荷，多场耦合作用显著。本文采用气动力、气动热、热响应和热应力/变形耦合求解的双向面/体耦合方法进行计算分析。迭代过程中，气动力计算使用更新外形后的流场网格，此网格由上一轮结构变形计算获得，从而获得更新的气动力载荷。再对考虑了气动加热的机翼结构进行热应力/应变计算，以更新结构变形，直到飞行器结构变形收敛为止。

图4 试验状态和飞行状态驻点温度对比Fig.4 Stagnation temperature of test status and flight status

图6为机翼结构最高温度随时间变化曲线。从图10中可以看出，在飞行初期结构温升较快，随时间推移，结构温升趋于平缓，并逐渐达到平衡。非均匀温升会造成结构热应力，同时温升也会造成材料物性改变，从而影响气动弹性特性。因此，在气动弹性分析中考虑气动加热的耦合影响十分必要。图7所示即为舵翼结构静止状态、不考虑热影响的飞行状态和考虑热影响的飞行状态(热气动弹性)比较图。可以发现，气动加热对模型气动弹性特性和结构安全性影响显著。因此，在结构设计中应充分考虑结构温升和应力/应变的影响，合理进行结构设计和材料使用。

图5 高超声速飞行器机翼模型Fig.5 Wing model of hypersonic vehicle

图6 结构最高/最低温度随时间变化曲线Fig.6 Maximum/minimum temperature versus time of wing model structure

图7 静止/气弹/热气弹状态机翼结构形变示意图Fig.7 Deformation diagram of wing structure in the ground stationary/ aeroelastic/ aerothermoelastic state

3.3 气动力/热/轨道综合优化设计

高超声速飞行器再入大气层飞行的特点决定了其气动布局的设计优化问题是一个多学科的设计优化问题。通过本文的多场耦合分析策略与方法，可形成飞行器气动力/热与轨道综合优化设计能力。以类X-37高超声速飞行器为例(如图8)，在考虑轨道和飞行特点的基础上，选取轨道航程最大和驻点沿轨道总加热量最小作为优化子目标。多目标优化时，采用加权函数方法确定各目标间的关系。以ω1和ω2表示两目标的加权参数。几何设计变量选取机身球形头部的半径Rhd、机身中段的宽度Bw、机身的总高度Bh以及机身尾部两个扩张段的开始位置x1和x2。

图8 类X-37参数化外形Fig.8 Parametric geometrical configurations

图9 最优布局与初始布局俯视图比较Fig.9 Top view of the fuselage

图10 驻点热流密度-时间曲线Fig.10 Heat flux of the stagnation point along the trajectory

通过上述分析可以看出，对于不同的子目标加权系数，通过优化所获得的最优气动布局和飞行轨道均具有非常明显的差别。飞行轨道总航程最大和沿轨道的驻点总加热量最小是两个相互矛盾的子目标，一个性能的改善通常会给另一个性能带来一定的损失。在设计过程中，究竟选取哪种情况作为最终的设计结果，还主要依赖于整个飞行任务的需求以及设计者对于整个任务的总体把握。

4 结束语

本文主要阐述了热防护系统中遇到的多场耦合问题，提出了多场耦合分析策略与方法。并采用本文提出的耦合分析方法，对实际工程中遇到的多场耦合问题进行了计算分析。气动力/热与结构的多场耦合问题是随着热防护系统设计水平的日益提高而逐步得到认识和发展的。热防护系统耦合分析方法的建立将对以下几方面的研究起到积极作用：

(1) 为飞行器防热结构的性能准确评估提供有效工具；

(2) 为提高防热结构设计精度奠定基础；

(3) 为飞行器热气动弹性问题研究提供技术支撑；

(4) 为飞行器舱内热管理系统设计提供较为精确的计算边界。

参考文献：

[1] THORNTON E A, DECHAUMPHAI P. Coupled Flow, Thermal, and Structural Analysis of Aerodynamically Heated Panels[J]. Journal of Aircraft ,1988, 25(11): 1052-1059.

[2] CHEN Y K, HENLINE W D, TAUBER M E, Mars Pathfinder Trajectory Based Heating and Ablation Calculations[J]. Journal of Spacecraft and Rockets, 1995, 32 (2): 3-4.

[3] Matthias C Haupt, Reinhold Niesner, Ralf Unger,et al. Computational Aero-Structural Coupling for Hypersonic Applications[C]∥9th AIAA/ASME Joint Thermophysics and Heat Transfer Conference. San Francisco, California,2006.

[4] 桂业伟, 陈兰, 余泽楚,等. 窗口区域烧蚀瞬态热响应耦合计算与温升规律研究[J].工程热物理学报, 1997, 5(3): 331-335.

GUI Ye-wei, CHEN Lan, YU Ze-chu, et al. The Numerical Simulation of Coupled Unsteady Ablation and Heat Transfer in the Porthole[J]. Journal of Engineering Thermophysics, 1997, 5(3): 331-335.

[5] 桂业伟, 袁湘江. 类前缘防热层流场与热响应耦合计算研究[J], 工程热物理学报, 2002, 23(6):733-735.

GUI Ye-wei, YUAN Xiang-jiang. Numerical Simulation on the Coupling Phenomena of Aerodynamic Heating with Thermal Response in the Region of the Leading Edge[J]. Journal of Engineering Thermophysics, 2002, 23(6): 733-735.

[6] 耿湘人, 张涵信, 沈清,等. 高速飞行器流场和固体结构温度场一体化计算新方法的初步研究 [J].空气动力学学报, 2002, 20(4): 422-427.

GENG Xiang-Ren, ZHANG Han-xin, SHEN Qing, et al. Study on an Integrated Algorithm for the Flowfields of High Speed Vehicles and the Heat Transfer in Solid Structure [J].Acta Aerodynamica Sinica, 2002, 20(4): 422-427.

[7] 黄唐, 毛国良, 姜贵庆,等. 二维流场、热、结构一体化数值模拟[J].空气动力学学报, 2000, 18(1):268-372.

HUANG Tang, MAO Guo-liang, JIANG Gui-qing, et al. Two Dimensional Coupled Flow-Thermal Structural Numerical Simulation[J]. Acta Aerodynamica Sinica, 2000, 18(1): 268-372.

[8] Van Leer B. Flux-Vector Splitting for the Euler Equations[R]. ICASE Report-82-30, 1982.

[9] 朱自强. 应用计算流体力学 [M]. 北京：北京航空航天大学出版社,1998.

ZHU Zi-qiang. The Application of Computational Fluid Dynamics [M].Beijing:Beijing University Press,1998.

[10] Van Leer B. Towards the Ultimate Conservative Difference Scheme V: A Second Order Sequel to Godunov’s Method[J].Journal of Computational Physics, 1979,32:101-136.

[11] 刘磊, 桂业伟, 耿湘人,等, 高超声速飞行器热气弹静态问题研究[J].空气动力学学报, 2013, 31(5): 559-563.

LIU Lei, GUI Ye-wei, GENG Xiang-ren, et al. Study on Static Aerothermoelasticity for Hypersonic Vehicle[J].Acta Aerodynamic Sinica, 2013, 31(5): 559-563.

[12] LIU Lei, GUI Ye-wei, DU Yan-xia, et al. Study on the Similarity Criteria of Aircraft Structure Temperature/Stress/Dynamic Response[J]. Journal of Thermal Science and Technology, 2012, 7(1):1-10。