0 引言

飞行器在大气中以高超声速飞行时，激波层内气体受到强激波的压缩作用而急剧升温，高温气体通过对流及辐射过程对飞行器进行加热，即所谓的气动加热。而在高温环境下，气体分子会振动激发、离解甚至电离，此时气体特性不再遵循完全气体条件，即所谓的高温真实气体效应^[1]。由于振动激发和分解的过程会消耗大量气流热量，可以显著降低波后的气体温度，进而缓解了气动加热。但是在接近飞行器表面附近的气流薄层内，原子、离子成分会受到飞行器表面材料化学特性作用而发生一定程度的复合反应，并释放大量热量，这种壁面效应称之为壁面催化效应，可以改变壁面附近的气体组分及能量分布状态^[2]。由浓度梯度引发的扩散效应会促进外层原子向壁面附近扩散释热，加剧气动加热。从20世纪50年代开始，关于壁面催化效应对气动热的影响已经成为重要的研究课题^[3-4]，同时飞行试验及地面风洞试验结果也不断验证了壁面催化效应的重要影响作用^[5-6]。

由于飞行试验及地面风洞试验成本高、测试难度大，能获得的有效试验数据比较匮乏，数值模拟是高超声速流动基础研究中不可或缺的工具^[7-10]。在对壁面催化效应的数值模拟中，目前有非催化、有限速率催化、完全催化和超催化等壁面条件^[11]。其中，非催化和完全催化壁面条件是两种极限情况：非催化壁面条件(Non-Catalytic Wall, NCW)下，壁面处原子、离子不发生复合反应，其组分百分比例不发生改变；而完全催化壁面条件(Fully Catalytic Wall, FCW)下，壁面处原子、离子完全复合成分子结构，即复合反应速率假设为无穷大，此时释放的热量远大于非催化壁面条件。实际飞行中，壁面材料特性应介于两种极限条件之间，应用何种壁面催化边界条件尚没有判断准则^[12-13]。而工程上对于飞行器的热环境预测，也常常采用完全催化壁面条件这种极不利的情况来保证热防护系统的可靠性，是一种比较保守做法，一定程度上会加大运行成本，因此对于不同飞行条件下壁面催化条件对于热流的影响亟需细致的机理和规律性研究。本文主要通过数值模拟方法来研究非催化及完全催化两种壁面催化条件对高超声速典型模型气动热的影响规律，为高超声速飞行器气动热预测提供数据支持。

1 数值方法

数值模拟采用了多组分N-S方程，在直角坐标系下，其轴对称控制方程如下：

$ \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial {\bf{F}}}}{{\partial x}} + \frac{{\partial {\bf{G}}}}{{\partial y}} + \delta \frac{{\bf{H}}}{y} = \frac{{\partial {{\bf{F}}_v}}}{{\partial x}} + \frac{{\partial {{\bf{G}}_v}}}{{\partial y}} + \delta \frac{{\partial {{\bf{H}}_v}}}{y} + {\bf{S}} $

(1)

其中，Q为位置变量，F、G和F_ν、G_ν为x、y方向的无黏、有黏通量，H、H_ν为无黏、有黏轴对称几何源项，S为化学反应源项，具体表达如下：

$ \begin{array}{l} \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} {{\rho _1}}\\ {{\rho _2}}\\ \vdots \\ {{\rho _{ns}}}\\ {\rho u}\\ {\rho v}\\ E\\ {\rho {e_{vib}}} \end{array}} \right], \;\;\;{\kern 1pt} {\bf{F}} = \left[ {\begin{array}{*{20}{c}} {{\rho _1}u}\\ {{\rho _2}u}\\ \vdots \\ {{\rho _{ns}}u}\\ {\rho {u^2} + p}\\ {\rho uv}\\ {(E + p)u}\\ {\rho {e_{{\rm{vib}}}}u} \end{array}} \right], \\ {\bf{G}} = \left[ {\begin{array}{*{20}{c}} {{\rho _1}v}\\ {{\rho _2}v}\\ \vdots \\ {{\rho _{ns}}v}\\ {\rho uv}\\ {\rho {v^2} + p}\\ {(E + p)v}\\ {\rho {e_{{\rm{vib}}}}v} \end{array}} \right], \;\;\;{\kern 1pt} {\bf{H}} = \left[ {\begin{array}{*{20}{c}} {{\rho _1}v}\\ {{\rho _2}v}\\ \vdots \\ {{\rho _{ns}}v}\\ {\rho uv}\\ {\rho {v^2}}\\ {(E + p)v}\\ {\rho {e_{{\rm{vib}}}}v} \end{array}} \right], \\ {\bf{S}} = \left[ {\begin{array}{*{20}{c}} {{{\dot \omega }_1}}\\ {{{\dot \omega }_2}}\\ \vdots \\ {{{\dot \omega }_{ns}}}\\ 0\\ 0\\ 0\\ {{{\dot \omega }_{{\rm{ve}}}}} \end{array}} \right], \\ {\mathit{\boldsymbol{F}}_v} = \left[ {\begin{array}{*{20}{c}} {{q_{x1}}}\\ {{q_{x2}}}\\ \vdots \\ {{q_{xns}}}\\ {{\tau _{xx}}}\\ {{\tau _{xy}}}\\ {u{\tau _{xx}} + \nu {\tau _{xy}} - {q_x} - {q_{{\rm{vib}}, x}} - \sum\limits_{i = 1}^{ns} \rho {h_i}{D_{im}}\frac{{\partial {c_i}}}{{\partial x}}}\\ { - {q_{{\rm{vib}}, x}} - \sum\limits_{i = 1}^{ns} \rho {e_{{\rm{vib}}, i}}{D_{im}}\frac{{\partial {c_i}}}{{\partial x}}} \end{array}} \right], \\ {\mathit{\boldsymbol{G}}_v} = \left[ {\begin{array}{*{20}{c}} {{q_{y1}}}\\ {{q_{y2}}}\\ \vdots \\ {{q_{yns}}}\\ {{\tau _{yx}}}\\ {{\tau _{yy}}}\\ {u{\tau _{yx}} + \nu {\tau _{yy}} - {q_y} - {q_{{\rm{vib}}, y}} - \sum\limits_{i = 1}^{ns} \rho {h_i}{D_{im}}\frac{{\partial {c_i}}}{{\partial y}}}\\ { - {q_{{\rm{vib}}, y}} - \sum\limits_{i = 1}^{ns} \rho {e_{{\rm{vib}}, i}}{D_{im}}\frac{{\partial {c_i}}}{{\partial y}}} \end{array}} \right], \\ {\mathit{\boldsymbol{H}}_v} = \left[ {\begin{array}{*{20}{c}} {{q_{y1}}}\\ {{q_{y2}}}\\ \vdots \\ {{q_{yns}}}\\ {{\tau _{yx}} - \frac{2}{3}y\frac{\partial }{{\partial x}}\left( {\mu \frac{v}{y}} \right)}\\ {{\tau _{yy}} - {\tau _{\theta \theta }} - \frac{2}{3}\frac{\mu }{y}v - \frac{2}{3}y\frac{\partial }{{\partial y}}\left( {\mu \frac{v}{y}} \right)}\\ {u{\tau _{yx}} + v{\tau _{yy}} - {q_y} - {q_{{\rm{vib}}.y}} - \sum\limits_{i = 1}^{ns} \rho {h_i}{D_{im}}\frac{{\partial {c_i}}}{{\partial y}} - \frac{2}{3}\frac{\mu }{y}{v^2} - \frac{2}{3}y\frac{\partial }{{\partial x}}\left( {\mu \frac{{uv}}{y}} \right) - \frac{2}{3}y\frac{\partial }{{\partial y}}\left( {\mu \frac{{{v^2}}}{y}} \right)}\\ { - {q_{{\rm{vib}}, y}} - \sum\limits_{i = mol} \rho {e_{{\rm{vib}}, i}}{D_{im}}\frac{{\partial {c_i}}}{{\partial y}}} \end{array}} \right], \\ \left[ {\begin{array}{*{20}{l}} {{\tau _{xx}}}&{{\tau _{xy}}}\\ {{\tau _{yx}}}&{{\tau _{yy}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {\frac{4}{3}\frac{{\partial u}}{{\partial x}} - \frac{2}{3}\frac{{\partial v}}{{\partial y}}}&{\frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}}}\\ {\frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}}}&{\frac{4}{3}\frac{{\partial v}}{{\partial y}} - \frac{2}{3}\frac{{\partial u}}{{\partial x}}} \end{array}} \right], \\ {\tau _{\theta \theta }} = \mu \left[ { - \frac{2}{3}\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}}} \right) + \frac{4}{3}\frac{v}{y}} \right] \end{array} $

(2)

以上各式中，$p, \rho = \sum\limits_{i = 1}^{ns} {{\rho _i}} = \sum\limits_{i = 1}^{ns} \rho {c_i}$为混合气体的压力和密度，c_i、ρ_i、ns分别为各组分的质量分数、组分密度和组分总数；u、v分别为x、y方向的速度分量；$E = \sum\limits_{i = 1}^{ns} {{\rho _i}} \left( {{e_{{\rm{tr}}, i}} + {e_{{\rm{rot}}, i}} + {e_{{\rm{vib}}, i}} + \Delta h_{{\rm{ref}}, i}^0} \right) + \sum\limits_{i = {\rm{mol}}} {{\rho _i}} {e_{{\rm{vib}}, i}}$为混合气体总能；D_im、h_i、e_{vib, i}分别为组分i的扩散系数、单位质量的焓值和振动能；q_x、q_y、q_{vib, x}、q_{vib, y}分别为平动-转动和振动热流在x、y方向的分量；${\bar \omega _i}, {\bar \omega _{{\rm{ve}}, i}}$分别表示组分i的化学反应速率和振动能源项。

数值模拟采用了热化学非平衡模型，化学反应模型采用Dunn & Kang的5组分17基元反应模型^[14]，热力学模型采用Park的T-T_v双温度模型，流动过程及热化学反应过程由平动温度T和振动温度T_v来描述^[15]。对流项的离散采用AUSMPW+格式，该格式将通量分为对流项和压力通量项，并根据马赫数确定分裂的形式^[16-17]。ASUM类格式在高超声速模拟中有广泛的应用，计算量小且对激波间断及边界层的捕捉具有一定优势，黏性项的离散采用中心差分格式^[18]。对于上述数值方法的可靠性，引用Wieting等在NASA兰利研究中心8英尺高焓风洞上开展的圆柱前缘高超声速气动热试验结果^[19]，圆柱模型半径3.81 cm，长度60.96 cm，试验来流参数如表 1所示，数值模拟与试验结果对比如图 1所示。从圆柱热流分布来看，数值模拟与实验结果吻合良好，验证了该数值方法的可靠性，可以满足本文的计算要求。

计算模型外形选取ELECTRE头部缩比的球锥模型^[20]，头部半径为35 mm，总长为70 mm，锥角为4.6°，如图 2所示。计算网格点数为200×300，壁面边界条件采用等温壁面条件(T=300 K)。

2 40 km高空Ma10-25气动热规律

来流条件选择40 km高空，具体的气流参数如表 2所示，飞行马赫数分别为10、15、20、25。

计算得到4组工况下模型壁面的热流分布如图 3所示。

从图 3(a~d)可以看出，模型头部驻点区为热流峰值，沿着壁面向尾部，热流逐渐减小，对于球锥模型，头部球形区域为主要的高热流区域，随着飞行马赫数提高，热流值变大。

对比研究两种催化条件对热流的影响。Ma=10时，两种催化条件得到的热流分布差别不明显，随着马赫数增加，两种催化条件差别逐渐增大，完全催化高于非催化条件热流值。提取5组马赫数工况下驻点热流值，并计算得到完全催化与非催化条件下驻点热流的比值，如图 4所示。图 4验证了之前的分析，在低马赫数时，两种催化条件下驻点热流率几乎重合，但在高马赫时差异明显。完全催化对比非催化条件下驻点热流的比值也随着马赫数增加而增加，说明完全催化条件在马赫数较高时，有显著增加壁面热流的效应，且随着马赫数提高，这种效应在加强。分析原因为：随着飞行马赫数的提高，模型头部激波的强度随之增大，导致波后气体离解程度加大，离解出的原子成分比重增加，造成壁面处有更多的原子进行催化复合释热，因此壁面催化效应的作用得到加强，对热流的提升也更为显著。

提取Ma=10、20模型驻点线气流粒子质量分数分布，如图 5所示。图 5(a)马赫数较低时，激波强度较弱，波后空气分子几乎没有离解，原子成分的质量比重可以忽略，壁面处原子复合释热效应微弱，因此壁面催化条件对热流的影响不大。随着马赫数提高，如图 5(b)中Ma=20时，由于强激波对波后气体压缩升温，O₂基本完全离解成O，N₂也有较大的离解程度，波后气体中存在着相当比重的原子成分，在壁面附近，即图中各组分质量分布曲线的最右侧位置，完全催化条件下原子完全复合成分子状态，N、O和NO原子的比重为0，而非催化条件下原子成分仍保持一定比重。表 3是非催化条件驻点处粒子质量分数表，可以看出，随着马赫数提高，驻点处的原子质量比重增加，而相对应的完全催化条件下原子质量比重则为0，即两种催化条件下原子质量比重的差异在不断加大，由于原子复合反应释热的效应不断得到增强，因此壁面催化效应对热流的影响随马赫数增加而加强，这再次验证了之前的分析。

上述分析是从总热流率着手进行的，而总热流率是由平动、振动及扩散热流组成，如公式(3)：

$ {q_w} = - k\frac{{\partial T}}{{\partial n}} - {k_v}\frac{{\partial {T_v}}}{{\partial v}} - \rho \sum\limits_{i = 1}^{ns} {{D_i}} {h_i}\frac{{\partial {c_i}}}{{\partial n}} $

(3)

这三部分热流分别由平动温度、振动温度及粒子浓度梯度引起。非催化条件下壁面处不存在粒子浓度梯度，因此扩散热流项为0。完全催化条件下驻点处三部分热流占总热流的比值如图 6所示，总热流主要来自平动热流及扩散热流，随马赫数增加，扩散热流对总热流的贡献在增加，在Ma=25时，总热流有63%来自于扩散热流，因此非催化条件下，认为扩散热流为始终为0是不合适的，尤其是在马赫数较高的时候。

提取模型驻点线平动、振动温度分布曲线如图 7所示，可以看出，两种催化条件下温度曲线的整体分布趋于一致，在马赫数较低时，以图 7(a)中Ma=10为例，平动、振动温度曲线的分布是完全重合的；马赫数增加，以图 7(b)中Ma=20为例，完全催化条件下两条温度曲线向壁面方向整体平移，激波层有变薄的趋势。分析原因，马赫数较高时，激波层内气体离解程度较大，完全催化条件下，壁面处原子完全复合成分子，形成了原子浓度梯度，促使壁面附近的原子成分向壁面处扩散，间接导致激波波后分子离解出的原子也会向壁面附近扩散，并且会促进激波层中气体离解的程度，形成一个由外层向内层的扩散趋势，最终会导致激波脱体距离减小，激波层变薄。因此壁面催化条件在马赫数较高时，不仅仅会影响壁面附近的流场特性及组分分布状态，而且间接对整个激波层都有影响作用。

3 40~70 km高空Ma 20气动热规律

选择固定飞行马赫数20，来流条件分布选择40、50、60、70 km高空气流参数，如表 4所示，模型外形不变，研究不同高空气流参数对气动热的影响规律。

图 8(a~c)分别为50、60、70 km高空两种催化条件下壁面热流的结果，并加入上节中已经得到的40 km结果一起进行比较分析。可以看出，随着飞行高度的增加，模型壁面热流值逐渐减小。

对比壁面催化条件对热流的影响，70 km高空两种催化条件对热流的影响不大，壁面热流分布线几乎重合，但是在其余三组高度均表现出显著的差异。提取四组工况驻点热流，并计算得到完全催化对比非催化条件驻点热流比值，如图 9所示。图 9中随着飞行高度增加，两种催化条件驻点热流的绝对差值减小，在70 km高空处趋近于重合，而驻点热流的比值线出现一个拐点，40 km高空的驻点热流率比值小于50 km高空，50 km高空以上该比值随高度增加而减小。分析具体原因，提取50、70 km高空模型驻点线气流粒子组分质量分数，如图 10所示，对比图 5(b)和图 10(a)可以看出，在40、50 km高空情况，激波层内O₂、N₂都有很大程度的离解，40 km高空气流密度较大，增加了粒子碰撞几率，激波层内流动更接近化学平衡，在气流到达壁面之前大部分原子已经完成复合，离解度相对较低，到壁面处只有少量原子因为壁面催化效应完成复合；而50 km高空气流密度较小，激波层内非平衡程度较强，离解度相对较大，大部分原子成分是一直保持到壁面处才由壁面催化效应完成复合的，即壁面催化效应对50 km高空热流的影响更为显著。而后随着飞行高度的继续增加，气流密度降低明显，导致激波层内离解程度微弱，如图 11(b)中，70 km高空处气流中几乎没有原子组分，则壁面催化效应对热流的影响程度在不断减弱，因此50 km高空以上两种催化条件下热流的比值随高度增加而减小。将非催化条件下模型驻点处粒子组分质量分数列成表 5，从表中可以更具体的看出，40 km原子组分质量比重小于50 km，50 km高空以上原子组分质量比重随高度增加而减小，再次验证了之前的分析。

采取上节的研究方法，分析完全催化条件下总热流中各部分热流的百分比，如图 11所示。总热流的主要贡献仍然来自于平动热流及扩散热流，在40、50 km高空扩散热流百分比高达59%和61%，之后比重开始下降，在70 km处基本可以忽略。图 12是50、70 km高空驻点线平动、振动温度分布曲线，激波层的离解程度主要受到密度影响，因此70 km高空流场特性及气流组分分布基本不受催化条件的影响，而在高度较低时，完全催化条件使两组温度曲线向壁面方向平移，激波脱体距离减小，激波层变薄，具体原因见上节的分析。

4 结论

本文采用数值模拟方法对高超声速球锥模型在非催化和完全催化两种壁面条件下的气动热进行了计算分析，研究壁面催化效应对气动热的影响规律，得到了以下结论：

1) 在40 km高空、马赫数10~25条件下，随着马赫数提高，球锥模型热流增加，完全催化对比非催化条件热流比值增大。这是因为激波层气体离解程度随马赫数提高而加大，导致壁面边界层有更多的原子成分，壁面处原子由催化效应复合释热程度加大，壁面催化效应对热流的影响加强, 而且马赫数较高时，完全催化条件还会影响到整个激波层内的流场特性及气体组分分布状态。

2) 在马赫数20、高度40~70 km条件下，随着高度提升，球锥模型热流减小，且50 km高空以上完全催化对比非催化热流比值减小。这是因为气流密度随高度增加减小，激波层气流离解程度受密度影响而减弱，壁面催化效应对热流的影响减弱。而40、50 km对比出现的拐点同样是受到密度的影响，在离解程度都相当的情况下，密度增加反而会使得原子到达壁面之前就完成复合，导致壁面催化效应的作用减弱。