飞行器气动加热烧蚀工程计算
序言
对于高速飞行器，气动加热及烧蚀和温度场是需要考虑的。但往往需要大规模数值计算，因此，需要一套简捷快速的工程估算方法。本文提出了集成气动热、材料烧蚀、瞬态温度场的耦合计算方法。通过算例对计算方法和程序进行了验证，表明该方法具有较高的效率和精度。在给定弹道条件下，实现了气动热、热防护材料烧蚀性能和弹体温度场耦合计算，为飞行器热防护层设计提供依据。
1 引 言
飞行器在高速飞行时，会产生较大的气动热。高温气体在飞行器飞行过程中持续向机体内部传热，从而使机体尤其是头锥、翼梢前缘等部位温度升高，进而影响结构强度，造成设备失效等不良后果。对于气动加热强烈的部位，常采用烧蚀材料进行热防护。热防护材料通常由多层材料组成，其防、隔热性能直接关系到能否完成既定任务，因此，在设计阶段，预估气动热环境和热防护材料性能是必不可少的环节。
要实现对飞行器热防护材料温度场及烧蚀情况的快速预估，需要实现在一定效率和精度下的气动热、烧蚀、多层材料温度场的耦合计算。国内已有研究人员对气动热[1-4]、烧蚀机理[5-8]等方面展开了不同程度的研究，但仍然缺乏一套高效并具有一定精度的，能够实现气动热环境下飞行器热防护层性能快速预报及辅助设计的方法。
气动热计算主要有两类方法：一是直接对N-S方程进行数值计算，二是采用工程算法。基于求解N-S方程的气动热算法虽然精度较高，但计算效率低，且对设计初期，气动外形还未完全确定的情况下，该方法并不适用。而工程算法具有较高的可信度和效率，只需给定特征外形和飞行参数即可预估气动热。为此，本文将飞行弹道作为输入参数，集成气动热计算的工程算法（也保留通过CFD或实验得到的热流数据输入接口）、材料烧蚀模型和传热计算，实现了具有较高效率和精度的气动热环境下热防护材料烧蚀传热的计算方法和软件。可对输入的不同飞行参数和多层材料组合进行快速评估，达到辅助设计的目的。
2 计算方法
以飞行弹道为输入参数的气动热烧蚀计算主要包括三部分：气动热计算，材料烧蚀计算，多层材料瞬态温度场计算。
2.1 气动热工程计算
目前常用的高超声速气动热工程算法综合考虑了外流场和实验数据，不仅计算效率高，而且对于飞行器的典型几何特征位置（如钝头旋成体，翼梢前缘等部位）的驻点、球头表面和锥体表面等具有较高的可信度。
驻点热流的工程计算方法采用Fay-Riddell公式[9]。该式是依据相似性假定，对高温气体的边界层方程进行简化得到的，适用于来流总焓hs在1549-24158kJ/kg之间，壁温Tw在300-3000K之间的计算：
\({{q}_{ws}}=0.763{{\Pr }^{-0.6}}{{\left( \frac{{{\rho }_{w}}{{\mu }_{w}}}{{{\rho }_{s}}{{\mu }_{s}}} \right)}^{0.1}}\sqrt{{{\rho }_{s}}{{\mu }_{s}}{{\left( \frac{d{{u}_{e}}}{dx} \right)}_{s}}}\left[ 1+\left( L{{e}^{0.52}}-1 \right)\frac{{{h}_{D}}}{{{h}_{s}}} \right]\left( {{h}_{s}}-{{h}_{w}} \right)\) (1)
球头表面热流密度按照Lees简化公式[10]进行归一化计算：
\(\frac{{{q}_{wl}}}{{{q}_{ws}}}=\frac{2\theta \sin \theta \left[ 1-\left( \frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}} \right) \right]{{\cos }^{2}}\theta +\left( \frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}} \right)}{D{{\left( \theta \right)}^{1/2\;}}+\frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}}}\) (2)
其中：
\(D\left( \theta \right)=\left( 1-\frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}} \right)\left( {{\theta }^{2}}-\frac{\theta \sin 4\theta }{2}+\frac{1-\cos 4\theta }{8} \right)+\frac{4}{{{\gamma }_{\infty }}Ma_{\infty }^{2}}\left( {{\theta }^{2}}-\theta \sin 4\theta +\frac{1-\cos 2\theta }{2} \right)\) (3)
式（2）中，qwl为球头表面热流，θ为球头圆心角。
锥体热流密度qwr同样通过归一化计算，其公式仅与驻点热流qws和来流参数有关：
\(\frac{{{q}_{wr}}}{{{q}_{ws}}}=A\left( {{\theta }_{c}} \right)\frac{{{{s}’}}/{{{R}_{N}}}\;}{{{\left[ B\left( {{\theta }_{c}} \right)+{{\left( {{{s}’}}/{{{R}_{N}}}\; \right)}^{3}} \right]}^{\frac{1}{2}}}}\) （4）
其中：
\(A\left( {{\theta }_{c}} \right)=\frac{\sqrt{3}}{2}{{\left[ \left( 1-\frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}} \right){{\sin }^{2}}{{\theta }_{c}}+\frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}} \right]}^{\frac{1}{2}}}\sqrt{\frac{\pi }{2}-{{\theta }_{c}}}\) （5）
\(B\left( {{\theta }_{c}} \right)=\frac{\frac{3}{16}}{{{\sin }^{2}}{{\theta }_{c}}{{\left[ \left( 1-\frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}} \right){{\sin }^{2}}{{\theta }_{c}}+\frac{1}{{{\gamma }_{\infty }}Ma_{\infty }^{2}} \right]}^{\frac{1}{2}}}}\times \frac{D\left( \theta \right)}{\theta }\) (6)
在给定飞行弹道后，通过飞行马赫数，攻角和大气参数计算出来流参数，即可快速计算出特征部位的气动热。
2.2 烧蚀计算
目前飞行器常用烧蚀防热材料有硅基材料和碳基材料两种。不同类型的烧蚀材料其烧蚀模型略有不同，但本质仍然是依据质量守恒和能量守恒。硅基材料的烧蚀模型如下。
本文算例采用的高硅氧/酚醛材料是典型的硅基材料。其烧蚀机理为：当硅基材料开始烧蚀后，表面不断向后退移，如图1所示。原始材料在气动加热下，温度逐渐升高，当温度达到热解温度后，材料中的树脂成分开始热解、蒸发并进入气体边界层；随着树脂材料的热解和蒸发，剩余材料的温度在热流的作用下持续升高，到达某一温度后开始碳化，形成碳化层；最靠近外部的壁面温度最高，材料中的SiO2等成分处于液体态，并形成液态层；在发生烧蚀的情况下，高速气流会对材料表面产生很高的剪切力，从而导致材料烧蚀表面被不断吹除，使烧蚀材料表面不断向内退移。通过材料的热容、熔化、蒸发、吸热反应、引射、热阻塞效应等吸热
或隔热。
硅基材料烧蚀机理主要由质量损失机理和吸热机理两部分构成。质量损失机理是指：随着温度升高，材料中的树脂成分首先发生热解反应，生成碳和气体，表面的碳与氧气反应燃烧，生成二氧化碳等气体带走质量；随着温度继续升高，材料中的高硅氧纤维开始软化、熔融，在材料表面形成一层薄SiO2液态层，一部分蒸发气化被带走，另一部分在气动剪切力的作用下被吹除。吸热机理是指：材料热容吸热；树脂成分热解吸热；高硅氧纤维熔化吸热；SiO2液态层蒸发吸热、热解气体和蒸发气体引射引起的热阻塞效应；碳燃烧放热；材料表面热辐射；材料被气动力吹除带走热量。
图1 高硅氧/酚醛复合材料烧蚀模型及X射线拍摄的三个区域图片[11]
根据上述两个机理得到烧蚀计算的能量守恒和质量守恒方程：
\(\Psi {{q}_{w}}-{{q}_{r}}+{{q}_{rs}}=\rho {{V}_{-\infty }}\left[ {{C}_{p}}\left( {{T}_{w}}-{{T}_{p}} \right)+{{C}_{p}}\left( {{T}_{p}}-{{T}_{0}} \right)+{{f}_{p}}\Delta {{h}_{p}}+{{f}_{c}}\Delta {{h}_{c}} \right]+{{\alpha }_{\text{Si}{{\text{O}}_{2}}}}{{f}_{\text{Si}{{\text{O}}_{2}}}}{{\rho }_{1}}{{V}_{-\infty }}\Delta {{h}_{v}}\) (7)
\(\left( 1\text{-}{{f}_{p}} \right){{\rho }_{1}}{{V}_{-\infty }}={{\alpha }_{\text{Si}{{\text{O}}_{2}}}}{{f}_{\text{Si}{{\text{O}}_{2}}}}{{\rho }_{1}}{{V}_{-\infty }}+B{{\rho }_{1}}\delta _{L}^{2}\left( \frac{\tau -2P{{\delta }_{L}}}{\mu w} \right)\) (8)
其中ψ为引射因子，δL为表面液态层厚度。
\(\Psi =1\text{-}\beta {{\left( \frac{{{M}_{a}}}{{{M}_{j}}} \right)}^{\eta }}\left( {{\alpha }_{\text{Si}{{\text{O}}_{2}}}}{{f}_{\text{Si}{{\text{O}}_{2}}}}{{\rho }_{1}}{{V}_{-\infty }}+{{f}_{p}}{{\rho }_{1}}{{V}_{-\infty }} \right)\) (9)
\({{\delta }_{\text{L}}}=\frac{{{\lambda }_{L}}{{T}_{W}}}{\mu \left( \Psi {{q}_{w}}-{{\alpha }_{\text{Si}{{\text{O}}_{2}}}}{{f}_{\text{Si}{{\text{O}}_{2}}}}{{\rho }_{1}}{{V}_{-\infty }}\Delta {{h}_{V}} \right)}\) (10)
\({{V}_{-\infty }}\)为烧蚀速率。用迭代法对方程（7）-（10）进行计算，得到烧蚀温度和烧蚀速率，然后以TW作为温度边界条件，\({{V}_{-\infty }}\)作为边界运动速度，根据材料内部的传热计算，就可计算出瞬态温度场。
2.3 多层材料温度场计算
对于大面积防热来说，沿物面方向温度梯度远小于厚度方向的温度梯度，因此采用一维传热可获得局部的传热情况。温度场计算控制方程为一维非稳态导热方程：
\(\frac{\partial T}{\partial x}=\lambda \frac{{{\partial }^{2}}T}{\partial {{y}^{2}}}\) (11)
其数值求解过程如下：首先是对求解区域进行离散化，即对于m层材料组成的组合壁面，将每层材料分别任意等分为N(i),i=1,m个节点，然后对每个节点列出对应的差分方程，最后联立求解各个节点的方程组，就可以得到各点的瞬时温度。
图2 多层材料离散
方程求解采用FTCS隐式差分格式，该格式构造简单，计算效率高，且无条件稳定。
对于任意一层内部材料i，内部节点j，差分方程写为：
\(\frac{{{\lambda }_{i}}}{\Delta {{y}_{i}}}\left( T_{j-1}^{n}-T_{j}^{n} \right)\Delta t-\frac{{{\lambda }_{i}}}{\Delta {{y}_{i}}}\left( T_{j}^{n}-T_{j+1}^{n} \right)\Delta t={{\left( {{C}_{p}}\rho \Delta y \right)}_{i}}\left( T_{j}^{n}-T_{j}^{n-1} \right)\) (12)
外边界条件包括热流、辐射换热和空气电离产生的热流：
\({{q}_{w}}-{{q}_{r}}+{{q}_{E}}=\frac{{{\lambda }_{i}}}{\Delta {{y}_{i}}}\left( T_{j}^{{}}-T_{j+1}^{{}} \right)+\frac{{{\left( {{C}_{p}}\rho \Delta y \right)}_{i}}}{2\Delta t}\left( T_{j}^{{}}-T_{j}^{0} \right)\) (13)
需要注意的是，当j位于两层材料的交界面时，控制体包含两种材料，其差分方程为：
\(\frac{{{\lambda }_{i}}}{\Delta {{y}_{i}}}\left( T_{j\text{-}1}^{{}}-T_{j}^{{}} \right)\text{=}\frac{{{\lambda }_{i}}}{\Delta {{y}_{i}}}\left( T_{j\text{+}1}^{{}}-T_{j}^{{}} \right)\text{+}\frac{{{\left( {{C}_{p}}\rho \Delta y \right)}_{i}}\text{+}{{\left( {{C}_{p}}\rho \Delta y \right)}_{i\text{+}1}}}{2\Delta t}\left( T_{j}^{{}}-T_{j}^{0} \right)\) (14)
3 计算验证
按照上述公式及数值方法编写了相应的耦合计算软件，该软件可以在输入弹道（或热流）后，计算飞行器特征部位任意指定位置的热流、多层材料温度场以及烧蚀率，并可实时修改防热层的层数、厚度等相关参数，便于优化设计。
首先将数值计算值与实验值进行对比，验证计算的准确性，再通过某导弹的假想飞行弹道，给出本文所述方法的具体应用。
3.1 热流计算
首先对几个不同高度和不同速度条件下的热流进行计算，并与参考实验数据[12]进行对比。图3给出了高度7.6km、21.4km和36.6km时对应不同速度的驻点热流，可见本文的计算与实验数据符合较好。
图3 驻点热流密度
图4给出了弹头RN=6.6mm球头热流密度的分布，在相同状态下与参考实验[13]数据进行对比，可见当\(\theta \prec 60{}^\circ \)符合较好，而在\(\theta \succ 70{}^\circ \)后，计算值误差较大。
图4 球头热流密度
3.2 传热计算验证
因商业软件Abaqus的无烧蚀传热分析精度得到业界认可，故本文采用该软件对多层材料温度场计算进行了验证。算例采用了热防护系统通用结构形式，由四层材料组成,最外层为第一层，第一层为高硅氧/酚醛复合材料，第二层为铬镍钢，第三层为微球形碳泡沫[14]，第四层为弹体材料铬钢。材料的物性参数和厚度见表1。
外边界条件为q0=200kW/m2，其他表面为绝热边界条件，初温300K，时长40s。
表1 材料物性参数和厚度
| Layer&
Name |
kg/m3 | Cp
J/(kg·K) |
W/(m·K) | thickness
mm |
| 1.silica/p-henolic | 1650 | 1100 | 0.5 | 6 |
| 2.chromi-um nickel steel | 7820 | 460 | 20 | 2 |
| 3.carbon form[14] | 600 | 700 | 0.3 | 4 |
| 4.chromium steel | 7750 | 470 | 36.1 | 2 |
图5对比了本文计算值和Abaqus计算值，符合较好。
图5 温度场对比
3.3 烧蚀模型验证
烧蚀模型验证以试验数据作为对比，计算中，热流、恢复焓、参考压力等与参考文献[15]的实验一致，高硅氧/酚醛部分物性参数取值见参考文献[16-19]。表2对比了计算结果与实验数据，其中材料烧蚀率计算值与实验值的平均误差为3.6%，说明计算和实验符合较好。
表2 烧蚀数值计算与电弧加热实验对比
| No | kW/m2 | kJ/kg | (MPa) | ||
| Cal | test[15] | ||||
| 1 | 7105.4 | 6740.8 | 0.161 | 0.29 | 0.3 |
| 2 | 11664 | 11053 | 0.161 | 0.5 | 0.477 |
| 3 | 15217 | 14235 | 0.163 | 0.61 | 0.59 |
| 4 | 6996.1 | 8667.1 | 0.154 | 0.31 | 0.32 |
| 5 | 12987 | 12058 | 0.152 | 0.52 | 0.5 |
| 6 | 15219 | 14987 | 0.155 | 0.62 | 0.6 |
4 应用算例
以某导弹的假想飞行弹道计算气动热和热防护材料在热环境下的烧蚀情况和温度场分布。导弹的初速度V0=2km/s，初始飞行高度为h=40km，终止飞行高度h=0km，飞行时长37s。在飞行过程中，俯仰角和攻角变化不大，俯仰角θ=30o，攻角α=0。图6所示为飞行速度和高度，由于飞行速度高，飞行高度基本呈线性变化。
图6 飞行速度和高度
导弹为钝头体，球头半径RN=6.6mm，半锥角=20o。弹头表面热防护材料为4mm高硅氧/酚醛复合材料，弹体材料为铬钢，材料物性参数见表1，初温25oC。计算中，大气参数根据飞行高度由标准大气模型确定。
通过本文集成的计算方法和程序，计算了弹头指定位置的热流、热防护材料烧蚀厚度和温度场分布随时间的变化。以下给出驻点、球头θ=30°和锥段（当地弹体半径）R1=10mm处的计算结果。
图7为以上三个部位沿飞行弹道的热流密度，可见随着飞行速度的加快和高度下降（大气压力及密度升高），热流密度也随之逐渐升高。飞行速度在飞行28s左右达最大值，随后由于阻力升高，速度迅速下降，热流密度也随之快速下降。
图7 计算部位热流密度
图8为对应位置的烧蚀情况。驻点的烧蚀厚度约为1.54mm，球头θ=30°处烧蚀厚度为0.65mm，锥段烧蚀厚度约为0.2mm。由图8可知，在飞行的初始阶段，由于飞行高度较高，空气稀薄，热流强度不足以使表面材料升高到烧蚀温度（约1400oC）；从20s左右开始，热防护材料发生烧蚀，随着时间推移，烧蚀厚度不断增加，在27s左右烧蚀率最大；随着飞行速度下降，热流密度减小，烧蚀率又逐渐降低，在35s左右，锥段烧蚀基本停止，随后驻点和球头的烧蚀也逐渐停止。对应热流密度可知，烧蚀发生在热流密度约4500kW/m2以上的热环境下。
图8 计算部位烧蚀厚度
图9为导弹飞行37秒后，弹头的温度分布。其中，驻点温度最高，约1700oC，球头30o处约1600oC，锥段则为1500oC左右。因为弹体材料为铬钢，其导热系数较大，在厚度方向上基本为等温分布，当热防护层的初始厚度为4mm时，弹体内壁温度升高到约300oC，满足不了热防护的要求。因此，将热防护材料厚度增加到6mm，重新进行计算，得到同一时刻的弹体内壁温度约为170oC，可满足设计要求。
图9 计算部位温度场分布
5 结论
本文集成了气动热工程算法、材料烧蚀计算和传热计算，形成了一种飞行器气动热及多层防热材料烧蚀温度场预报的方法。通过算例验证了该方法具有较高的效率和可信度。采用该方法可以快速进行给定弹道条件下的气动热、热防护材料烧蚀性能和弹体温度场计算，为飞行器热防护层设计提供依据。
参考文献
[1]吕丽丽, 高超声速气动热工程算法研究 [D]. 西北工业大学,2005.2.
[2] 杨凯, 高效伟, 高超声速气动热工程算 法[J]. 导弹与载人航天运载技术. 2010.4.
[3] 薛鹏飞, 龚春林, 谷良贤, 基于轴对称比拟的气动热计算方法研究[J]. 计算机仿真, 2012.6.
[4] 李建林, 唐乾刚, 霍林等[J]. 国防科技 大学学报, 2012.12.
[5] 欧东斌, 陈连忠, 陈海群等, 硅基复合材料烧特性试验研究[J]. 宇航材料工艺, 2009.1.
[6] 孙冰, 林小树, 刘小勇等. 硅基材料烧蚀模型研究[J]. 宇航学报, 2003.5.
[7] 张俊华, 李锦文, 魏化震等, 高性能酚醛树脂基烧蚀复合材料的研究[J]. 纤维复合材料, 2009.3.
[8] 易法军, 猛松鹤, 梁军等, 防热复合材料高温体积烧蚀模型[J]. 哈尔滨工业大学学报, 2001.12.
[9] Fay,J.A. and Riddell,F.R, Theory of stagnation point heat transfer in dissociated air[J],Journal of Aeronautic Science,1958.25(2):73-85.
[10] Lees, L., Laminar heat transfer over blunt nosed bodies at hypersonic flight speed, Jet Propulsion,1956,26(4):259-259.
[11] 时圣波，高硅氧/酚醛复合材料的烧蚀 性能与热力学行为研究[D]，哈尔滨工业大学，2013.5.
[12] Rose, P.H and Stark,W.I., Stagnation point heat transfer measurements in dissociated air[J],Journal of Aerospace Science,1959:421-430.
[13] Kemp, N.H.,Rose, R.H. and Detra R.W., Laminar heat transfer around blunt bodies in dissociated air, Journal of the Aerospace Science,1959:421-430
[14]Bruneton E, Tallaron C, Gras-Naulin N,et al. Evolution of the structure and mechanical behaviour of a carbon foam at high temperatures[J]. Carbon,2002.40:1919
[15]姜贵庆，李仲平，俞继军，硅基复合材 料粘性系数的参数辨识[J]，空气动力学学报，2007.8
[16]姜贵庆，刘连元，高速气流传热与烧蚀热防护[M]，北京，国防工业出版社，2003:52-66
[17]张志成，潘海林，刘初平等，高超声速气动热和热防护[M],北京，国防工业出版社，2003:222-248
[18] Henderson J B, Verma Y P, Tant M R, Moore G R.Measurement of the Thermal Conductivity of Polymer Composites to High Temperatures Using the Line Source Technique[J]. Polymer Composites. 1983, 4(4): 219-224.
[19] Henderson J B, Tant M R. A Study of the Kinetics of High-temperature Carbon-Silica Reactions in an Ablative Polymer Composite[J]. Polymer Composites. 1983, 4(4): 233-237.
作者简介：
孙得川，大连理工大学，教授；所学专业：航空宇航推进理论与工程。主要从事以下领域研究：
- — 多孔介质中流动与传热的直接模拟方法
- — 喷嘴喷雾过程的建模及仿真
- — 多孔介质中推进剂分解反应过程建模
- — 低温推进剂相变过程的建模及仿真
- — 气动热环境下的材料烧蚀过程建模
- — 喷气推进系统设计的新思路、新方法
- — 气动热环境模拟的天地一致性
- — 计算流体力学的新方法及其实现
