1 引言

高超声速飞行器所处的环境，空气密度低、飞行速度高，周围的空气经过激波强烈压缩而急剧升温，空气中的氧分子与氮分子发生离解、电离等复杂的化学反应，此时，量热气体与平衡气体已经不再适用，需要使用考虑化学反应的非平衡气体模型进行模拟。高超声速热环境的准确预测对于飞行器的设计起到了至关重要的作用^[1]。目前进行高超声速气动热预测的方法主要三类：（1）纯工程方法；（2）纯数值方法；（3）基于Prandtl边界层理论的边界层外无粘外流解与边界层内工程方法相结合的方法^[2~4]。针对连续流区域气动热的计算，纯数值方法是通过求解Navier-Stokes方程来实现的，高超声速气动热的数值模拟一直是CFD（Computational Fluid Dynamics）研究的难点^[5~7]。

目前对高超声速飞行器进行气动热环境预测的CFD程序主要是基于有限体积法。在实际工程计算中，应用广泛的是多块结构网格程序，国外程序主要包括NASA Langly研究中心开发的LAURA^[8]，NASA Ames研究所开发的DPLP^[9]；在国内，主要是中国空气动力研究与发展中心的董维中、高铁锁、李海燕等^[10~12]、国防科技大学的李桦、潘沙、柳军等^{[13, 14]}都开发了相应的程序，并对气动热环境计算当中存在的数值离散格式、物理模型，边界条件等问题进行了相应的研究。对于复杂外形飞行器，相比较于结构网格，非结构网格在网格前处理、并行计算以及网格自适应等^[15~17]方面更具优势，但是目前针对高超声速气动热环境进行预测的非结构的程序很少，比较成熟的是NASA Langly研究中心开发的FUN3D^[18]程序。在国内，南京航空航天大学的王江峰，伍贻兆，余奇华^{[19, 20]}采用空间上为Jameson格式、时间推进上为五步Runge-Kutta显式格式的格心型非结构网格有限体积法对高超声速化学非平衡流以及热化学非平衡流开展了流场计算以及并行效率的研究。因而进一步研究非结构网格技术在高超声速气动热预测中的应用，以及开发相应的程序是很有必要的。

格心型与格点型格式是有限体积法构造控制体的两种方法，两种方法各有特点（优劣仍没有确定结论）^{[21, 22]}。格心型非结构有限体积法计算时所需的内存空间和时间约为格点型格式的两倍^[23]。对于考虑化学反应的气动热环境模拟，计算量很大，本文选用格点型非结构网格对高超声速化学非平衡流动进行模拟。相对于格心型格式，格点型格式要面临物理边界处理的问题，当控制体仅有一半左右留在边界上时，格心型网格的边界处理的误差储存于网格的中心，而格点型网格的边界处理产生的误差寄存于顶点上，也即将会直接在边界上产生新的离散误差。

对于高超声速化学非平衡粘性流体的定常计算，化学反应源项会带来刚性以及稳定型问题^{[13, 24~26]}，为此，本文选择隐式时间推进方法进行计算。对于格点型非结构网格，边界处理一直是棘手处理的问题，尤其是边界需要隐式处理的时候，不能像格心型那样使用虚网格技术进行处理，为了进行隐式强制处理边界条件，需要对隐式矩阵进行特殊处理。因此，对非结构格点型有限体积法中催化壁面隐式处理的研究很有必要。

本文针对格点型格式在隐式处理完全催化壁面上存在的问题，提出了强制边界条件矩阵的“松弛”形式，使得边界上满足完全催化的强制条件而不出现非物理解。基于此方法，开发了非结构网格格点型有限体积法的程序。在此基础上，分别计算了高超声速圆柱在等温壁与辐射平衡壁条件下完全催化壁面下的热流，并且与结构网格计算结果进行了比较分析，证明了本文提出的“松弛”方法是可行的、有效的。

2 控制方程及物理模型 2.1 控制方程

高超声速连续流的主要控制方程为考虑化学非平衡的三维可压缩Navier-Stokes方程，其守恒型积分形式具体可以写为^[23]

$ \iiint {\frac{{\partial \mathit{\boldsymbol{q}}}}{{\partial t}}{\rm{d}}\mathit{\Omega }} + \iint {\left( {\mathit{\boldsymbol{g}} - \mathit{\boldsymbol{h}}} \right){\rm{d}}\sigma } = \iiint {\mathit{\boldsymbol{\dot \omega }}{\rm{d}}\mathit{\Omega }} $

(1)

式中Ω为控制体，dσ为控制体单元面；q是守恒变量，g和h分别为对流通量与粘性通量，$ \mathit{\boldsymbol{\dot \omega }} $为源项，其表达式为

$ \mathit{\boldsymbol{q}} = \left[ \begin{array}{l} {\rho _s}\\ \rho u\\ \rho v\\ \rho w\\ \rho E \end{array} \right];\mathit{\boldsymbol{g}} = \left[ \begin{array}{l} {\rho _s}U\\ \rho Uu + p{n_x}\\ \rho Uv + p{n_y}\\ \rho Uw + p{n_z}\\ \rho UH \end{array} \right] $

(2)

$ \mathit{\boldsymbol{h}} = \left[ \begin{array}{l} {n_x}{\mathit{\Phi }_{x,s}} + {n_y}{\mathit{\Phi }_{y,s}} + {n_z}{\mathit{\Phi }_{z,s}}\\ {n_x}{\tau _{xx}} + {n_y}{\tau _{xy}} + {n_z}{\tau _{xz}}\\ {n_x}{\tau _{yx}} + {n_y}{\tau _{yy}} + {n_z}{\tau _{yz}}\\ {n_x}{\tau _{zx}} + {n_y}{\tau _{zy}} + {n_z}{\tau _{zz}}\\ {n_x}{\mathit{\Theta }_x} + {n_y}{\mathit{\Theta }_y} + {n_z}{\mathit{\Theta }_z} \end{array} \right] $

(3)

$ \mathit{\boldsymbol{\dot \omega = }}\left[ \begin{array}{l} {{\dot \omega }_s}\\ 0\\ 0\\ 0\\ 0 \end{array} \right] $

(4)

式中$ {\dot \omega _s} $为组分s的质量生成率，ρ_s为组分s的密度，ρ为混合气体的密度，u, v, w分别为x, y, z三个方向上的速度分量，E为单位质量单位体积内的总能量，U为控制体面上的法向速度，p为压强，H为单位质量单位体积内的总焓，τ_ij表示应力张量的9个分量，n_x, n_y, n_z表示控制体面上法向矢量的三个分量，其它变量的表达式如下所示

$ {\mathit{\Phi }_{x,s}} = \rho {D_s}\frac{{\partial {y_s}}}{{\partial x}} $

(5)

$ {\mathit{\Phi }_{y,s}} = \rho {D_s}\frac{{\partial {y_s}}}{{\partial y}} $

(6)

$ {\mathit{\Phi }_{z,s}} = \rho {D_s}\frac{{\partial {y_s}}}{{\partial z}} $

(7)

$ {\mathit{\Theta }_x} = u{\tau _{xx}} + v{\tau _{xy}} + w{\tau _{xz}} + \eta \frac{{\partial T}}{{\partial x}} + \rho \sum\limits_{s = 1}^{{N_s}} {{h_s}{D_s}\frac{{\partial {y_s}}}{{\partial x}}} $

(8)

$ {\mathit{\Theta }_y} = u{\tau _{yx}} + v{\tau _{yy}} + w{\tau _{yz}} + \eta \frac{{\partial T}}{{\partial y}} + \rho \sum\limits_{s = 1}^{{N_s}} {{h_s}{D_s}\frac{{\partial {y_s}}}{{\partial y}}} $

(9)

$ {\mathit{\Theta }_z} = u{\tau _{zx}} + v{\tau _{zy}} + w{\tau _{zz}} + \eta \frac{{\partial T}}{{\partial z}} + \rho \sum\limits_{s = 1}^{{N_s}} {{h_s}{D_s}\frac{{\partial {y_s}}}{{\partial z}}} $

(10)

在式（5）~（10）中，y_s表示组分s的摩尔分数，T表示混合气体的温度，η是空气混合物热传导系数，D_s是组分s的摩尔扩散系数，h_s是单位质量组分s的静焓。

2.2 热力学关系

本文所考虑的多组分混合气体可以认为是有多种单一组分热完全气体组成的，混合气体的焓值以及组分焓值根据以下式子进行计算

$ {h_s}\left( T \right) = \frac{R}{{{M_s}}}\left( {\sum\limits_{k = 1}^5 {\frac{{A_k^s{T^k}}}{k}} + A_6^s} \right) $

(11)

$ h = \sum\limits_{i = 1}^{{N_s}} {{h_s}\left( T \right)} $

(12)

式（11）当中的拟合系数可以在文献[27]中找到。

2.3 输运系数

为了使得Navier-Stokes方程能够封闭，需要计算混合气体的输运系数（粘性系数、热传导系数和扩散系数）。考虑包含多组分的非平衡反应的影响，含混气体中的输运过程与量热完全气体有所不同。本文计算方法与文献[24, 27]中计算方法是一致的。

2.4 化学反应源项

对于式（4）中单位体积组分s的质量生成率$ {\dot \omega _s} $，使用有限速率模型的计算公式为

$ {{\dot \omega }_s} = {M_s}\sum\limits_{r = 1}^{{N_r}} {\left( {{\beta _{s,r}} - {\alpha _{s,r}}} \right)\left( {{\beta _{{\rm{f}},r}} - {\alpha _{{\rm{b}},r}}} \right)} $

(13)

式中N_r表示基元反应的总个数，α_{s, r}, β_{s, r}分别表示组分s在第r个反应中反应物与生成物的化学计量数，R_{f, r}, R_{b, r}表示第r个反应中的正反应生成率与逆反应生成率，表达式为

$ {R_{{\rm{f}},r}} = {k_{{\rm{f}},r}}\prod\limits_{s = 1}^{{N_s}} {{{\left( {\frac{{{\rho _s}}}{{{M_s}}}} \right)}^{{\alpha _{s,r}}}}} $

(14)

$ {R_{{\rm{b}},r}} = {k_{{\rm{b}},r}}\prod\limits_{s = 1}^{{N_s}} {{{\left( {\frac{{{\rho _s}}}{{{M_s}}}} \right)}^{{\beta _{s,r}}}}} $

(15)

式中k_{f, r}, k_{b, r}分别表示正反应速率与逆反应速率。本文计算采用5组分（N, N₂, O, O₂, NO）化学反应模型，基元反应的具体描述以及反应速率的计算详见文献[27]。

3 数值离散 3.1 计算格式

本文使用格点型非结构网格算法求解NavierStokes方程，控制体如图 1中的虚线所示（如ABCDEF）。数据结构采用的是基于边的数据结构，在计算对流通量的时候只需要对所有的边（如图 1中蓝色虚线所示）进行循环即可。

对流项采用具有二阶精度的迎风型STVD格式^[16]，具体见式（16）~（26）所示，式子中的下标L，R表示的是控制面两侧，图中蓝色虚线所表示的控制面（AB）上的通量的计算时使用的q_R，q_L即是节点1、3上的守恒变量值。粘性项梯度采用Green-Gauss公式进行计算，熵修正函数使用的是Harten-Yee^[28]型修正方法。

由于模拟的是定常流动，为了避免对组分生成源项采用显式求解时带来的刚性问题，本文采用隐式求解方法，源项采用点隐进行处理。使用基于线性欧拉向后时间差分的时间推进方式^[18]，每一时间步上的线性方程组如式（25）所示，使用Gauss-Seidel迭代求解该线性方程组。式（25）与（26）中的Rhs表示的是控制体上的通量，V表示的是控制体Ω的体积，I是单位阵，Δt为时间步长；式（26）中的[D]ⁿ与[O]ⁿ表示[A]ⁿ矩阵中的对角部分与非对角部分。

$ \mathit{\boldsymbol{g}} = \frac{1}{2}\left[ {{\mathit{\boldsymbol{g}}_{\rm{L}}} + {\mathit{\boldsymbol{g}}_{\rm{R}}} - \mathit{\boldsymbol{R}}\left| \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} \right|\left( {{\rm{d}}\mathit{\boldsymbol{\tilde q}} - {\rm{d}}{{\mathit{\boldsymbol{\tilde q}}}_{\lim }}} \right)} \right] $

(16)

$ {\rm{d}}\mathit{\boldsymbol{\tilde q}} = {\mathit{\boldsymbol{R}}^{ - 1}}\left( {{\mathit{\boldsymbol{q}}_{\rm{R}}} - {\mathit{\boldsymbol{q}}_{\rm{L}}}} \right) $

(17)

$ {\rm{d}}{{\mathit{\boldsymbol{\tilde q}}}_{\lim }} = 2\varphi \left( {{p_{\rm{L}}},{p_{\rm{R}}}} \right)\mathit{\Psi } $

(18)

$ \mathit{\Psi } = \min \bmod \left( {{\rm{d}}{{\mathit{\boldsymbol{\tilde q}}}_{\rm{L}}},{\rm{d}}{{\mathit{\boldsymbol{\tilde q}}}_{\rm{R}}},{\rm{d}}\mathit{\boldsymbol{\tilde q}},\frac{1}{4}\left( {{\rm{d}}{{\mathit{\boldsymbol{\tilde q}}}_{\rm{L}}} + {\rm{d}}{{\mathit{\boldsymbol{\tilde q}}}_{\rm{R}}}} \right)} \right) $

(19)

$ {\rm{d}}{{\mathit{\boldsymbol{\tilde q}}}_{\rm{M}}} = {R^{ - 1}}\left[ {{\rm{d}}\mathit{\boldsymbol{s}} \cdot \nabla {\mathit{\boldsymbol{q}}_{\rm{M}}}} \right] $

(20)

$ {\rm{d}}\mathit{\boldsymbol{s}} = \left( {{x_{\rm{L}}} - {x_{\rm{R}}}} \right)\mathit{\boldsymbol{i}} + \left( {{y_{\rm{L}}} - {y_{\rm{R}}}} \right)\mathit{\boldsymbol{j + }}\left( {{z_{\rm{L}}} - {z_{\rm{R}}}} \right)\mathit{\boldsymbol{k}} $

(21)

$ \varphi \left( {{p_{\rm{L}}},{p_{\rm{R}}}} \right) = \frac{1}{2}\left[ {1 - \cos \left( {{p_{{\rm{rat}}}}{\rm{ \mathsf{ π} }}} \right)} \right] $

(22)

$ \hat \varphi = \max \left( {\frac{{{p_{\rm{L}}}}}{{{p_{\rm{R}}}}},\frac{{{p_{\rm{R}}}}}{{{p_{\rm{L}}}}}} \right) $

(23)

$ {p_{{\rm{rat}}}} = \min \left[ {1,\max \left( {0,\frac{{{{\hat \varphi }_{\max }} - \hat \varphi }}{{{{\hat \varphi }_{\max }} - {{\hat \varphi }_{\min }}}}} \right)} \right] $

(24)

$ {\left[ \mathit{\boldsymbol{A}} \right]^n}{\left\{ {\Delta \mathit{\boldsymbol{q}}} \right\}^n} = - {\left\{ {\mathit{\boldsymbol{Rhs}}} \right\}^n} $

(25)

$ {\left[ \mathit{\boldsymbol{A}} \right]^n} = {\left[ \mathit{\boldsymbol{D}} \right]^n} + {\left[ \mathit{\boldsymbol{O}} \right]^n} = \frac{V}{{\Delta t}}\mathit{\boldsymbol{I}} + \frac{{\partial {{\left\{ {\mathit{\boldsymbol{Rhs}}} \right\}}^n}}}{{\partial \mathit{\boldsymbol{q}}}} $

(26)

3.2 边界处理

本文所研究的范围设置在连续流领域，因而固体壁面上采用无滑移壁面。在不考虑烧蚀与固体热传导的情况下，常常给定等温壁或辐射平衡壁面，分别见式（27），（28）所示

$ T = {T_{\rm{w}}} $

(27)

$ {q_{\rm{w}}} = \varepsilon \sigma T_{\rm{w}}^4 $

(28)

同样地，对于组元连续方程，在壁面处需要给出边界条件，文献[16, 24]中的描述，一般给定的是完全催化壁面、非催化壁面还有有限催化壁面。本文着重考虑的是完全催化壁面边界条件的施加。在文献[15, 16, 24]中，完全催化壁面定义的是在表面处各个组元的质量分数与自由来流一致，即

$ {c_{i{\rm{w}}}} = {c_{i\infty }} $

(29)

式中$ {c_i} = \frac{{{\rho _i}}}{\rho } $表示的是组元i的质量分数。

对于格心型有限体积法，式（27）~（29）提出来的壁面边界条件可以通过虚网格技术来实现；对于格点型有限体积法，在边界上只有半个控制体，很多时候是将边界点当作内点处理，同时也没有做隐式处理^[29]。图 2中由红色线条与绿色线条围成的多边形（M₁ -M₂ -M₃ -M₄ -M₅）即为边界节点P_b所在的控制体。

本文使用隐式时间推进方式，由于P_b不是内部节点，在边界上需要对落在半个控制体上的点P_b进行特殊处理以达到边界强制的目的。换句话说，就是需要对式（26）中的[A]ⁿ进行特殊处理，对于内部点，一般情况下[O]ⁿ ≠ 0，但是对于边界点，[O]ⁿ = 0。格点型格式下，根据文献[18]中对于量热气体，等温壁面的隐式强制处理方式如式（30）所示。式（30）表示的是Δu = Δv = Δw = ΔT = 0，其中D₁₁，D₁₂，D₁₃，D₁₄，D₁₅表示的是式（26）中[D]ⁿ的元素，其计算方法与内点一致。换句话说，式（30）是式（25）的特殊形式，边界初始化后，无滑移等温壁面的条件在边界上每个点都满足。对于化学非平衡气体，等温壁面上需要同时满足式（27），（29）两个式子；边界初始化后为了使得在迭代过程中边界上满足式（29），需要在每个时间步上组元的质量分数的增量为0，即Δc_i = 0，进而推导出组元质量增量的关系，如式（31）所示。将式（31）写成矩阵形式，并结合式（30）中关于动量方程以及能量方程的约束条件，可以得到边界上强制隐式方程形式，见式（32）。

$ \left[ {\begin{array}{*{20}{c}} {{D_{11}}}&{{D_{12}}}&{{D_{13}}}&{{D_{14}}}&{{D_{15}}}\\ { - u}&1&0&0&0\\ { - v}&0&1&0&0\\ { - w}&0&0&1&0\\ { - {C_{\rm{V}}}{T_{\rm{w}}} + \frac{1}{2}\left( {{u^2} + {v^2} + {w^2}} \right)}&{ - u}&{ - v}&{ - w}&1 \end{array}} \right]\\ \left[ \begin{array}{l} \Delta \rho \\ \Delta \rho u\\ \Delta \rho v\\ \Delta \rho w\\ \Delta \rho E \end{array} \right] = \left[ \begin{array}{l} Rh{s_1}\\ 0\\ 0\\ 0\\ 0 \end{array} \right] $

(30)

$ \Delta {c_i} = \Delta \left( {\frac{{{\rho _i}}}{\rho }} \right) = \frac{1}{\rho }\left( {{\delta _{ij}}\Delta {\rho _j} - {c_{i\infty }}\Delta {\rho _i}} \right) = 0 $

(31)

$ {\left[ \mathit{\boldsymbol{A}} \right]^n}{\left\{ {\Delta \mathit{\boldsymbol{q}}} \right\}^n} = \left[ {\begin{array}{*{20}{c}} {1 - {c_{i\infty }}}&0&0&0&0\\ { - u}&1&0&0&0\\ { - v}&0&1&0&0\\ { - w}&0&0&1&0\\ {\frac{{\partial T}}{{\partial {\rho _i}}}}&{\frac{{\partial T}}{{\partial \rho u}}}&{\frac{{\partial T}}{{\partial \rho v}}}&{\frac{{\partial T}}{{\partial \rho w}}}&{\frac{{\partial T}}{{\partial \rho E}}} \end{array}} \right]\left[ \begin{array}{l} \Delta {\rho _i}\\ \Delta \rho u\\ \Delta \rho v\\ \Delta \rho w\\ \Delta \rho E \end{array} \right] = \left[ \begin{array}{l} 0\\ 0\\ 0\\ 0\\ 0 \end{array} \right] $

(32)

对于边界上的节点，式（32）表示壁面节点上每个时间步上线性方程组，也就是式（25）在等温、完全催化壁上的具体体现。可以清楚地看出，对于边界节点而言，在每个时间步上都准确满足Δρ_i = Δρu = Δρv = Δρw = ΔT = 0。这个结果与真实物理情况是相悖的，因为壁面处的密度跟压强是与来流不同的，并且在时间推进过程中是不断变化的。得到这种非物理现象的原因在于式（32）当中的[A]ⁿ是非奇异的，并且右端项为零向量，因此，无论内点如何计算，边界上的点都是不满足物理现象的。为了使得点隐格式的强制方式能够有效使用，需要对式（32）进行特殊处理一下。

因为所有组元的质量分数满足以下式子

$ \sum\limits_{s = 1}^{{N_s}} {{c_s} = 1} $

(33)

即所有组元的质量分数之和为1，故式（32）中对组元质量分数的约束是过约束，这样就造成了计算得到的组元质量是非物理的，换句话说，式（32）是壁面边界条件式（29）成立的充分但是不必要条件。为此，需要将式（19）中多余的约束去掉，物理条件的约束是N_s个，而式（32），（33）两个式子表示的是N_s + 1个约束条件，因此，需要将式（32）中对于Δρ_i的约束去掉一个。对于式（32）中的第j个组分方程的强制条件进行“松弛”，具体形式见式（34）所示，其中，$ j = \mathop {\max }\limits_{i = 1}^{{N_s}} \left( {{c_{i\infty }}} \right) $。

式（34）与（32）两个式子不同之处在于对第j个组元方程的边界上的条件不是强制的，其中，D_j，D_u，D_v，D_w，D_e表示的是式（26）中[D]ⁿ的元素，其计算方法与内点一致，Rhs_j表示的是第j个组元的组分扩散方程计算得到的残差，其具体做法与式（30）当中连续方程的处理是一致的。经过“松弛”处理后的式（34）既能够满足无滑移、等温度壁面，同时也能满足N_s -1组分的质量分数满足式（29）；由于式（33），在最后收敛的时候，对于第j个组分也是满足式（29）的。

$ {\left[ \mathit{\boldsymbol{A}} \right]^n}{\left\{ {\Delta \mathit{\boldsymbol{q}}} \right\}^n} = \left[ {\begin{array}{*{20}{c}} {1 - {c_{k\infty }}}&0&0&0&0\\ {{D_j}}&{{D_{\rm{u}}}}&{{D_{\rm{v}}}}&{{D_{\rm{w}}}}&{{D_{\rm{e}}}}\\ { - u}&1&0&0&0\\ { - v}&0&1&0&0\\ { - w}&0&0&1&0\\ {\frac{{\partial T}}{{\partial {\rho _i}}}}&{\frac{{\partial T}}{{\partial \rho u}}}&{\frac{{\partial T}}{{\partial \rho v}}}&{\frac{{\partial T}}{{\partial \rho w}}}&{\frac{{\partial T}}{{\partial \rho E}}} \end{array}} \right]\left[ \begin{array}{l} \Delta {\rho _k}\\ \Delta {\rho _j}\\ \Delta \rho u\\ \Delta \rho v\\ \Delta \rho w\\ \Delta \rho E \end{array} \right] = \left[ \begin{array}{l} 0\\ Rh{s_j}\\ 0\\ 0\\ 0\\ 0 \end{array} \right] $

(34)

对于辐射平衡壁面条件的施加，可以类比于等温壁，能量方程项的强制施加可以保留式（34）的形式，在每个时间步上满足ΔTⁿ⁺¹ = 0，只需要在计算时，依据壁面热流给出壁面温度

$ {T_{{\rm{new}}}} = \sqrt[4]{{\frac{{q_{\rm{w}}^n}}{{\varepsilon \sigma }}}} $

(35)

$ {T^{n + 1}} = \delta {T_{{\rm{new}}}} + \left( {1 - \delta } \right){T^n} $

(36)

式中ε是表面材料辐射系数，本文中给定ε = 0.89；σ是Stefan-Boltzmann常数，σ = 5.67 × 10^-8W/(m²∙K⁴)；δ为迭代松弛因子，一般情况下取δ = 0.01。

4 算例验证与讨论

基于上一节中的理论方法，开发了非结构有限体积法程序UNFLOW。为了验证本文所提出处理完全催化壁面边界条件的方法是可行的，计算半径为1m的圆柱的化学非平衡绕流，得到壁面热流与壁面上的各个组元的质量分数，将得到的结果与结构化网格程序AHEAT的计算结果进行对比。AHEAT是可以求解量热气体、平衡气体、化学非平衡以及热化学非平衡气体的结构化网格程序，之前的文章中已经对其正确性进行了验证^{[30, 31]}。所以在本文中选用AHEAT计算的结果进行对比验证。

4.1 等温条件下圆柱扰流问题

高超声速飞行器的翼前缘、舵面等受热严重的区域可以近似为圆柱面，为了更为简捷地验证方法的正确性，选择半径为1m的圆柱，在H = 53 km高空处的来流条件（T_∞ = 265.0K，p_∞ = 53.0 Pa）下，速度为V_∞ = 4230m/s，攻角为0°时的气动热环境的预测。在该条件下，空气发生化学反应，选择五组分（N, N₂, O, O₂, NO）进行计算，壁面边界条件为等温壁T_w = 343K，完全催化壁。自由来流中各个组元的质量分数见表 1所示。

首先，使用结构网格程序AHEAT进行计算，分别采用热平衡（单温）与热非平衡（双温）条件进行计算。计算热流时壁面法向网格尺寸的选择有不同的方法^{[13, 32, 33]}，为了使得计算得到的壁面热流可靠准确，选择文献[13]中所述准则，本文中所有算例的壁面处网格雷诺数均保持在2~6；在验证网格无关性的基础上（网格无关性过程在此不列出）选择环向网格为79，垂直于壁面方向的网格为78的网格进行计算，需要指出的是，采用的网格在对称方向上的网格数为1，计算网格见图 3所示。

计算结果见图 4所示，图中“Structured，1-T”表示的是结构网格程序AHEAT采用热平衡（单温）计算的结果；而“Structured，2-T”表示的是结构网格程序AHEAT采用热非平衡（双温）计算的结果。为方便起见，本文曲线图中“1-T”表示计算采用单温模型进行计算，而“2-T”表示计算采用双温模型进行计算。可以看出，采用单温与双温计算得到的壁面热流是一致的，说明圆柱所处流场中的热非平衡很不明显，可以使用热平衡（单温）方程进行计算。

使用本文所提出的施加边界条件的方法进行计算，采用的网格见图 5所示。需要说明的是，非结构所采用的网格与结构网格的区别仅仅是网格类型不同，非结构网格在对称方向上有一层网格，非结构网格全部为四面体，节点坐标是完全一致的。结构与非结构使用单温进行计算的温度场云图见图 6所示。图 4中，“Unstructured，1-T”表示的是非结构程序UNFLOW采用热平衡（单温）计算的结果，从图中可以看出，除了在x = 0处的热流峰值有较小差别外，其他各处的热流值与结构网格计算结果是相吻合的。本文对完全催化边界条件进行强制时，采用的方法是对自由来流当中质量分数最大组分的方程进行“松弛”处理的，对于收敛后的结果，壁面处的组分质量分数是否与自由来流一致主要取决于质量分数最大的N₂，图 7给出了壁面上N₂的质量分数以及N₂所在组分扩散方程的残差沿壁面的分布。可以清楚地看出，壁面处N₂的质量分数与来流的是吻合的，N₂所在组分扩散方程的残差虽然不像其他组分那样强制为0，但是在迭代计算后，它的残差也达到了10^-8量级，说明了该方程也收敛了。

4.2 辐射平衡条件下圆柱扰流问题

本节算例采用与上一节同样的网格、工况，不同之处在于壁面不是等温壁面，而是辐射平衡壁面。分别使用结构网格程序AHEAT和非结构网格程序UNFLOW进行计算。结构网格程序AHEAT分别采用热平衡（单温）与热非平衡（双温）计算。在图 8中可以清楚地看出，壁面温度与热流在单温与双温条件下的结果是一致的，说明在该工况下，热力学是平衡的，可以采用单温进行模拟热环境。

非结构网格下采用单温五组分进行计算，结构网格与非结构网格计算得到的温度云图如图 9所示，从图中可以看出，非结构网格计算得到的温度场与结构网格所得到的结果吻合很好。与结构网格计算得到的壁面热流和温度进行比较，从图 8中可以看出，吻合得很好，说明本文提出的方法可以有效地处理辐射平衡壁面下的完全催化壁面条件。同样地，可以绘出壁面上N₂的质量分数还有N₂组分方程的残差曲线，见图 10所示。

4.3 等温条件下球头扰流问题

本节算例选择高超声速飞行器中典型部位——球头进行计算，球头半径为R = 1m。来流速度V_∞ = 5 km /s，密度为ρ_∞ = 0.001kg/m³，温度T_∞ = 200K，攻角为0°。该条件下，空气发生化学反应，选择五组分（N, N₂, O, O₂, NO）进行计算，壁面边界条件为等温壁T_w = 500K，完全催化壁。自由来流中各个组元的质量分数见表 1所示。

选择轴对称网格进行计算，将二维平面网格沿对称轴旋转4.5°后得到计算所用网格。AHEAT所使用的结构网格如图 11所示，该网格环向网格量为64，垂直于壁面方向的网格量为100，垂直于轴对称方向的网格量为1；UNFLOW所使用的网格是将AHEAT使用的网格分割为四面体，如图 12所示。

分别使用AHEAT与UNFLOW进行计算，得到的温度云图如图 13所示，可以看出，UNFLOW计算得到的温度场与AHEAT计算得到的温度场是相吻合的。图 14给出了球头壁面上热流的计算结果，从图中可以清楚地看出，二者计算得到的热流相差很小，UNFLOW计算得到的热流是可信的。与前两节类似，图 15给出壁面上各个组元的质量分数以及各个组元方程的收敛残差。因为表 1中质量分数最大的是N₂，最后计算得到的N₂的质量分数与表 1中一致，并且可以看出N₂的组元方程的残差达到了10^-8量级，表明该方程也收敛了。

4.4 讨论

通过上一节的讨论可以看出，本文所提出的边界处理方法对高超声速化学非平衡流动下考虑完全催化壁面条件的热流计算是准确的，壁面处经过“松弛”的组元的质量分数也满足边界条件。3.1与3.2中的自由来流的组元质量分数都如表 1中所示，N₂的质量分数最大，其值达到了0.766，远远大于其他组分的质量分数之和，通过对N₂的组元方程在壁面处的隐式矩阵进行“松弛”处理，计算得到的壁面组元的质量分数满足边界条件，并且热流与结构网格程序计算得到的热流是相吻合的。

如果自由来流中质量分数最大的组元并没有像表 1那样具有明显的质量“占优”，本文提出的“松弛”处理方法是否还能够有效，为了验证自由来流组元质量分数的组成对本文所提方法是否敏感，通过改变自由来流中组元的质量分数来观察壁面处理情况是否能够正确有效，具体算例的自由来流的组元组成见表 2所示。

表 2中三个算例的条件与4.2中算例相比，除了来流的组元成分组成不一样以外，其它条件均一致。这三个算例中N，O，NO的质量分数保持不变的，通过改变N₂与O₂的质量分数，来检测最大质量分数组元所占比重是否影响计算结果。逐渐增大N₂的质量分数，使其与比重最大的O₂的质量分数相差逐步减小，直至最后二者的质量分数均为0.50。

需要关注的是组分在壁面上是否能够满足完全催化壁的条件，图 16中绘出了表 2中各个工况下的壁面处各个组元的质量分数的分布以及三个工况中被“松弛”的组元连续性方程的残差。从图中可以清楚地看出，在前两位组元质量分数相差很小甚至相等的情况下，壁面边界条件仍然能够充分保证，并且经“松弛”处理的组元方程能够收敛。同时，从图 17中可以看出，对于完全催化壁面，在保证壁面质量分数与来流条件一致的情况下，壁面的热流与温度是一定的。

5 结论

本文提出了一种边界上隐式矩阵的“松弛”方法，用以处理高超声速化学非平衡流动中的完全催化壁面条件。在此基础上使用基于格点型非结构有限体积法对圆柱、球头在高超声速下化学非平衡流动进行了模拟，通过对比分析，得到以下结论：

（1）针对高超声速化学非平衡流，本文提出的“松弛”方法可以正确地表征等温与辐射平衡条件下的完全催化壁面，并且计算得到的壁面热流是合理、准确的。

（2）通过改变来流中组元的质量分数来测试本文提出的强制隐式矩阵对来流组元质量分数是否敏感。数值实验表明，本文提出的边界处理方法对来流中比重最大组元的质量分数是否“占优”并不敏感。即使在两个组元质量分数各占50%的情况下，边界“松弛”处理方法也可以准确处理，不会造成算法的失效。