热化学非平衡模型和表面温度对气动热计算影响分析

董维中，丁明松，高铁锁，江 涛

（中国空气动力研究与发展中心 计算空气动力研究所，四川 绵阳 621000）

0 引 言

高超声速飞行器在大气层飞行过程中，不仅存在高热流的气动热环境问题，还存在严重的高温气体效应／非平衡效应问题，这就给高超声速飞行器热防护设计和气动设计增加非常大的困难。在这种情况下，准确预测高超声速飞行器气动力特性和表面气动热环境就非常关键。

对高温气体效应／非平衡效应问题研究越来越得到重视，国内外建立了有关的计算软件和发表了相当多论 文［1－8］。2009 年 前 后，意 大 利 航 天 研 究 中 心（CIRA）的 Votta R等人［7－8］，考虑高温空气的化学非平衡和热力学非平衡效应、滑移效应和表面温度分布，表面温度取不同的等温度值或表面热辐射平衡条件确定的温度分布，通过数值求解非平衡Navier－Stokes方程的CFD方法（H3NS），对比分析了热化学非平衡效应、表面催化效应和表面温度分布对飞行器气动力／热特性的影响。

在数值计算气动热中，影响准确预测高超声速飞行器气动热环境的因素非常多，不仅有计算网格和计算格式及其边界处理方法，还有气体模型、表面温度和表面催化特性等等。在气体模型中，不仅有完全气体模型和化学反应混合气体模型，还有气体热力学模型。在化学反应混合气体模型中，常用的有5组分、7组分和11组分空气化学反应模型。在这些影响准确预测高超声速飞行器气动热环境因素中，目前国内外对空气化学反应的组分模型、热力学非平衡模型和表面温度对气动热的计算结果的影响规律还不是非常清楚，在高马赫数和热化学非平衡条件下，气动热数值随表面温度变化有什么特征，这些问题都需要深入分析研究。

本文考虑高温空气化学反应效应和热力学非平衡效应，利用自主开发的气动物理计算软件系统（AEROPH）中的高超声速飞行器高温气体／非平衡效应流场计算软件（AEROPH＿Flow），开展热化学非平衡气动热数值计算，分析高温空气化学反应的组分模型、热力学非平衡模型和表面催化特性对气动热的计算结果的影响，考察在高马赫数和热化学非平衡条件下气动热数值随表面温度的变化特征，提升准确预测高超声速飞行器气动热环境的能力。

1 非平衡流场控制方程和计算方法

控制方程是三维热化学非平衡Navier－Stokes方程，无量纲化形式如下［9－10］：

其中，

这里i＝（1，2，…，N），N是所考虑的组分方程数。ρi是组分i的密度，EV是分子组分的总振动能量，Re是Reynolds数，wi和wV是化学非平衡源项和振动非平衡能量源项。

对于热化学非平衡流动情况，状态方程和能量关系式如下［9－10］：

这里ei，tr和ei，e是组分i的平动和转动能、束缚电子的激发能。θVj是组分j的振动特征温度，g0i和g1i是组分i的束缚电子第0和第1能级的简并因子，θei是第一能级的特征温度。Δ是组分生成能。空气组分的物理化学数据见文献［9－10］。

对于化学非平衡流动情况，在控制方程中去掉振动非平衡能量方程和相关项，状态方程和能量关系式如下［9－10］：

完全气体的控制方程和状态方程见文献［9］。方程（1）采用全隐式的对称型TVD格式进行差分离散，粘性项用中心差分格式离散。

2 空气化学反应模型

较常用的空气化学反应模型有5组分（O2，N2，NO，O，N）、7组分（O2，N2，NO，O，N，NO＋，e）和11组分（O2，N2，NO，O，N，NO＋，e，O＋2，N＋2，O＋，N＋）的模型。在本文研究中，选用了5组分、7组分和11组分的Dunn－Kang空气化学反应模型。

3 热力学非平衡模型

对于热力学非平衡，只考虑振动非平衡效应，考虑振动能量非平衡的组分有O2，N2，NO，NO＋，O＋2和等分子组分。为了减少控制方程数量，振动非平衡效应考虑一个振动温度。振动非平衡能量源项有如下形式：

这里、EVj和τVj分别是分子组分j的平衡振动能、非平衡振动能和振动能松弛的特征时间，τVj具体计算公式计算见文献［9－10］。

4 气动热计算方法

对于完全气体模型，表面热流计算公式为：

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

对于热化学非平衡气体模型，表面热流由平动项的对流热流、振动项对流热流和组分扩散热流三部分组成，计算公式为：

这里k和kV是平动项和振动项热传导系数。hi和Di是组分i的焓和扩散系数。对于完全非催化条件，有＝0；对于完全催化条件，有ciw＝ci∞（来流条件）。

5 空气化学反应模型和热力学模型对气动热计算影响的分析

针对高超声速飞行器，我们建立了气动物理计算软件系统，并开展了大量的验证与确认工作，在高温气体／非平衡效应流场计算软件方面，开展了大量的与国外文献［11－13］以及国内外的飞行试验、地面试验和理论计算的结果对比分析工作［9－10，14－17］。

计算模型：球头模型，RN＝35mm；计算状态：H＝55km，M＝20；热力学非平衡模型：一温度模型（T）和两温度模型（T－Tv）；空气化学反应模型：5组分、7组分和11组分化学反应模型；表面条件：Tw＝300K，完全催化（FCW）和完全非催化（NCW）。

图1给出了不同组分空气化学反应模型的热流分布，从图中可以看到，7组分模型和11组分模型计算的热流分布几乎没有差别，5组分模型计算的热流值偏小，在M＝20和完全催化条件下热流偏小13%左右，原因可能是5组分模型没有考虑电离效应。图2给出了不同热力学模型的热流分布，从图中可以看到，不同热力学模型计算的热流分布有一定差异，最大差异在10%左右，两温度模型计算的热流值高于一温度模型计算的热流值。从总的来看，在同样催化条件下，不同组分空气化学反应模型和不同热力学模型计算的热流分布差异主要集中在球头模型的驻点区域，后面差别较小。

6 表面温度对气动热计算影响的分析

计算模型：球锥模型，RN＝35mm，θ＝5°，L＝2m；计算状态：H＝60km和55km，M＝20；热力学非平衡模型是一温度模型，空气化学反应模型是7组分模型；表面条件：Tw＝300K，800K，1500K，2000K，完全催化和完全非催化。

图1 不同组分空气化学反应模型的热流分布（NCW和FCW）Fig.1 Distribution of heat transfer rate on the hemisphere for different species chemical reaction models（NCW and FCW）

图2 不同热力学模型的热流分布（NCW和FCW）Fig.2 Distribution of heat transfer rate on the hemisphere for different thermodynamic nonequilibrium models（NCW and FCW）

图3 球锥模型H＝55km的压力和表面热流分布云图Fig.3 Contours of pressure and heat transfer rate for the sphere－cone body at the altitude of 55km

图3给出Rn＝3 5mm了球锥模型H＝5 5km、M∞＝18、Tw＝300K、Ns＝7、NCW 状态热流分布云图，从图可以看出：高热流区域和空气化学反应最强区域在球锥体的头部区域，在后身部区域热流低和空气化学反应弱。

图4给出了球锥模型驻点热流和球锥末端热流及其成分随表面温度变化的分布。从图可以看出：（1）不同表面温度计算的热流分布有一定的差异，在Tw＝300K～800K之间，表面温度越高，表面热流越高，在Tw＝800K～2000K之间，随表面温度增加，表面热流先缓慢升高，Tw＝1500K之后，有降低趋势；（2）表面温度增加，组分扩散项热流一直增加，表面温度从300K到2000K，组分扩散项热流在总热流中的百分比从33%左右到40%左右；（3）表面温度变化对球锥头部区域热流值影响大，但对球锥身部区域影响小；（4）表面催化效应对球锥体的头部区域热流影响大，对后身部区域热流影响小。

图4 球锥模型驻点和末端热流随表面温度变化的分布Fig.4 Variation of the heat transfer rate at the stagnational point and the end of the sphere－cone body for different surface temperatures

从分析可以看出：气动热随着表面温度的变化规律在高马赫数和非平衡条件下就变得非常复杂，不能再认为气动热遵从随着表面温度的升高而降低的规律。

7 高升阻比外形不同表面温度气动热计算分析

图5 高升阻比外形驻点热流随飞行高度变化分布Fig.5 Variation of the heat transfer rate at the stagnational point of the high lift－to－drag hypersonic vehicle for different flight altitudes and surface temperatures

对于典型简化高升阻比外形，头部是RN＝35mm的半球，计算状态：H＝45km～70km，M＝11～23，图5给出了表面温度Tw＝600K和Tw＝1000K时高升阻比外形驻点热流随飞行高度变化分布。由图可以看出：（1）从H＝70km到H＝45km，驻点热流最大值在H＝55km；（2）表面材料催化效应对热流值影响较大，不同表面温度热流值有所不同，当表面温度Tw＝600K时，对于完全非催化条件，驻点热流在1650kW／m2至3000kW／m2之间，而对于完全催化条件，驻点热流在1900kW／m2至4600kW／m2之间，当表面温度Tw＝1000K时，对于完全非催化条件，驻点热流在1500kW／m2至3200kW／m2之间，而对于完全催化条件，驻点热流在1800kW／m2至5000kW／m2之间；（3）完全气体（PG）模型计算的热流在完全非催化和完全催化热流之间；（4）表面温度对热流计算存在一定影响，表面温度由600K上升到1000K时，H＝45km时驻点热流略微下降，而其他高度热流都有升高。

图6给出了H＝55km不同表面温度和不同飞行马赫数高升阻比外形热流变化分布。由图可以看到：表面温度变化对头部区域的热流影响比较明显，当M＝10时，热流值随着表面温度的增加而下降；当M＝14时，在Tw＝300K～600K时，热流值随着表面温度的增加而增加，在Tw＝600K～1500K时，热流值随着表面温度增加变化缓慢，且有下降趋势；当M＝18时，在Tw＝300K～900K时，热流值随着表面温度的增加而增加，在Tw＝900K～1500K时，热流随着表面温度增加变化缓慢，且有下降趋势。从分析结果看，高马赫数和热化学非平衡条件下，在计算气动热时，表面温度最好取保守上限值，这样就可以找到最大的驻点热流值。

图6 高升阻比外形不同表面温度和飞行马赫数条件下驻点热流变化分布（FCW）Fig.6 Variation of the heat transfer rate at the stagnational point of the high lift－to－drag hypersonic vehicle for different flight Mach number and surface temperature

8 结束语

通过本文的研究工作，可得到如下结论：

（1）对高温气体效应／非平衡效应严重的区域，为了准确预测高超声速飞行器气动热环境，我们要根据飞行环境的热化学机制或空气化学反应和非平衡效应的强弱，选择适当的空气化学反应模型和热力学模型，不能随意用完全气体模型、少组分的空气化学反应模型以及不考虑热力学非平衡效应，特别是飞行速度非常高的情况。在本文的计算条件下，数值分析表明：（a）空气化学反应模型7组分模型和11组分模型计算的热流分布几乎没有差别，而5组分模型计算的热流值偏小达10%以上；（b）不同热力学模型计算的热流分布差异较小，最大差异在10%左右，两温度模型计算的热流值高于一温度模型计算的热流值。

（2）高马赫数和热化学非平衡条件下，气动热的数值随着表面温度的变化规律变得非常复杂，不能再认为气动热的数值遵从随着表面温度的升高而降低的规律；在计算气动热时，表面温度最好取接近真实飞行情况分布和不同的固定值，这样就可以找到最大的或准确的热流值，在高超声速飞行器的热防护设计中更有应用价值。

本文只是分析了在再入或飞行速度在7km／s以下的条件，还需要开展系统的分析研究，如更高再入或飞行速度、不同模型尺寸和变表面温度等情况。

［1］GNOFFO P A.Applications of program LAURA to threedimensional AOTV flowfields［R］.AIAA－86－0565.

［2］HARTUNG L C.Development of a nonequilibrium radiative heating prediction method for coupled flowfield solutions［R］.AIAA 91－1406.

［3］KEENAN J A，CANDLER G V.Simulation of graphite ablation and oxidation under re－entry conditions［R］.AIAA－94－2083.

［4］KEENAN J A，CANDLER G V.Simulation of ablation in earth atmospheric entry［R］.AIAA 93－2789，1993.

［5］曾明，杭建，林贞彬，等.不同热化学非平衡模型对高超声速喷管流场影响的数值分析［J］.空气动力学学报，2008，24（3）：346－349.

［6］柳军，刘伟，曾明，等.高超声速三维热化学非平衡流场的数值模拟［J］.力学学报，2003，35（6）：730－735.

［7］PEZZELLA G，VOTTA R.Finite rate chemistry effects on the high altitude aerodynamics of an apollo－shaped reentry capsule［R］.AIAA 2009－7306.

［8］VOTTA R，et al.Hypersonic low－density aerothermodynamics of orion－like exploration vehicle［J］.Journal ofSpacecraftandRockets，2009，46（4）：781－787.

［9］董维中.热化学非平衡效应对高超声速流动影响的数值计算与分析［D］.［博士学位论文］.北京航空航天大学，1996.

［10］董维中.气体模型对高超声速再入钝体气动参数计算影响的研究［J］.空气动力学学报，2001，19（2）：197－202.

［11］NETTERFIELD M P.Validation of a Navier－Stokes code for thermochemical nonequilibrium flows［R］.AIAA 92－2878.

［12］MUYLAERT J，et al.Standard model testing in the European high facility F4and extrapolation to flight［R］.AIAA－92－3905.

［13］CANDER G V，MACCORMACK R W.The computation of hypersonic ionized flows in chemical and thermal nonequilibium［R］.AIAA－88－0511.

［14］董维中，高铁锁，丁明松.高超声速非平衡流场多个振动温度模型的数值研究［J］.空气动力学学报，2007，25（1）：1－6.

［15］董维中，高铁锁，张志成.再入体实验模型热化学非平衡绕流流场的数值分析［J］.空气动力学学报，2008，26（2）：163－166.

［16］董维中，乐嘉陵，刘伟雄.驻点壁面催化速率常数确定的研究［J］.流体力学实验与测量，2000，14（3）：1－6.

［17］董维中，乐嘉陵，高铁锁.钝体标模高焓风洞试验和飞行试验相关性的数值分析［J］.流体力学实验与测量，2002，16（2）：1－8.