高超声速飞行器表面温度分布与气动热耦合数值研究

董维中, 高铁锁, 丁明松, 江涛, 刘庆宗

中国空气动力研究与发展中心， 绵阳 621000

高超声速飞行器表面温度分布与气动热耦合数值研究

董维中*, 高铁锁, 丁明松, 江涛, 刘庆宗

中国空气动力研究与发展中心， 绵阳 621000

针对高超声速飞行器热防护设计中的高温气体非平衡效应问题和气动热环境精确预测问题，基于流场的非平衡Navier-Stokes方程、表面的能量守恒方程和内部的热传导方程，考虑流场的非平衡效应、表面的热辐射效应、催化效应和烧蚀效应以及热防护层内部的热传导效应，建立了初步的表面温度分布与气动热的耦合计算方法，完善了高超声速飞行器气动物理流场计算软件(AEROPH_Flow)。在表面材料为碳-碳(C-C)条件下，对飞行高度为65 km和飞行速度为8，10 km/s的半球以及飞行高度为50 km和飞行速度为8 km/s的球锥模型，开展了表面温度分布与气动热的耦合计算，验证了计算方法和计算软件，分析了表面温度分布对气动热环境的影响。研究结果表明：表面温度分布对气动热的计算结果有较大影响，在气动热环境的预测中，不仅要考虑热化学非平衡效应和表面催化效应的影响，还要考虑表面温度分布的影响，最好是采用表面温度分布与气动热耦合计算的方法，以减小表面温度分布对气动热计算结果的影响。为此,需要发展完善非平衡流场/表面催化和烧蚀/热传导温度场(气/表/固)的计算模型、耦合求解技术和计算软件，实现对高超声速飞行器的真实飞行条件下高温气体非平衡效应和气动热环境的精确模拟。

高超声速飞行器； 热化学非平衡效应； 表面温度分布； 气动热环境； 数值模拟

高超声速飞行器在临近空间飞行(或再入或滑翔)过程中，如果飞行马赫数在10以上，周围绕流流场的空气最高温度在4 000 K以上。在这样的条件下，周围绕流流场中高温空气中的氧分子O2和氮分子N2存在不同程度的离解和电离等化学反应，分子和原子的内能模式被不同程度的激发，这种化学-物理现象就是高温气体效应。高温气体效应包括化学效应和热力学效应，在高速流动情况下，化学效应分为化学冻结流、化学非平衡流和化学平衡流，热力学效应分为热力学冻结流、热力学非平衡流和热力学平衡流，即出现非平衡效应。同时，飞行器表面的防热材料会对高温空气的化学反应产生催化效应，使流场中高温化学反应生成的氧原子和氮原子在飞行器表面附近发生复合反应，而流场中的氧组分和氮组分也会使飞行器表面的防热材料发生氧化和氮化等热化学现象。若飞行马赫数非常高或外形特征尺寸非常小，气动加热严峻，飞行器表面的防热材料还会发生微烧蚀或烧蚀现象，烧蚀产生的气体组分又会与流场中的气体组分发生化学反应。以上所有这些气动-化学-物理现象，统称为高温气体非平衡效应。高温气体非平衡效应不仅对飞行器气动力特性和气动热环境产生严重影响，还会对飞行器目标特性和通信传输特性等产生影响。

近十年来，临近空间飞行器、飞船和深空探测器等高超声速飞行器的研制越来越受到高度关注，并得到了较大发展。特点是飞行器外形越来越复杂、再入或飞行速度越来越高、控制要求更精确和有效载荷更大等，在这种情况下，准确预测气动力特性和气动热环境就非常关键，高温气体非平衡效应的研究也就越来越得到重视。在工程应用中，一般是通过飞行验证试验、地面风洞试验和数值模拟与理论分析相结合的方法，解决高超声速飞行器飞行中的高温气体非平衡效应问题。目前，数值求解热化学非平衡Navier-Stokes方程的CFD数值方法是解决飞行器高温气体非平衡效应问题比较好的数值模拟与理论分析的方法，因为它能够自动模拟流场中任何区域的化学冻结/平衡/非平衡和热力学化学冻结/平衡/非平衡的流动特征。因此，需要建立精确的高温气体非平衡气动力特性和气动热环境CFD预测技术，以此弥补地面试验研究的局限性，为高超声速飞行器设计提供关键数据和技术支持，提升高超声速飞行器设计水平。

从20世纪90年代开始，关于高温气体非平衡效应和气动热环境的数值模拟与理论分析研究就有了较大发展，到目前为止，国内外发表了相当多的论文，建立了有关的计算软件。1990年前后，文献[1]～文献[3]对高温空气非平衡效应的流动控制方程和物理化学模型作了大量的研究和总结。1992年，Netterfield[4]对弹道靶试验的球模型开展了热化学非平衡流场的数值模拟研究，对比了计算和试验的激波位置，给出了驻点线的温度和组分分布，成为非平衡流场计算的典型对比算例。Muylaert等[5]针对飞行试验标模ELECTRE，考虑完全非催化和完全催化壁面条件，开展了风洞试验条件和飞行试验条件的热化学非平衡流场数值模拟，对比分析了计算和飞行试验测量的气动热数据以及风洞试验和飞行试验相似准则，成为非平衡气动热计算的典型对比算例。1993年开始，Keenan和Candler[6-7]针对表面防热材料为碳-碳(C-C)的半球和钝锥体模型，考虑流场的两温度热力学非平衡效应、表面催化和烧蚀效应、内部热传导效应，开展了轴对称情况下的非平衡流场与气动热环境耦合的数值模拟研究和理论分析。1998年，Hassan等[8]针对IRV-2的头部球锥段及其飞行弹道的大量特征点，考虑表面防热材料碳-碳的烧蚀边界变化，采用热化学非平衡气动热力学研究和分析软件SACCARA、材料热响应计算软件COYOTE和表面烧蚀组分计算软件ACE，并补充边界计算模型，开展了非平衡流场与气动热环境的耦合研究，给出了不同弹道的驻点线温度和组分分布、表面热流和表面温度分布。2005年，Edquist[9]基于有限催化特性和辐射平衡表面温度分布，采用LAURA计算软件，研究了火星探测器MSL的气动热环境。2007年，Hash等[10]采用DPLR、LAURA和US3D 3个非平衡流场计算软件，开展了FIRE II 的非平衡气动热环境计算和对比分析。Subrahmanyam和Papadopoulos[11]介绍了计算软件SPARTA在气动热环境分析方面的发展和应用情况，详细给出了计算软件的湍流模型、高温气体的物理化学模型、表面催化计算模型、无烧蚀情况的辐射平衡表面温度分布计算模型、烧蚀情况的一维热传导方程和表面能量守恒计算表面温度分布的计算模型以及材料特性数据库等。2009年前后，意大利航天研究中心(CIRA)的Pezzella等[12-13]考虑高温空气的化学非平衡和热力学非平衡效应、滑移效应和表面温度分布，表面温度取不同的等温度值或表面热辐射平衡条件确定的温度分布，通过数值求解非平衡Navier-Stokes方程的CFD方法(H3NS)，对比分析了热化学非平衡效应、表面催化效应和表面温度分布对飞行器气动力/热特性的影响。2012年，Surzhikov[14]采用多个振动温度非平衡流动控制方程和高温气体辐射传输方程耦合的非平衡辐射热环境计算软件NERAT和高温气体特性计算软件 ASTEROID，考虑辐射平衡表面温度分布，开展了ISS-CEV的气体辐射加热研究。Chen和Gokcen[15]利用非平衡流场计算软件DPLR和热响应与烧蚀计算软件TITAN，考虑表面的组分质量和能量守恒(表面温度分布)，研究了表面为碳基材料PICA的非平衡表面热化学效应。

从以上研究可以看出，在气动热环境数值研究方面，国外的发展趋势是考虑表面温度分布的影响，发展热化学非平衡流场/表面催化或烧蚀/热传导温度场(气/表/固)耦合的数值模拟技术。在国内，近十多年来，董维中[16-23]、李海燕[24]、柳军[25-26]、曾明[27-28]、樊菁[29]、马辉[30]、苗文博[31]、王逸斌[32-33]和范绪箕[34]等在高温气体非平衡效应和气动热环境的数值模拟与理论分析方面做了大量研究工作，取得了非常大的发展，研制了成熟的CFD高温气体非平衡流场计算软件。但表面计算模型还处于表面等温壁条件的情况，气动热环境与固体热传导温度场是独立求解的，考虑表面防热材料烧蚀情况的研究甚少，有些研究还是采用国外计算软件开展的。近年来，笔者开展了热化学非平衡模型和表面温度对气动热计算的影响研究[35]，发现在高马赫数和热化学非平衡条件下，表面温度对气动热计算结果有较大影响，气动热数值随着表面温度的变化规律变得非常复杂，不能再认为气动热遵从随着表面温度的升高而降低的规律，表面温度最好取接近真实飞行情况的分布。

要实现准确模拟高超声速飞行器的高温气体非平衡气动热环境的数值模拟技术，要解决的关键技术如下：①飞行器绕流流场的数值模拟，主要解决高温气体流动的非平衡效应和湍流效应问题，包括高温空气和烧蚀产物；②飞行器内部热传导温度场的数值模拟，主要解决表面防热材料热传导问题以及微烧蚀或烧蚀问题；③飞行器表面的计算模型，主要解决表面防热材料催化问题、飞行器绕流流场与飞行器内部热传导温度场之间的相互耦合问题，即气/表/固耦合的问题。

本文在研制高超声速飞行器气动物理流场计算软件(AEROPH_Flow)的基础上，围绕以上关键技术，基于流场的非平衡Navier-Stokes方程、表面的能量守恒方程和内部的热传导方程，考虑流场的非平衡效应、表面的热辐射效应、催化效应和烧蚀效应以及内部的热传导效应，发展表面温度分布与气动热的耦合计算方法，并对计算方法和计算软件进行验证，开展表面温度分布对气动热环境的影响分析。

1 非平衡绕流流场计算方法

1.1 控制方程和求解方法

控制方程是三维化学非平衡Navier-Stokes方程，其无量纲化形式为[16,19]

(1)

式中：Q为守恒变量；F、G和H为对流项：FV、GV和HV为黏性项；Re为雷诺数；W为非平衡源项。其中，

ρi(i=1, 2, …,N)为组分i的密度，N为所考虑的组分方程数；wi为组分i的化学反应源项；ρ为总密度；u、v和w分别为3个方向的速度；E为总内能。

采用结构网格的有限差分方法离散式(1)，对流项采用AUSMPW+格式离散，黏性项采用中心格式离散，隐式离散方程采用LU-SGS方法求解。为了克服方程的非平衡刚性问题，非平衡源项、对流项和黏性项采用全隐式处理。

1.2 高温气体化学反应模型

对于自由空气来流，O2和N2是主要组分。在高温气体流动中，就会有NO、N、O、NO+和e等。当表面材料为碳-碳时，就会发生烧蚀现象，在流场中可能有12种组分：O2、N2、NO、N、O、NO+、e、CO、CO2、CN、C和C2。表1给出了考虑表面碳-碳材料烧蚀效应中12种组分的31个高温空气/碳-碳烧蚀化学反应模型的具体反应式。

表1 高温空气/碳-碳烧蚀化学反应模型Table 1 Chemical reaction models of high-temperature air/C-C ablation

Notes: M1-M8are collider.

正逆反应速率常数为

式中：T为温度；i=1,2,…,Nr,Nr为高温气体化学反应个数；系数Af,i、Bf,i、Cf,i、Ab,i、Bb,i和Cb,i见参见文献[6]、文献[7]、文献[16]、文献[18]和文献[36]。第i个化学反应为

式中：ρj为组分j的密度；Mj为组分j的分子量。第j个组分的化学反应的生成源项为

(2)

1.3 热力学模型

热力学平衡模型是一个温度模型，考虑组分的束缚电子能量的激发效应和分子组分振动能量的激发效应。状态方程和能量关系式[16,19]为

表2 高温空气/碳-碳烧蚀组分物理化学数据[6,7,18]

1.4 表面烧蚀计算模型和催化条件

对于表面碳-碳材料情况，主要考虑由氧化反应和催化反应引起的热化学烧蚀[6,7]:

反应是不可逆的，前2个反应为氧化反应，第3个反应是由于表面催化引起的氧分子复合反应，这里没有考虑生成CO2的有关反应以及和氮相关的表面反应，因为它们的反应速率较低。反应的速率常数为

式中：Tw为飞行器表面温度；Mi为原子组分分子量；R为通用气体常数；αi为表面反应的催化效率，来源于试验拟合数据，具体形式为

由上述反应引起的组分质量通量为

(3)

在表面有质量引射的情况下，在表面处法向速度不为零，组分浓度与表面烧蚀机理有关。表面边界条件为

(4)

式中：uτ和vn分别为表面切向和法向速度；h为气体的焓；V为速度矢量；n为法向单位矢量，下标w表示参数在飞行器表面的当量值。速度无滑移和引射的条件为

(5)

表面组分质量分数ci,w满足表面组分质量守恒方程

(6)

式中：Di为组分扩散系数。表面无引射、表面完全非催化条件(NCW)为

(7)

表面无引射、表面完全催化条件(FCW)为

ci，w=ci，∞

(8)

式中：ci，∞为来流组分条件。

1.5 气动热计算公式

对于化学非平衡气体模型，表面热流由对流热流和组分扩散热流2部分组成，计算公式为

(9)

式中：κ和hi分别为热传导系数和i组分的焓值。

2 热传导温度场计算方法

在飞行器内部热传导温度场数值模拟方面，控制方程是无内热源的热传导方程，无量纲化控制方程形式[6,7,37]为

(10)

式中：

E为能量项；qx、qy和qz为3个方向的热通量；ρ、T、cp和κ分别为材料密度、温度、热容和热传导系数；L∞、T∞和t∞为长度、温度和时间的无量纲化参数；Rec为雷诺数无量纲量；ac∞为热扩散系数的无量纲量。在一般情况下，材料密度ρ可以认为是常数，材料热容cp和热传导系数κ是温度T的函数。

从飞行器表面到温度场的热流计算公式为

(11)

热传导温度场控制方程式(10)的时间离散采用隐式双时间步法,空间离散可以采用中心差分方法，通过表面计算模型、qw和qc与非平衡绕流流场控制方程式(1)耦合求解。对于材料传热有4类边界条件，第1类边界是给出表面瞬时温度分布，适用于耦合流场边界；第2类边界是给出表面热流分布，适用于耦合固体边界；第3类边界是给出对流放热情况，适用于不同固态介质的边界；第4类边界是固定温度，其他边界条件是绝热和对称等。由于流动时间尺度与传热时间尺度差异很大，为适应瞬时传热计算需求，可引入适当参数，实现松弛同步计算。在流场区域中，气-固传热边界由求解质量和能量守恒方程计算壁面参数，在固体区域，固-气传热边界采用第一类边界条件。

3 表面温度分布计算方法

对于高超声速飞行器气动热环境数值模拟，必须给定表面温度分布。最常规最简单的方法是等表面温度方法，任何表面位置的温度都是假设的一个恒温值，即等壁温方法,与实际的飞行情况差别较大；接近实际情况的方法是考虑表面热辐射效应、催化效应、烧蚀效应和热传导效应，建立表面能量守恒方程及其数值求解表面温度Tw的方法，并与流场控制方程耦合求解，流场稳定和收敛时，就可以得到表面温度分布。

利用准定常假设，由交界面一直延伸到物体深部的控制体的能量守恒方程[6,7,18]为

(12)

若不考虑热传导效应，只考虑表面热辐射效应、催化效应和烧蚀效应，表面能量守恒方程简化为

(13)

若不考虑烧蚀效应，只考虑表面热辐射效应、催化效应和热传导效应，表面能量守恒方程简化为

(14)

若不考虑烧蚀效应和热传导效应，只考虑表面热辐射效应和催化效应，表面能量守恒方程进一步简化为

(15)

表面能量守恒方程式(12)～式(15)是4种不同情况的表面温度分布计算方法的控制方程，可以应用于不同飞行器和不同飞行条件。在本文中，把求解不考虑热传导效应的表面能量守恒方程式(13)和式(15)的计算方法称为辐射平衡表面温度分布计算方法。

4 辐射平衡表面温度分布和气动热计算分析

为了验证辐射平衡表面温度分布计算方法，分析表面温度分布和防热材料烧蚀效应对气动热环境的影响，不考虑热传导效应，选择文献[6]的算例：计算模型为半径Rn=1.0 m的半球模型，飞行速度V=8 km/s和10 km/s，飞行高度H=65 km，飞行攻角α=0°。为了对比分析，考虑了如下5种情况。

1) 完全非催化等温壁计算模型：表面为完全非催化条件，表面温度Tw=1 500 K，流场高温气体化学反应模型为7个组分的空气化学反应模型，记为NCW_Tw=1 500 K。

2) 完全催化等温壁计算模型：表面为完全非催化条件，表面温度Tw=1 500 K，流场高温气体化学反应模型为7个组分的空气化学反应模型，记为FCW_Tw=1 500 K。

3) 完全非催化辐射平衡表面温度分布计算模型：表面为完全非催化条件，表面温度分布由只考虑表面材料热辐射效应的计算模型给出(式(15))，流场高温气体化学反应模型为7个组分的空气化学反应模型，记为NCW_Tw_Radiation。

4) 完全催化辐射平衡表面温度分布计算模型：表面为完全非催化条件，表面温度分布由只考虑表面材料热辐射效应的计算模型给出(式(15))，流场高温气体化学反应模型为7个组分的空气化学反应模型，记为FCW_Tw_Radiation。

5) 碳-碳防热材料辐射平衡表面温度分布计算模型：表面防热材料为碳-碳(C-C)，考虑碳-碳防热材料的实际催化特性和氧化特性，表面温度分布由考虑表面材料热辐射效应、实际催化特性和氧化特性的计算模型给出(式(13))，流场高温气体化学反应模型为12个组分(O2、 N2、NO、N、O、NO+、e、CO、CO2、CN、C、C2)的空气/碳-碳化学反应模型。

图1和图2给出了V=8 km/s的半球驻点线组分质量分数分布、表面温度和表面热流分布与文献[6]的比较。图1中横坐标r和X含义类似，其绝度值均为离开驻点的距离，纵坐标Cs为组分质量分数。图2中横坐标S为弧长。从图中可以看到，计算的主要高温气体组分分布、表面温度和热流分布与文献的结果一致，验证了本文建立的碳-碳防热材料表面氧化和催化特性计算模型以及辐射平衡表面温度分布计算方法在气动热环境数值模拟中的正确性。

图3给出了H=65 km、V=8 km/s，Rn=1 m时表面为完全非催化辐射平衡、完全催化辐射平衡和碳-碳防热材料辐射平衡3种情况下计算得到的表面热流和辐射平衡表面温度分布云图比较，从图中可以看出3种情况计算的表面热流和辐射平衡表面温度分布有较大差异。图4所示为V=8 km/s和10 km/s时5种情况下的表面热流和3种情况下的辐射平衡表面温度分布比较。从图中可以看到：①飞行速度从V=8 km/s到V=10 km/s，热流增加一倍左右，表面温度增加300～500 K；②不同表面催化特性计算模型得到的表面热流和温度分布有明显差异，完全非催化计算模型计算的表面热流(NCW_Tw=1 500 K和NCW_Tw_Radiation)和表面温度(NCW_Tw_Radiation)较低，完全催化计算模型计算的表面热流(FCW_Tw=1 500 K和FCW_Tw_Radiation)和表面温度(FCW_Tw_ Radiation)较高，而考虑实际的碳-碳防热材料有限催化特性和烧蚀效应的表面计算模型(C-C)计算的表面热流和温度在两者之间，在头部驻点区域，最高的完全催化条件的热流(FCW_Tw_ Radiation)是完全非催化条件的近3倍，即催化热流占总热流的65%左右，而考虑实际的碳-碳防热材料有限催化特性和烧蚀效应的总热流只是完全非催化条件的近2倍，即催化热流占总热流的50%左右，相应的表面温度相差150～600 K；③表面温度分布对表面热流分布有影响，特别是完全催化条件下的热流分布影响较大，而且飞行速度越高影响越大，辐射平衡表面温度分布下的驻点热流高于等温壁条件(Tw=1 500 K)下的热流，V=8 km/s时两者相差8%左右，V=10 km/s时两者相差15%左右。

图1 半球驻点线组分质量分数分布(V=8 km/s)Fig.1 Mass fraction distribution along stagnation line for semi-sphere case at V=8 km/s.

图2 半球表面热流和表面温度分布(V=8 km/s)Fig.2 Distribution of surface heat flux and surface temperature for semi-sphere at V=8 km/s

图3 3种情况下半球表面热流和表面温度云图(V=8 km/s)Fig.3 Contour of surface heat flux and surface temperature for semi-sphere at V=8 km/s with three conditions

计算结果分析表明：防热材料烧蚀效应或表面计算模型对气动热环境的预测有明显的影响，要准确预测气动热环境，必须考虑实际防热材料的表面辐射特性、表面催化特性和氧化等烧蚀特性，建立适合应用的非平衡流场与防热材料耦合的表面温度分布和气动热计算方法。

图4 不同飞行速度下半球表面热流和表面温度对比Fig.4 Comparison of surface heat flux and surface temperature of semi-sphere at different flight speed

5 考虑热传导效应的表面温度分布和气动热计算分析

为了验证考虑热传导效应的表面温度分布计算方法，选择了文献[7]的算例：计算模型为典型钝锥再入体，头部半径为Rn=0.127 m,飞行速度V=8 km/s，飞行高度H=50 km，飞行攻角α=0°；热力学模型为一温度模型，即不考虑热力学非平衡效应。考虑C-C烧蚀材料固体传热，传热参数采用的曲线拟合公式[7]为

公式的适用范围为300～3 000 K，如果温度超过3 000 K，取3 000 K时的数值，材料密度为1 775 kg/m3。流场和温度场计算网格如图5(a)所示，红色为流场部分网格，绿色为温度场网格。固体内表面尖点进行倒圆设计，倒圆半径为1 mm，见图5(b)。

图6所示为钝锥体流场和固体温度场的温度云图，流场中温度较高，最高温度达到15 000 K左右，而固体温度场温度相对较低，球头部分区域达到3 000 K以上。图7所示为钝锥表面热流、表面温度和轴线温度的分布，横坐标S为表面长度(以驻点为零点)，X为轴向坐标(以驻点为原点)。可以看出，与文献[7]结果符合良好，验证了考虑热传导效应的表面温度分布和气动热耦合计算方法。

图5 钝锥体流场和温度场的计算网格Fig.5 Flow field and temperature field computational mesh of sphere-cone

图6 钝锥体流场和固体温度场的温度云图Fig.6 Temperature contour of flow field and solid temperature field of sphere-cone

图7 钝锥体表面热流、表面温度和轴线温度分布Fig.7 Distribution of surface heat flux, surface temperature and ablator centerline temperature of sphere-cone

为了分析飞行器固体热传导对气动热环境预测的影响，对比计算了有无固体传热情况下烧蚀非平衡流场。图8给出了钝锥体有无热传导情况下表面热流和温度分布对比。可以看出：考虑固体热传导时，表面温度有较为明显的下降，较大区域温度下降幅值在100～200 K左右，但表面热流相差不大。其原因是气动热随表面温度变化的规律：对于流场计算来说，钝锥体温度场传热带走部分热量，使表面温度有一定程度下降；在头部驻点区域，表面温度较高，下降的相对幅度较小，因而对流场相对影响不大，热流变化较小；在锥身，虽然表面温度下降的相对幅度较大，但身部区域的表面温度值范围对热流影响相对较小，因此，表面热流变化同样较小。此外，值得指出的是，本文采用的计算条件，其固体内部采用的绝热壁条件，如果采用某种冷却的方式降温，弹体内部热传导对表面温度和热流的影响可能会更加明显。

图8 钝锥体有无热传导的表面热流和温度分布对比Fig.8 Comparison of surface heat flux and surface temperature distribution with and without heat conduction of sphere-cone

通过上述计算分析可以看出，在某些情况下，固体热传导对飞行器表面温度计算精度影响较为明显，要准确预测气动热环境，必须考虑飞行器内部固体结构的传热特性。

6 结束语

表面温度的确定方法越来越重要，它不仅严重影响气动力特性和气动热环境的预测和评估，还影响飞行器目标特性等问题的预测和评估。通过本文的研究工作，初步建立了表面温度分布与气动热的耦合计算方法，完善了高超声速飞行器高温气体非平衡流场计算软件，不仅为解决高马赫数和热化学非平衡条件下高超声速飞行器气动热环境精确预测问题奠定了坚实基础，还可应用于气动力特性和目标特性等方面的研究。针对不同飞行条件或飞行器外形或流场和表面的气动物理化学现象，可以采用简化的或完整的表面温度分布与气动热的耦合计算方法，模拟接近真实飞行条件的高温气体非平衡效应和气动热环境。

为了进一步完善表面温度分布与气动热的耦合计算方法，还需要开展大量研究工作：

1) 开展飞行器表面的物理化学计算模型和内部热响应计算模型研究，建立完善的非平衡流场/表面催化和烧蚀/热传导温度场(气/表/固)的计算模型、耦合求解技术和计算软件，实现模拟高超声速飞行器真实飞行的高温气体非平衡效应和气动热环境。

2) 开展高温非平衡气体物理化学模型研究。国内没有建立完整的大气空气参数和高温气体反应动力学模型及其物理化学参数，现有数据基本上是引用国外的数据，更没有防热材料烧蚀产物、控制和冷却喷流气体与空气的高温化学反应模型及其物理化学参数。热力学非平衡的参数更加复杂，分子组分的振动松弛特征时间和特征温度等物理化学参数缺乏，特别是防热材料烧蚀产物的分子组分的数据。另外，各种高温气体组分输运系数也较缺乏，在飞行器第二宇宙速度再入时，周围高温气体几乎是高电子密度的等离子体，即使是纯空气不考虑表面烧蚀，也存在输运系数适用性的问题。

3) 开展高温气体非平衡流场数值计算软件的试验验证研究。由于试验设备和测量技术的限制，现有计算分析手段还没有得到充分的验证，需要发展高焓试验设备和测量技术，开展飞行试验验证研究。

4) 开展缝隙和干扰区气动热环境数值模拟与计算网格生成技术研究。飞行器外形越来越复杂，既要捕捉流动的化学反应和非平衡特性，又要精确模拟缝隙和干扰区气动热环境，对计算方法和计算网格的要求更加严格和精细。计算网格不仅有绕流流场的计算网格，还有热响应温度场的计算网格和烧蚀动表面的动态计算网格。

5) 开展高温气体非平衡、烧蚀和热响应非定常计算方法研究。国内开展比较少，需要重视这方面的研究。

随着国家高超声速飞行器技术和计算机技术的发展，通过大力发展高超声速飞行器飞行试验验证技术、地面风洞试验技术和数值模拟技术，解决高超声速飞行器的高温气体非平衡效应和气动热环境耦合问题的能力会得到进一步提升。

[1] Cander G V, Maccormack R W. The computation of hypersonic ionized flows in chemical and thermal nonequilibrium, AIAA-1988-0511[R]. Reston: AIAA, 1988.

[2] Gupta R N, Yos J M, Thompson, et al. A review of reaction rates and thermodynamic and transport properties for 11-species air model for chemical and thermal nonequilibrium calculation to 30 000 K, NASA-TM-101528[R]. Washington, D.C.: NASA, 1989.

[3] Park C. Nonequilibrium hypersonic aerothermodynamics [M].New York: John Wiley &Sons, 1990: 145-218.

[4] Netterfield M P. Validation of a Navier-Stokes code for thermochemical nonequilibrium flows, AIAA-1992-2878[R]. Reston: AIAA, 1992.

[5] Muylaert J, Walpot L, Hauser J, et al. Standard model testing in the European high facility F4 and extrapolation to flight, AIAA-1992-3905[R]. Reston: AIAA, 1992.

[6] Keenan J A , Candler G V. Simulation of ablation in earth atmospheric entry, AIAA-1993-2789[R]. Reston: AIAA, 1993.

[7] Keenan J A, Candler G V. Simulation of graphite sublimation and oxidation under re-entry conditions, AIAA-1994-2083[R]. Reston: AIAA, 1994.

[8] Hassan B, Kuntz D W, Potter D L, et al. Coupled fluid/thermal prediction of ablating hypersonic vehicles, AIAA-1998-0168[R]. Reston: AIAA, 1998.

[9] Edquist K T. Afterbody heating predictions for a Mars science laboratory entry vehicle, AIAA-2005-4817[R]. Reston: AIAA, 2005.

[10] Hash D, Olejniczak J, Wright M J, et al. FIRE II calculations for hypersonic nonequilibrium aerothermodynamics code verification: DPLR, LAURA, and US3D, AIAA-2007-0605[R]. Reston: AIAA, 2007.

[11] Subrahmanyam P, Papadopoulos P. Development of a parallel CFD solver SPARTA for aerothermodynamic analysis, AIAA-2007-2976[R]. Reston: AIAA, 2007.

[12] Pezzella G, Votta R . Finite rate chemistry effects on the high altitude aerodynamics of an Apollo-shaped reentry capsule, AIAA-2009-7306[R]. Reston: AIAA, 2009.

[13] Votta R, Schettino A, Ranuzzi G, et al. Hypersonic low-density aerothermodynamics of Orion-like exploration vehicle[J]. Journal of Spacecraft and Rockets, 2009, 46(4): 781-787.

[14] Surzhikov S. The effect of non-equilibrium dissociation on radiative heating of entering space vehicle, AIAA-2012-1146[R]. Reston: AIAA, 2012.

[15] Chen Y K, Gokcen T. Effect of non-equilibrium surface thermochemistry in simulation of carbon based ablators[J]. Journal of Spacecraft and Rockets, 2013, 50(5): 917-926.

[16] Dong W Z. Numerical simulation and analysis of thermo-chemical non-equilibrium effects at hypersonic flows[D].Beijing: Beijing University of Aeronautics and Astronautics, 1996 (in Chinese). 董维中. 热化学非平衡效应对高超声速流动影响的数值计算与分析[D]. 北京： 北京航空航天大学， 1996.

[17] Dong W Z, Le J L, Liu W X. The determination of catalytic rate constant of surface materials of testing model in the shock tube[J]. Experiments and Measurements in Fluid Mechanics, 2000, 14(3): 1-6 (in Chinese). 董维中， 乐嘉陵， 刘伟雄. 驻点壁面催化速率常数确定的研究[J]. 流体力学实验与测量, 2000, 14(3): 1-6.

[18] Dong W Z, Gao T S, Zhang Q Y. Numerical simulation of 3-D graphite ablation flow fields over a hypersonic reentry blunt body[J]. Acta Aerodynamica Sinica, 2001, 19(4): 388-394 (in Chinese). 董维中， 高铁锁， 张巧云. 高超声速三维碳-碳烧蚀流场的数值研究[J]. 空气动力学学报, 2001, 19(4): 388-394.

[19] Dong W Z. Thermal and chemical model effect on the calculation of aerodynamic parameter for hypersonic reentry blunt body[J]. Acta Aerodynamica Sinica, 2001, 19(2): 197-202 (in Chinese). 董维中. 气体模型对高超声速再入钝体气动参数计算影响的研究[J]. 空气动力学学报， 2001， 19(2): 197-202.

[20] Dong W Z, Le J L, Gao T S. Numerical analysis for correlation of standard model testing in high enthalpy facility and flight test[J]. Experiments and Measurements in Fluid Mechanics, 2002, 16(2): 1-8 (in Chinese). 董维中， 乐嘉陵， 高铁锁. 钝体标模高焓风洞试验和飞行试验相关性的数值分析[J]. 流体力学实验与测量， 2002, 16(2): 1-8.

[21] Dong W Z, Gao T S, Ding M S. Numerical studies of the multiple vibrational temperature model in hypersonic non-equilibrium flows[J]. Acta Aerodynamica Sinica, 2007, 25(1): 1-6 (in Chinese). 董维中， 高铁锁， 丁明松. 高超声速非平衡流场多个振动温度模型的数值研究[J]. 空气动力学学报, 2007, 25(1): 1-6.

[22] Dong W Z, Gao T S, Zhang Z C. Numerical analysis of thermal and chemical nonequilibrium testing flow-field around the reentry blunt body[J]. Acta Aerodynamica Sinica, 2008, 26(2): 163-166 (in Chinese). 董维中， 高铁锁， 张志成. 再入体实验模型热化学非平衡绕流流场的数值分析[J]. 空气动力学学报, 2008, 26(2): 163-166.

[23] Dong W Z, Gao T S, Ding M S, et al. Numerical analysis for the effect of Silicon based material ablation on the flow field around reentry blunt body[J]. Acta Aerodynamic Sinica, 2010, 28(6): 708-714 (in Chinese). 董维中， 高铁锁， 丁明松，等. 硅基材料烧蚀产物对再入体流场特性影响的数值计算[J]. 空气动力学学报， 2010, 28(6): 708-714.

[24] Li H Y. Numerical simulation and analysis of high temperature real gas effects at hypersonic aerothermodynamics[D]. Mianyang: China Aerodynamics Research and Development Center, 2007 (in Chinese). 李海燕. 高超声速高温气体流场的数值模拟[D]. 绵阳：中国空气动力研究与发展中心， 2007.

[25] Liu J, Liu W, Zeng M, et al. Numerical simulation of 3D hypersonic thermochemical nonequilibrium flow[J]. Acta Mechanica Sinica, 2003, 35(6): 730-735 (in Chinese). 柳军， 刘伟， 曾明, 等. 高超声速三维热化学非平衡流场的数值模拟[J]. 力学学报, 2003, 35(6): 730-735.

[26] Liu J. Experimental and numerical research on thermo- chemical nonequilibrium flow with radiation phenomenon [D]. Changsha: National University of Defense Technology, 2004 (in Chinese). 柳军. 热化学非平衡流及其辐射现象的实验和数值计算研究[D]. 长沙： 国防科学技术大学， 2004.

[27] Zeng M, Hang J, Lin Z B， et al. Numerical analysis of the effects of different thermo-chemical nonequilibrium models on hypersonic nozzle flow[J].Acta Aerodynamica Sinica, 2008, 24(3): 346-349 (in Chinese). 曾明， 杭建， 林贞彬, 等. 不同热化学非平衡模型对高超声速喷管流场影响的数值分析[J]. 空气动力学学报, 2008, 24(3): 346-349.

[28] Gao B, Hang J, Lin Z B, et al. The experiment exploration of catalyst effects on aerodynamic heat in real gas effects[J].Experiments and Measurements in Fluid Mechanics, 2004, 18(2): 55-58 (in Chinese). 高冰, 杭建, 林贞彬, 等. 高温真实气体效应中催化效应对气动热影响的实验探索[J]. 流体力学实验与测量, 2004, 18(2): 55-58.

[29] Fan J. Criteria on high-temperature gas effects around hypersonic vehicles[J]. Chinese Journal of Theoretical and Applied Mechanics， 2010, 42(4): 591-596 (in Chinese). 樊菁. 高超声速高温气体效应判据[J]. 力学学报, 2010, 42(4): 591-596.

[30] Ma H , Zhao L, Wang F M. Numerical simulations of hypersonic three-dimensional chemical nonequilibrium flows[J]. Journal of Beijing University of Aeronautics and Astronautics, 2004, 30(2): 168-172 (in Chinese). 马辉, 赵烈, 王发民. 高超声速三维化学非平衡流动的数值模拟[J]. 北京航空航天大学学报, 2004, 30(2): 168-172.

[31] Miao W B, Cheng X L, Ai B C. The influence of catalyze condition on thermal environment predicting[J]. Spacecraft Environment Engineering, 2009, 26(Sup.): 45-49 (in Chinese). 苗文博, 程晓丽, 艾邦成. 壁面催化条件对热环境预测的影响[J]. 航天器环境工程, 2009, 26(增刊): 45-49.

[32] Wang Y B, Wu Y Z, Liu X Q. Simulation of 3D hypersonic ionized and radiating flows in thermal and chemical nonequilibrium[J]. Chinese Journal of Computational Physics, 2008, 25(4): 421-426 (in Chinese). 王逸斌, 伍贻兆, 刘学强. 耦合辐射的三维热化学非平衡流数值模拟[J]. 计算物理, 2008, 25(4): 421-426.

[33] Zhang X H, Wu Y Z, Wang J F. Aero-heating numerical simulation of axisymmetric reentry vehicle body[J]. Chinese Journal of Applied Mechanics, 2012, 29(3): 284-290 (in Chinese). 张向洪, 伍贻兆, 王江峰. 轴对称再入舱模型气动热特性数值模拟研究[J]. 应用力学学报, 2012, 29(3): 284-290.

[34] Fan X J. On the mathematical simulation of non-equilibrium hypersonic around the flying vehicle[J]. Advances in Mechanics, 2004, 34(2): 224-236 (in Chinese). 范绪箕. 关于飞行器高超声速不平衡气体绕流的数值模拟[J]. 力学进展, 2004, 34(2): 224-236.

[35] Dong W Z, Ding M S, Gao T S， et al. The influence of thermo-chemical nonequilibrium model and surface temperature on heat transfer rate[J]. Acta Aerodynamica Sinica, 2013， 31(6)： 692-698 (in Chinese). 董维中， 丁明松， 高铁锁， 等. 热化学非平衡模型和表面温度对气动热计算影响分析[J]. 空气动力学学报, 2013, 31(6): 692-698.

[36] Blottner F G. Prediction of electron density in the boundary layer on entry vehicle with ablation, N71-21113[R]. 1971.

[37] Tao W Q. Numerical heat transfer[M]. 2nd ed. Xi’an: Xi’an Jiao Tong University Press, 2008: 78-90 (in Chinese). 陶文铨. 数值传热学[M]. 2版. 西安： 西安交通大学出版社， 2008: 78-90.

Tel: 0816-2463298

E-mail: dongwz1966@163.com

*Corresponding author. Tel.： 0816-2463298 E-mail: dongwz1966@163.com

Numerical study of coupled surface temperature distribution and aerodynamic heat for hypersonic vehicles

DONG Weizhong*, GAO Tiesuo， DING Mingsong, JIANG Tao, LIU Qingzong

ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

In order to predict precisely aeroheating and effects of high-temperature non-equilibrium gas during design the thermal protection system of hypersonic vehicles, the simulation methods of coupled surface temperature and heat transfer rate are developed based on the Navier-Stokes equations of non-equilibrium flow fields, the energy conservation equation at the surface with radiation, catalytic action and ablation, and the unsteady heat conduction equations of heat-shield, and the computational code AEROPH_Flow is perfected ,which are developed by us for numerically simulating the flow over hypersonic vehicles and predicting the aero-physical characteristics . The numerical simulation results are presented, including a semi-sphere geometries at the altitude of 65 km with the free stream velocity of 8 km/s and 10 km/s, a sphere-cone geometries at the altitude of 50 km with 8 km/s, and the polycrystalline graphite is selected as ablative material on the surface. The distributions of surface temperature and heat transfer rate are obtained and the analysis is done for the influence of the surface temperature distribution on heat transfer rate. The results show that the surface temperature distribution has a more important influence on the computational results of heat transfer rate, the factors considered in the high-precision prediction of aero-thermal environment are not only thermo-chemical non-equilibrium effect and surface catalytic effect, but also the surface temperature distributions, so the best method for high-precision prediction of aero-thermal environment is the coupling of surface temperature and heat transfer rate, and it is essential to develop perfect physical models, solving methods and numerical simulation codes of coupled non-equilibrium flow field, surface catalytic action and ablation, and heat conduction of heat-shield for the high-precision prediction of aero-thermal environment of hypersonic vehicles under the real flight condition.

hypersonic vehicle; thermo-chemical non-equilibrium effect; surface temperature distribution; aerodynamic heat environment; numerical simulation

2014-07-04； Revised： 2014-09-10； Accepted： 2014-10-13； Published online： 2014-10-14 14:10

National Natural Science Foundation of China (91216204)

2014-07-04; 退修日期： 2014-09-10; 录用日期： 2014-10-13; 网络出版时间： 2014-10-14 14：10

www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0241.html

国家自然科学基金 (91216204)

Dong W Z, Gao T S, Ding M S, et al. Numerical study of coupled surface temperature distribution and aerodynamic heat for hypersonic vehicles[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 311-324. 董维中， 高铁锁， 丁明松， 等．高超声速飞行器表面温度分布与气动热耦合数值研究[J]． 航空学报, 2015, 36(1): 311-324.

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2014.0241

V434+.11

A

1000-6893(2015)01-0311-14

董维中 男， 博士， 研究员， 博士生导师。主要研究方向： 气动物理学和计算高温气体动力学。

*通讯作者.Tel.： 0816-2463298 E-mail: dongwz1966@163.com

URL： www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0241.html