热防护系统分区协调耦合推进方法

黄 杰，姚卫星，陈 炎,3，孔 斌,3

(1. 南京航空航天大学飞行器先进设计技术国防重点学科实验室，南京 210016；2. 南京航空航天大学机械结构力学及控制国家重点实验室，南京 210016；3. 中国航空工业集团公司成都飞机设计研究所，成都 610091)

0 引 言

空天飞行器再入阶段受到巨大的气动加热作用，为保证内部机体结构在可承受的温度范围内，需要在机体表面附加热防护系统(Thermal protection system,TPS)[1-3]，热防护系统的设计为高超声速飞行器关键技术之一，而陶瓷防热瓦是应用最广泛的隔热结构之一，它通过应变隔离垫粘接在机体表面，其热控性能直接影响到飞行器的安全，必须进行热防护系统的热分析。

传统的热防护系统热控分析方法将气动热与结构传热人为地分割开来。文献[4-7]在已知热流密度的情况下进行了飞行器头部防热瓦、盖板式热防护系统以及多层隔热结构的热分析，获得了热防护系统和机体结构的温度-时间历程。其求解过程为先进行气动加热的分析，再将热流密度作为边界条件直接分析传热方程求解热防护系统的温度场，而忽略气动热与结构传热之间的耦合效应，这必然会影响热防护系统热控性能的分析精度。

由于空天飞行器再入过程中产生巨大的气动加热会造成结构温度急剧升高，而结构温度升高后，激波层和边界层内气体与热防护系统外表面之间的温度梯度将会减小，造成壁面热流密度的降低，即气动加热与结构温度场存在强烈的耦合效应，应给予足够的重视。

耦合分析的关键在于气动热的准确计算、耦合迭代的策略以及耦合面之间的数据插值。以往的研究常常采用Eckert[8]提出的参考焓法等工程算法进行热流密度的分析，其分析精度有限，现代气动热分析采用计算流体力学(Computational fluid dynamics,CFD)方法求解Navier-Stokes方程[9-10]，可借助CFD++或FLUENT等软件或自编程序准确分析驻点及其他区域的热流密度分布。而耦合分析常采用迭代求解法，其核心技术是要保证瞬态耦合的时间精度和迭代的协调性。此外由于结构传热网格单元的尺度远大于流体单元，耦合面节点非一一对应，故耦合面之间的数据交换需要采用插值法完成，其关键技术是要保证插值精度和耦合面热流通量的守恒性。

本文在充分考虑气动热与结构传热耦合的基础上，提出了热防护系统热控性能分析的分区协调耦合推进方法，对比了圆管算例的计算结果与试验结果，最后进行了防热瓦与主动冷却系统同时存在时的热防护系统热控性能及其影响因素分析。

1 流体和结构传热的数值算法

1.1 流体分析数值算法

不考虑体积力和内热源情况下，直角坐标系下的流体动力学N-S方程的积分形式为：

(1)

式中：W为守恒向量，Fc为对流通量，Fv为黏性通量，dS为控制体V的边界面，n为边界面dS的外法线单位向量。将式(1)按有限体积法(Finite volume method,FVM)进行空间离散可得：

(2)

式中：Wi和Vi分别为控制体i的守恒向量和体积，NF为控制体边界面的数目，ΔSN为第N个边界面的面积。由于对流通量Fc具有高度非线性特点，并且集中体现了流场的对流特征，采用具有TVD(Total Variation Diminishing)性质的NND格式[11]对其进行空间离散，其中半离散化的上风型NND格式为：

(3)

(4)

(5)

(6)

式中：sgn为符号函数，min表示取较小值。

控制方程中的黏性项采用中心差分格式离散，湍流模型采用Menter’s SST两方程模型[12]，而针对非定常问题的时间离散，在n+1时刻采用时间二阶精度的隐式三点后差离散可得到二阶精度的离散方程：

(7)

(8)

式中：Δτ和Δt分别为虚拟时间步长和物理时间步长，称为双时间步长法[13]；p和n分别为虚拟时间迭代步和物理时间迭代步。虚拟时间步上的内迭代可采用LU-SGS格式[14]求解，当p→∞时虚拟时间项趋近于零，式(8)的定常解即为二阶精度的非定常解。

1.2 结构传热数值算法

在无体积热源的假设下结构瞬态热传导的控制方程为：

(9)

式中：ρ0为结构材料密度，c为材料比热容，kx，ky和kz分别为材料三个方向的导热系数。其中比热容和导热系数一般为温度的函数。

针对本文热防护系统的热分析问题，其外表面边界条件为壁面热流密度Qaero和壁面热辐射量Qrad，其表达式分别为：

(10)

(11)

式中：∂T/∂n为壁面法向温度梯度，ε为壁面热辐射率，玻尔兹曼常数σ=5.67×10-8W/(m2·K4)，Twall为壁面温度，Tat为大气环境温度。

对式(9)进行有限元离散可得到总体合成矩阵求解方程：

(12)

(13)

求解式(13)即可得到结构各个时刻节点的瞬时温度。

2 耦合分析模型和策略

空天飞行器再入过程中热防护系统外壁面巨大的气动加热效应会造成结构温度急剧升高，而结构温度升高后，激波层和边界层内气体与壁面之间的温度梯度将会减小，造成壁面热流密度的降低，即气动加热与结构温度场存在强烈的相互耦合效应，如图1所示。

本文采用分区协调耦合推进方法进行气动热与结构传热的分析，其耦合推进策略如图2所示，即将流场与结构分为两个不同的区域，分开建模与求解，其中流场采用FVM求解，而结构热传导采用有限元法(Finite element method,FEM)求解，其基本假设和特点为：

1) 在任意n到n+1时间步内FVM(或FEM)求解过程中壁面温度(或壁面热流密度)不变，即各子学科求解时边界条件冻结。

2) 各子学科在相同的时间节点进行数据交换，以保证耦合分析的协调性以及耦合时间精度。

根据图2所示的耦合推进策略，本文提出了适用于瞬态和稳态的热防护系统热分析流程(见图3)，主要步骤为：

1) 首先建立相应的FVM和FEM数值分析模型，定义来流条件和初始温度T0，进行定常气动热分析，结果作为耦合分析的热流密度初始条件。

2) 进行第n到n+1时间步的流场求解，将计算获得的热流密度Qn+1通过插值方法传递给结构FEM模型外表面。

3) 进行第n到n+1时间步的结构FEM传热分析，输出分析获得的结构温度场Tn+1。

4) 判断时间是否达到分析的总时间ttotal(针对瞬态分析)，或壁面热流密度和结构温度场是否收敛(针对稳态分析)，若达到ttotal或收敛，结束耦合分析，否则进行热流密度和壁面温度的数据传递。

5) 返回2)，进行下一个迭代步的求解，直到分析时间达到ttotal或壁面热流密度和结构温度场收敛，结束分析。

3 耦合面数据传递

以上气动热和结构传热的耦合分析中涉及到热流密度和壁面温度的数据传递，若这两个学科之间的数值模型网格节点一一对应，那么通过节点之间的对应关系就可以方便地进行热流密度和壁面温度的数据传递。但实际分析中结构网格尺度通常远大于流体网格，即FVM和FEM耦合面节点不一致，如图4所示，故需要在耦合面进行插值以实现数据传递。

本文的分析涉及到热流密度和结构壁面温度的数据插值，其中的关键技术是要保证耦合面上热流密度通量的守恒性。

(14)

式中：Q为FVM网格面的热流密度，q为与FVM网格节点对应的结构传热FEM网格面的热流密度。当Scoup为耦合面时，式(14)体现了热流密度通量的总体守恒性。本文采用如下基于控制面的双向映射插值方法进行热流密度的数据传递，算法主要步骤为：

1) 将FVM和FEM耦合面上节点从物理空间(x,y,z)通过坐标转换映射到参数空间(u,v)中，即u=u(x,y,z)和v=v(x,y,z)，它将三维曲面上的空间节点转换到二维的控制面上。

2) 在物理空间中采用Kdtree算法搜索任意FEM网格节点ζi(x,y,z)附近的FVM网格节点ηi(x,y,z)，并将其转换到二维的控制面得到ζi(u,v)和ηi(u,v)。

3) 将FVM网格节点坐标ηi(u,v)和相应的热流Qi(u,v)代入到如下三次插值函数Q(u,v)，采用最小二乘法求解插值函数的系数bj(j=1,2,…,10)。

Q(u,v)=b1u3+b2v3+b3u2v+b4v2u+

b5u2+b6v2+b7uv+b8u+b9v+b10

(15)

4) 将所有FEM网格节点ζi(u,v)代入到已知系数bj(j=1,2,…,10)的插值函数Q(u,v)中，即可求得FEM网格节点插值热流密度qi(u,v)，并通过式(14)进行热流密度通量的守恒性检验。

4 算例1：耦合方法的验证

本文选取NASA高超声速圆管风洞试验模型[15]进行算例分析。圆管内径R1=25.4 mm，外径R2=38.1 mm，圆管材料为不锈钢，其密度ρ=8030 Kg/m3，导热系数k=16.72 W/(m·K)，比热容c=502.48 J/(kg·K)，此外圆管初始壁面温度T0=294.4 K。高超声速来流马赫数Ma=6.47，来流温度T=241.5 K，来流压强P=648.1 Pa，攻角α=0°。建立了二维分析模型，如图5所示，其中流场网格量为60×100，壁面第一层网格高度Δh=1×10-5m，并且在激波处进行了局部加密。圆管结构有限元网格量为40×25，采用四节点平面单元模拟。此外分析类型为瞬态分析，时间步长Δt=1×10-4s，总时间ttotal=2 s，并且将定常流场作为耦合分析的初始条件。

分析获得了圆管驻点温度随时间的变化曲线，如图6所示，从图中可观察到曲线斜率在初始阶段很大，随时间的推移逐渐减小，说明结构开始阶段升温很快，随后温度上升速度逐渐减缓，最后趋于一个稳定的温度值(稳态解)。

图7为2 s时刻圆管外壁面的相对温度分析情况,从图中可知壁面耦合分析结果分布曲线与试验结果分布曲线吻合良好，其中驻点处温度耦合分析结果Tstag为442 K，而试验结果为465 K，相对误差为4.95%。从而验证了分区协调耦合推进方法的正确性与分析精度，其可用于分析热防护系统气动热与结构传热的耦合问题。

5 算例2：空天飞行器头锥热防护系统耦合分析

现代空天飞行器同时采用防热瓦及主动冷却系统以保证蒙皮内表面温度稳定在设定值，如图8所示。本文进行了空天飞行器头锥热防护系统的稳态耦合分析。其中防热瓦外表面涂层热辐射率ε=0.85，LI-900防热瓦厚度h1=60 mm，应变隔离垫(SIP)厚度h2=4.5 mm，铝合金蒙皮厚度h3=1.5 mm，蒙皮导热系数k=180 W/(m·K)。防热瓦和SIP的导热系数均随温度呈非线性变化趋势[16]。飞行马赫数Ma=7，高度H=60 km，攻角α=25°，初始结构温度T0=300 K，大气环境温度Tat=247.1 K。建立了三维分析模型，流场和结构网格如图9所示，结构热分析采用八节点六面体单元模拟，蒙皮内表面温度约束在325 K，为此主动冷却系统需要持续的提供功率以保证蒙皮内表面维持在此温度。

利用本文的耦合方法进行了空天飞行器头锥热防护系统的分析。图10为非耦合和耦合方法分析获得的对称壁面处的热流密度分布曲线，可知非耦合分析方法获得的驻点热流密度比耦合方法高了53.8%。这是由于非耦合方法未考虑壁面温度升高造成激波层和边界层与壁面温度梯度减小。

图11和和12分别为非耦合和耦合方法分析得到的防热瓦及SIP外表面温度云图。其中非耦合和耦合方法获得的防热瓦最高温度分别为1108.5 K和994.1 K，相差了114.4 K；而非耦合和耦合方法获得的SIP最高温度分别为472.9 K和440.3 K，相差了32.6 K。故非耦合方法获得的防热瓦和SIP温度偏高，这是由于非耦合方法考虑壁面温度升高对气动热的反馈作用造成的。

研究了飞行马赫数Ma对防热瓦和SIP最高温度、主动冷却系统最大功率的影响，如图13所示，从图中可观察到马赫数从5增加到9，防热瓦最高温度上升了69.2%，SIP最高温度上升了36.1%，而最大冷却功率上升了266.1%。图14为以上响应随热防护系统外表面涂层热辐射率ε的变化情况，热辐射率从0.25增加到1，防热瓦最高温度降低了24.5%，SIP最高温度降低了18.2%，而最大冷却功率降低了50.1%，故应尽量采用高辐射率的涂层以降低热防护系统的厚度和主动冷却系统的功率，达到减重的目的。

在进行热防护系统的设计时，防热瓦的导热系数和厚度的选取直接影响到热控性能，是热防护系统最重要的设计参数。图15和图16分别为防热瓦和SIP最高温度、主动冷却系统最大功率随防热瓦导热系数比k/k0以及厚度h1的变化曲线，k0代表防热瓦原始导热系数。从图中可观察到防热瓦导热系数比从0.128增加到8时，防热瓦最高温度降低了3.2%，而防热瓦厚度从20 mm增加到80 mm，防热瓦最高温度上升了1.5%，即防热瓦导热系数和厚度对其最高温度影响很小。随着防热瓦导热系数的增加和厚度的降低，外部气动热能量能够更容易地传递到SIP及蒙皮，造成SIP温度升高以及主动冷却系统功率增加，并且影响显著。故在防热瓦重量约束条件下应尽量采取较厚和较低导热系数的防热瓦，以提高热防护系统的隔热性能和降低主动冷却系统的重量。

6 结 论

1) 本文针对空天飞行器热防护系统热控性能研究，提出了分区协调耦合推进分析方法，其中气动热和结构传热分别采用FVM和FEM求解，且在耦合面进行数据插值以实现数据传递。

2) 进行了圆管算例分析，圆管壁面热流密度和温度分布与试验结果吻合良好，从而验证了本文的耦合推进方法的正确性与分析精度。

3) 研究了空天飞行器头锥热防护系统的热控性能，非耦合方法获得的防热瓦和SIP温度相对耦合方法偏高，这是由于非耦合方法未考虑壁面温度升高对气动热的反馈作用，而耦合方法充分考虑了此影响。

[1] Olynick D. Trajectory-based thermal protection system sizing for an X-33 winged vehicle concept [J]. Journal of Spacecraft and Rockets, 1998, 35(3): 249-257.

[2] Shideler J L, Webb G L, Pittman C M. Verification tests of durable thermal protection system concepts [J]. Journal of Spacecraft and Rockets, 1985, 22(6): 598-604.

[3] Oscar A M, Anurag S, Bhavani V S, et al. Thermal force and moment determination of an integrated thermal protection system [J]. AIAA Journal, 2010, 48(1): 119-128.

[4] 刘双, 张博明. 发汗式主动冷却金属热防护系统主动冷却效率研究[J]. 宇航学报, 2011, 32(2):433-438. [Liu Shuang, Zhang Bo-ming. Investigation on transpiration active cooling metallic thermal protection systems [J]. Journal of Astronautics, 2011, 32(2): 433-438.]

[5] 孟松鹤, 杨强, 霍施宇, 等. 一体化热防护技术现状和发展趋势[J]. 宇航学报, 2013, 34(10):1295-1302. [Meng Song-he, Yang Qiang, Huo Shi-yu, et al. State-of-arts and trend of integrated thermal protection system [J]. Journal of Astronautics, 2013, 34(10): 1295-1302.]

[6] 魏东, 杜雁霞, 石友安, 等. 非均匀气动加热下隔热层结构的优化设计[J]. 宇航学报, 2015, 36(10):1108-1113. [Wei Dong, Du Yan-xia, Shi You-an, et al. Optimization design of heat insulation layer with non-uniform heat flow loads [J]. Journal of Astronautics, 2015, 36(10): 1108-1113.]

[7] 朱言旦, 刘伟, 曾磊, 等. 飞行器结构部件导热/辐射耦合传热特性预测方法[J]. 宇航学报, 2016, 37(11):1371-1377. [Zhu Yan-dan, Liu Wei, Zeng Lei, et al. Method of predicting conduction-radiation coupled heat transfer characteri-stics for vehicle structural component [J]. Journal of Astron-autics, 2016, 37(11): 1371-1377.]

[8] Eckert E R G. Engineering relations for friction and heat transfer to surfaces in high velocity flow [J]. Journal of Aerospace Sciences, 1955, 22(8): 585-587.

[9] Ju Y. Lower-upper scheme for chemically reacting flow with finite rate chemistry [J]. AIAA Journal, 1985, 33(8): 1418-1425.

[10] 杨肖峰, 唐伟, 桂业伟, 等. 火星环境高超声速催化加热特性[J]. 宇航学报, 2017, 38(2):205-211. [Yang Xiao-feng, Tang Wei, Gui Ye-wei, et al. Hypersonic catalytic aeroheating characteristics for Mars entry process [J]. Journal of Astronautics, 2017, 38(2): 205-211.]

[11] 沈清, 顾钢民, 高树椿, 等. NND格式在航天飞机头部段N-S方程求解中的应用[J]. 空气动力学学报, 1989, 7(2):145-155. [Shen Qing, Gu Gang-min, Gao Shu-chun, et al. Applications of the NND-scheme to the solving of Navier-Stokes equations on the fore-head of space-shuttle-like body [J]. Atca Aerodynamica Sinica, 1989, 7(2): 145-155.]

[12] Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[13] 王保国, 刘淑艳, 张雅, 等. 双时间步长加权ENO-强紧致高分辨率格式及在叶轮机械非定常流动中的应用[J]. 航空动力学报, 2005, 20(4):534-539. [Wang Bao-guo, Liu Shu-yan, Zhang Ya, et al. Pseudo-time method of weighted ENO-strongly compact high resolution schemes and its application to turbomachinery unsteady flow [J]. Journal of Aerospace Power, 2005, 20(4): 534-539.]

[14] Yoon S, Jameson A. Low-upper Gauss-Sediel method for the Euler and Navier-Stokes equations [J]. AIAA Journal, 1988, 26(9): 1025-1026.

[15] Wieting A R. Experimental study of shock wave interference heating on a cylindrical leading edge [R]. NASA-TM-100484, 1987.

[16] Ng W H, Friedmann P P, Waas A. Thermomechanical analysis of a thermal protection system with defects and heat shorts [R]. AIAA-2006-2212, 2006.