再入飞行器端头烧蚀的耦合计算方法

龙丽平 韩俊 万田 仲峰泉, 田保未

再入飞行器端头烧蚀的耦合计算方法

龙丽平1韩俊2万田1仲峰泉1,2田保未2

（1 中国科学院力学研究所高温气体动力学国家重点实验室，北京 100190；2 中国科学院大学工程科学学院，北京 100049）

端头烧蚀是超高速飞行器再入过程中非常关心的问题。端头材料在高温高压环境中，因烧蚀其原有气动外形和结构传热边界不断变化，而气动外形和传热边界的变化又反过来影响端头热流、温度分布和烧蚀量。它们之间表现出复杂的强耦合、非线性特征。本文以碳基材料端头帽烧蚀过程为例，发展了端头帽绕流、烧蚀和结构传热耦合计算方法。通过气动、烧蚀和结构热响应计算程序的耦合和迭代，实现了对端头帽再入烧蚀过程的实时动边界模拟，并在飞行试验条件下，得到了与测量数据基本吻合的结果。

再入过程；端头帽；烧蚀；移动边界；多学科耦合计算方法

0 引言

热防护一直是高超声速飞行器及各类飞行器研制的核心技术之一。对于高马赫数再入大气层或在大气层内长时间巡航的飞行器来说，结构的热防护及其对飞行器气动性能的影响始终是一个关键的问题[1-3]。端头帽作为钝头体飞行器的头部，其气动加热环境最为恶劣，也是飞行器热防护的重点区域之一。端头帽的选材经历了难熔金属、陶瓷、石墨、C/C以及在C/C材料中添加金属材料的低烧蚀C/C等一系列材料应用过程[4-6]。国内外许多学者针对各种材料进行了烧蚀性能的研究[7-14]和烧蚀模型的发展[15-18]，对端头结构设计及优化起到了重要的指导作用。但以往报道的再入烧蚀分析大多以理论分析和解析关系为主，并且大多研究针对某一来流条件（飞行速度、高度）研究烧蚀特性，缺乏对再入烧蚀全程以及对烧蚀参数多维分布的系统性研究。因此本文旨在针对端头帽结构和环境条件，考虑烧蚀与结构传热、烧蚀外形与气动力/热的耦合作用关系，编制了气动—烧蚀—热传导耦合计算程序，实现了对再入过程中端头帽烧蚀与结构演化的全程模拟。该方法为高超声速飞行器端头帽的防热设计与优化以及评估烧蚀对飞行器气动性能的影响提供了有力工具。

1 计算流程

本文发展的端头帽外流、烧蚀以及结构传热耦合计算方法的思路如下：首先将再入过程按照飞行轨迹离散成若干准定常的飞行状态，对于每一个飞行状态及其时间，将外流与烧蚀/结构计算解耦。首先进行外流计算，将外流计算获得的边界层传热系数、摩阻系数、压力等参数作为输入条件，求解碳烧蚀速率、组分等。烧蚀与端头结构传热实时耦合计算，两者自动完成温度、热流和烧蚀量等参数的传递并实时更新烧蚀外形进行结构热传导计算。烧蚀/结构计算获得的烧蚀结构外形（由烧蚀后退量决定）、温度、碳蒸汽引射速率等参数再作为外流计算的边界条件，重新计算外流。如此，形成了外流与烧蚀/结构的耦合计算过程。

图1 绕流—烧蚀—结构耦合计算流程图

2 计算模块

2.1 气动计算

连续介质计算方法基于N-S方程求解，包含热化学非平衡效应。质量守恒方程为

振动能守恒方程为

总能守恒方程为

化学反应机理采用Park 1990模型（REF），逆向反应速率常数由公式（5）给出。

其中k为反应平衡常数，由NASA CEA程序求出。混合气体的输运系数由Blottner公式和Wilke混合模型给出，其中扩散系数采用定Lewis数假设，Lewis数采用标准值0.7，振动松弛模型采用Landau-Teller公式

控制方程采用有限体积方法离散，其中无粘项采用三阶迎风MUSCL TVD格式，限制器为minmod。粘性项采用二阶中心差分。整个方程组采用DPLR方法[19]隐式求解。湍流模型采用k-omega SST双方程模型。图2为典型再入钝头体壁面的湍流热流分布图，通常层流的热流峰值在头部驻点，而层流转捩为湍流后，热流峰值则向下游移动到肩部，这会导致钝头体肩部出现局部烧蚀量增大的后果。

2.2 烧蚀计算

本文建立了碳材料的热化学烧蚀模型，给出不同温度、压力下烧蚀速率与产物组分，发展了烧蚀与端头结构非稳态热传导耦合的求解方法。

图2 典型再入钝头体壁面的湍流热流分布图

图3 烧蚀表面热流、质量交换模型

进一步的，我们根据烧蚀表面的质量守恒条件，建立壁面的元素浓度与质量损失因子的相容关系

有关烧蚀花纹与机械剥蚀的预测模型仍在发展中，随着对烧蚀机理的深入认识与测试手段的改进，烧蚀花纹与机械剥蚀模型的准确性及适用范围将会有很大改善。

2.3 结构传热计算

结构热传导控制方程为瞬态热传导方程

边界条件

基于ABAQUS二次开发平台，同时还实现了烧蚀—热传导模块与气动计算模块的非实时耦合计算，根据烧蚀量的情况，每隔一定时间间隔进行数据交换。以下为利用本文发展的烧蚀多学科耦合计算程序测试的一个算例。

3 算例

以某飞行试验为例，端头为石墨材料，其外形如图4所示，从高度70km，以5.8km/s的速度、零度攻角再入飞行，计算总时间约28.25s。图5为轴对称结构计算模型。=0.032m，锥角9.8°。经本文的流/固耦合烧蚀程序计算，给出图6的端头帽温度随时间变化云图，图7为驻点附近各个时间段的烧蚀量云图。图8为该算例端头帽烧蚀最终外形与初始外形对比。如图8所示，烧蚀外形近似为典型的双锥结构[24]。其中烧蚀锥角为51.72度。这与该端头帽烧蚀的测量数据（锥角约48度）吻合的较好（误差为7.75%）。计算给出的驻点烧蚀后退量为0.306R（约9.79mm），测量数据为0.328R（约10.51 mm）。针对端头帽再入烧蚀的复杂问题，计算数据与试验数据取得了令人满意的相符程度。

图4 端头帽几何模型

图5 计算模型

图6 随时间变化的温度云图（单位: K）

图8 端头帽烧蚀外形与初始外形对比

从细节看，计算和实验结果的偏差主要出现在壁面附近速度大的地方，这可能是因为在肩点之后，边界层与壁面相互作用，出现了湍流的缘故，端头烧蚀花纹对于气动热流产生了放大作用，导致这部分的计算后退量相较于实验值偏小。经过分析，后退量的差异主要是表面粗糙度模型以及碳基材料高温热力学参数（如导热系数、比热容）的误差所导致的。这说明烧蚀花纹与气动计算的相互作用影响烧蚀后退量，也是本文计算方法今后需要重点改进之处。

4 结论

本文从端头烧蚀的多学科耦合物理过程出发，将自编气动计算程序、烧蚀计算程序和商业软件ABAQUS相结合，实现了基于移动网格技术的多学科耦合烧蚀计算程序，并把该程序应用来模拟某飞行试验再入过程中端头帽的烧蚀，顺利实现了气动、烧蚀、传热三个模块的数据传输和联合求解，给出了整个再入过程中端头帽烧蚀与结构传热的时间、空间演化过程，并且具有较好的计算效率和计算精度。同时，结合流动稳定性理论、材料属性等，该方法初步考虑了湍流转捩、碳—碳表面花纹演化以及机械剥落等因素。计算结果总体上与试验数据相符，基本模拟了再入过程中烧蚀外形影响气动特性，气动特性又返回来影响端头热流和烧蚀量的耦合效应。本文建立的高超声速飞行器端头烧蚀多学科耦合计算方法，较之前的单学科计算烧蚀方法有较大的进步，为高超声速飞行器端头防热设计与优化提供了有力工具。

[1] 贾洲侠, 吴振强, 吴建国, 等.飞行器气动热与结构传热双向耦合研究[J]. 强度与环境, 2019, 46(6): 16-23.[JIA Zhouxia, WU Zhenqiang, WU Jianguo, et al. Study on two-way coupled fluid-structure-thermal analysis for hypersonic vehicles[J].Structure & Environment Engineering, 2019, 46(6): 16-23.]

[2] 何钦华, 刘佳, 田玉坤, 等. 高超声速飞行器低气压环境结构热试验控制技术研究[J]. 强度与环境, 2019, 46(6): 9-15.[HE Qinhua, LIU Jia, TIAN Yukun, et al. Research on thermal test control technology of hypersonic vehicle in low atmosphere pressure environment[J]. Structure & Environment Engineering, 2019, 46(6): 9-15.]

[3] 张东, 任东鸿, 赵春章, 等. 高速飞行器流固界面热力载荷转换方法及其应用[J]. 强度与环境, 2019, 46(4): 54-58.[Zhang Dong, Ren Donghong, Zhao Chunzhang, et al. The high speed aircraft fluid-solid interface thermal load conversion and its application [J]. Structure & Environment Engineering, 2019, 46(4): 54-58.]

[4] 杨雯, 霍浩亮, 李海波, 等. 航天多功能热控材料及结构研究进展[J]. 强度与环境, 2020, 47(2): 1-12.[Yang Wen, Huo Haoliang, Li Haibo, et al. Research progress of multifunctional thermal control materials and structures of aerospace vehicles[J]. Structure & Environment Engineering, 2020, 47(2): 1-12.]

[5] 巴晓蕾, 郭文婧, 罗嘉, 等. 航天复合材料机翼疲劳试验加载技术研究[J]. 强度与环境, 2019 (3): 14-19.[Ba Xiaolei, Guo Wenjing, Luo Jia, et al. Research on loading technology of the fatigue life testing for composite wings[J]. Structure & Environment Engineering, 2019 (3): 14-19.]

[6] 高魁垠, 李海波, 吴建国, 等. 2D-C/SiC复合材料螺栓连接结构可靠性分析[J].强度与环境, 2020, 47(1): 33-40.[Gao Kuiyin, LI Haibo, Wu Jianguo, et al. Reliability analysis of 2D-C/SiC composite bolted joints[J]. Structure & Environment Engineering, 2020, 47(1): 33-40.]

[7] 黄海明, 杜善义, 施惠基. C/TaC/C复合材料的烧蚀性能研究[J]. 航空兵器, 2002, 5:5-7. [Huang H M, Du S Y, Shi H J. The ablation property of C/TaC/C [J]. Hang Kong Bin Qi, 2002, 5: 5-7.]

[8] 黄海明, 杜善义, 吴林志, 等. C/C复合材料烧蚀性能分析[J]. 复合材料学报, 2001, 18(3):76-80. [Huang H M, Du S Y, Wu L Z, et al. Analysis of the ablation of C/C composites [J]. Acta Materiae Compositae Sinica, 2001, 18(3): 76-80.]

[9] 王德文, 杨月诚, 查柏林, 等. 炭/炭复合材料耐高温和高速粒子烧蚀/侵蚀研究[J]. 弹箭与制导学报, 2013, 33(1): 132-134. [Wang D W, Yang Y C, Cha B L., et al. The study on anti-ablative property of C/C Composites [J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2013, 33(1): 132-134.]

[10] R R Williams. ATJ-S Graphite Ablation at Angle of Attack[C]. AIAA 8th Thermophysics Conference, Palm Spring, California, July 16-18, 1973.

[11] W R Grabowsky, H P Eldridge, R C Wheaton. Reentry vehicle instrumentation survey [C]. AIAA 1st Flight Testing Conference, Las Vegas, Nevada, November 11-13, 1981.

[12] J W Metzger, M J Engel, N S Diaconis. Oxidation and sublimation of graphite in simulated re-entry environments [J]. AIAA Journal, 1976, 5(8): 451-460.

[13] 尹健, 张红波, 熊翔, 等. 烧蚀条件对混合基体C/C复合材料烧蚀性能的影响[J]. 材料科学与工程学报, 2007, 25(1): 26-29. [Yin J, Zhang H B, Xiong X, et al. Influences of ablation conditions on the ablation performances of Mix-matrix C/C composites [J]. Journal of Materials Science and Engineering. 2007, 25(1): 26-29.]

[14] L W Hunter, L L Perini, D W Conn, et al. Calculation of carbon ablation on a re-entry body during supersonic/subsonic flight [J]. Journal of Spacecraft, 1986, 23(5): 487-492.

[15] 王臣, 梁军, 吴世平, 等. 高温烧蚀条件下C/C材料热力耦合场模拟[J]. 复合材料学报, 2006, 23(5): 143-148. [Wang C, Liang J, Wu S P, et al. Numerical simulation of C/C composites coupled thermos-mechanical field under the condition of high temperatures and ablation[J]. Acta Materiae Compositae Sinica, 2006, 23(5): 143-148.]

[16] 王中原, 史金光. 超高速飞行弹箭气动烧蚀数值模拟研究[J]. 南京理工大学学报, 2003, 27(5): 595-602. [Wang Z Y, Shi J G. The simulation of aerodynamic heating and ablation for a hypervelocity projectile [J]. Journal of Nanjing University of Science and Technology, 2003, 27(5): 595-602.].

[17] 王松柏. 耐热烧蚀材料烧蚀和侵蚀速度的计算方法[J]. 兵工学报, 1996, 17(4): 374-377. [Wang Songbai. The calculated method of the property and velocity of erosion with hot preservative materials [J]. Acta Armamentarii, 1996, 17(4): 374-377]

[18] 杨琼梁, 史晓鸣, 许斌, 等. 烧蚀防热层与结构耦合热传导分析的交替计算法[J]. 宇航学报, 2011, 32(8): 1854-1858. [Yang Q L, Shi X M, Xu B, et al. Alternative algorithm for heat transfer analysis of ablative protection layer coupled with structure [J]. Journal of Astronautics, 2011, 32(8): 1854-1858.]

[19] Michael J W, Graham V C, Deepak B. Data-parallel line relaxation method for the navier-stokes equations[J]. AIAA Journal, 1998, 36(9): 1603-1609.

[20] 姜贵庆, 刘连元. 高速气流传热与烧蚀热防护[M]. 北京：国防工业出版社, 2003.

[21] 刘志刚, 张巨生, 韩杰才. 炭基复合材料超高温热化学烧蚀产物分析[J]. 炭素, 2005, 124(4): 21-25.[Liu Z G, Zhang J S, Han J C. Results analyzed of super-high temperature thermal chemical ablation for carbon-based composite material[J]. CARBON, 2006, 124(4): 21-25.]

[22] R M Grabow, C O White. Surface roughness effects on nose tip ablation characteristics [J]. AIAA Journal, 1975, 13(5): 605-610.

[23] M B Ziering, V Di Cristina, Thermo-mechanical erosion of ablative plastic composites [C]. AIAA 7th Thermo-physics Conference, San Antonio, April 10-12, 1972.

[24] G J Hochrein, G F Wright. Analysis of the tater nose-tip boundary layer transition and ablation experiment [C]. AIAA 14th Aerospace Sciences Meeting, Washington D.C., 1976.

Research of Coupling Computational Method about the Nose-Tip Ablation of Reentry Vehicle

LONG Li-ping1HAN Jun2WAN Tian1ZHONG Feng-quan1,2TIAN Bao-wei2

（1 Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China；2 School of Engineering Science, University of Chinese Academy of Sciences, Beijing 100049, China)

Nose-tip ablation is a very important problem in the reentry process of hypersonic vehicle. In the high temperature and high pressure environment, the original aerodynamic shape and structure heat transfer boundary are constantly changing because of the ablation, and the change of aerodynamic shape and heat transfer boundary affects nose-tip heat flux, temperature distribution and ablation quantity in turn. They exhibit complex strong coupling and non-linear characteristics. Taking the ablation process of nose-tip of carbon based materials as an example, the coupling calculation method of flow, ablation and structure heat transfer of nose-tip are developed. Through the coupling and iteration of the calculation program of the aerodynamic, ablation and structure thermal response, the real time moving boundary simulation of the nose-tip reentry ablation process is realized, and the results are basically in agreement with the measured data under the flight test conditions.

Reentry process; nose tip; ablation; moving boundary; coupling simulation

V416.5

A

1006-3919(2021)02-0008-07

10.19447/j.cnki.11-1773/v.2021.02.002

2020-12-03；

2021-01-05

中国科学院战略性先导科技专项（XDA17030100）

龙丽平（1977—），女，博士，副研究员，研究方向：高超声速飞行器结构优化设计，可变形高超飞行器总体优化设计；（100190）北四环西路15号中国科学院力学研究所.