近年来，随着高超声速飞行器技术的迅速发展，其关键技术之一的推进技术及燃烧流场的准确预测得到高度重视。超声速燃烧室中流场推进过程与燃烧过程复杂，并具有强耦合相互作用行为，燃料和空气掺混并发生化学反应时具有显著的三维流动效应，各种复杂的化学反应在时间尺度上差异较大，同时包含了剧烈的激波间断、多组元快速化学反应等一系列复杂流动现象。为详细精确的描述超声速燃烧室流动，就对数值模拟方法的精度提出了较高的要求。国内外在这方面做了大量的研究，将一系列高分辨率格式推广应用于求解可压缩化学反应流场，如TVD^[1]，NND^[2]，AUFS^[3]，HLLE+^[4]及CUSP^[5]等格式。Roe格式^[6]是P L Roe在1981年提出的一种基于Riemann近似解的通量差分分裂格式，该格式以其数值耗散低、接触间断分辨率高、激波捕捉能力强等优点，广泛应用于完全气体流动的求解中。

以理论推导为基础，基于非结构网格数据格式，将Roe格式由量热完全气体扩展至热完全化学反应混合气体。验证算例采用经典的H₂-Air预混绕钝头体激波诱导燃烧流场^[7-10]，比较分析了3种经典化学动力学模型对超声速条件下激波诱导燃烧流场的影响，得到了与参考文献相符的数值模拟和试验结果，表明所发展的方法可以应用于多组元燃烧流场的数值模拟与分析。

1 控制方程

含化学反应源项的守恒型三维多组元Euler方程为

${{\partial U} \over {\partial t}} + {{\partial E} \over {\partial x}} + {{\partial F} \over {\partial y}} + {{\partial G} \over {\partial z}} = {\rm{ }}S$

(1)

式中：U为守恒变量;E，F和G为对流通量;S为化学反应源项，其具体形式见文献^[11]。

2 数值方法

Roe格式^[6]以数值耗散低、接触间断分辨率高、激波捕捉能力强等优点广泛应用于完全气体流动的求解中。本文研究的重点是将Roe格式扩展至化学反应流，用以提高复杂多组元燃烧流场的数值模拟能力。

2.1 Roe通量分裂格式

对于多组元化学反应气体，完全气体状态方程不再适用,需要重新推导。为导出化学反应条件下Roe格式的具体形式，取参数向量

$W = {\rho ^{1/2}}{({Y_1}, \ldots ,{Y_s}, \ldots ,{Y_n},u,v,w,H)^T}$

(2)

式中：Y_s表示组元s的质量分数;n表示组元个数;H为单位质量总焓;u,v,w表示速度矢量在x,y,z 3个方向上的分量。这里仅对x方向对流通量项进行详细描述，其余方向上的结果可以类推得到。对于控制体面I上对流通量项E，有

${E_I} = {1 \over 2}({E_L} + {E_R}) - {1 \over 2}(\Delta {F_A} + \Delta {F_B} + \Delta {F_C})$

(3)

式中

$\eqalign{ & \Delta {F_A} = \left( {\Delta \rho - {{\Delta {p^2}} \over {{{\tilde a}^2}}}} \right)\left| {\tilde u} \right|\left[ {\matrix{ {{{\tilde Y}_1}} \cr \vdots \cr {{{\tilde Y}_s}} \cr \vdots \cr {{{\tilde Y}_n}} \cr {\tilde u} \cr {\tilde v} \cr {\tilde w} \cr {\tilde H - {{{{\tilde a}^2}} \over {\tilde \gamma - 1}}} \cr } } \right] + \tilde \rho \left| {\tilde u} \right|\left[ {\matrix{ {\Delta {Y_1}} \cr \vdots \cr {\Delta {Y_s}} \cr \vdots \cr {\Delta {Y_n}} \cr 0 \cr {\Delta v} \cr {\Delta w} \cr \Phi \cr } } \right] \cr & \Delta {F_{B,C}} = {1 \over {2{{\tilde a}^2}}}(\Delta p \mp \tilde \rho \tilde a\Delta u)\left| {\tilde u \mp \tilde a} \right|\left[ {\matrix{ {{{\tilde Y}_1}} \cr \vdots \cr {{{\tilde Y}_s}} \cr \vdots \cr {{{\tilde Y}_n}} \cr {\tilde u \mp \tilde a} \cr {\tilde v} \cr {\tilde w} \cr {\tilde H \mp \tilde u\tilde a} \cr } } \right] \cr} $

(4)

其中

$\eqalign{ & \Phi = - \sum\limits_{s = 1}^n {{{{{\tilde \varepsilon }_s}} \over \beta }} \Delta {Y_s} + \tilde v\Delta v + \tilde w\Delta w \cr & \beta = {{{R_u}} \over {{C_v}\bar M}},{C_v} = \sum\limits_{s = 1}^n {{{{{\bar W}_s}} \over {({{\bar W}_1} + \ldots + {{\bar W}_n})}}} {C_{v,s}} \cr & {{\tilde a}^2} = \sum\limits_{s = 1}^n {{{\tilde Y}_s}{{\tilde \varepsilon }_s}} + \beta \left[ {{{\tilde H}^2} - \left( {{{\tilde u}^2} + {{\tilde v}^2} + {{\tilde w}^2}} \right)} \right] \cr & {{\tilde \varepsilon }_s} = {{{R_u}} \over {{M_s}}}\tilde T + {1 \over 2}\beta \left( {{{\tilde u}^2} + {{\tilde v}^2} + {{\tilde w}^2}} \right) - \beta {{\tilde e}_s} \cr & \Delta \left( \cdot \right) = {\left( \cdot \right)_L} - {\left( \cdot \right)_R} \cr} $

(5)

式中:矩阵中“^-”表示算术平均，“~”表示Roe平均；R_u表示通用气体常数；M_s表示组元s的分子量；表示混合气体的平均分子量。上述推导得到的控制方程通量的Jacobian矩阵特征值在声速点附近会趋近于0，违反了熵条件，产生非物理解。为解决这一问题，采用Harten-Yee型熵修正方法^[12]，引入额外的耗散来处理奇异的特征值。

2.2 数值离散方法

控制方程离散采用基于非结构网格的有限体积方法，空间离散采用Roe通量分裂格式。远场边界采用基于Riemann不变量的无反射边界条件，物面边界采用无穿透条件与完全非催化壁条件，对称面采用对称边界条件。

2.3 点隐式方法

由于化学反应流场中存在多种时间尺度，在反应剧烈区域，差别可能达到数个量级，不可避免地会出现刚性问题，为解决数值计算中出现的刚性问题，加速收敛，时间推进采用常用的点隐式格式，对化学反应源项做隐式处理^[13-14]。

流场控制方程空间离散后，形式为

${\left( {{{dU} \over {dt}}} \right)^{(n{\rm{ }})}} = - \nabla F_v^{(n)} + {\rm{ }}{S^{(n)}}$

(6)

式中(n)表示第n个时间步的流场参数值。对化学反应源项S隐式处理，即

${\left( {{{dU} \over {dt}}} \right)^{(n)}} = - \nabla F_v^{(n)} + {\rm{ }}{S^{(n + 1)}}$

(7)

对化学反应源项进行Taylor展开，取到Δt项，即

${S^{(n + 1)}} = {S^{(n)}} + {\left( {{{\partial S} \over {\partial U}}} \right)^{(n)}}{\left( {{{\partial U} \over {\partial t}}} \right)^{(n)}}\Delta t + O(\Delta {t^2})$

(8)

忽略Δt 的高阶项，将式(8) 代入式(7) 整理后得到

${\left( {{{dU} \over {dt}}} \right)^{(n)}} = {\left[ {I - {{\left( {{{\partial S} \over {\partial U}}} \right)}^{(n)}}\Delta t} \right]^{ - 1}}( - \nabla F_v^{(n)} + {S^{(n)}})$

(9)

式(9) 即为点隐式处理后的控制方程形式，其中I为单位矩阵，F_v为对流通量项。经过这样的处理后，化学反应源项所引起的刚性问题就可以得到很好地解决。

3 计算结果与分析

当按理想化学当量比预混的H₂-Air以超声速或高超声速流过钝头体时，激波后的高温会引起混合气体在短距离内发生剧烈的化学反应、释放大量能量，这就是激波诱导燃烧流动，该流动对数值计算方法的精度提出了更高的要求。Lehr等人在20世纪70年代针对激波诱导燃烧现象做了大量试验^[7]，得到激波与燃烧波位置、形状的数据，Soetrisno等人采用TVD格式对相应工况下激波诱导燃烧现象进行了数值模拟^[10]。这些试验数据与数值模拟结果已经成为测试化学反应流模拟能力的常用参考算例。

采用的计算模型与Lehr试验模型一致，即半径7.5 mm的球头，计算网格尺度为120×90×120，网格信息按照非结构网格处理；混合气体来流的组元摩尔分数比为：H₂∶O₂∶N₂=2∶1∶3.76。

3.1 验证算例

当来流马赫数为4.512，静温250 K，静压42 662 Pa时，Lehr通过试验观察到激波和燃烧阵面解耦的爆轰流场，为验证Roe格式在多组元燃烧流场数值模拟中的可行性，采用Shang提出的7组元7反应模型^[15]模拟了该条件下的爆轰流场。图 1给出了流场中温度等值线分布。与试验观察到的现象一致，激波阵面与燃烧波阵面之间存在一个诱导区域。这是由于流场在经过激波后的压强和温度未达到自动点火的条件，预混气体发生燃烧反应的诱导期增加。图 2(a，b)分别给出了驻点线上压强、温度分布以及H₂,O₂,H₂O质量分数分布。从图中可以看出，驻点线上H₂质量分数在诱导区内几乎不变，燃烧波后H₂和O₂反应生成H₂O的质量分数逐渐增加。表 1给出了激波脱体距离与试验数据的比较结果，误差在可接受范围之内，认为Roe格式可以扩展应用于三维多组元燃烧流场的数值模拟与分析。

此外，从图 2(a)中可以发现燃烧阵面所在位置压强出现一个尖峰，这是由于燃烧反应剧烈，能量释放快，难以维持等压燃烧所致，而文献^[10]计算结果并未体现这一现象，说明了Roe格式可以精确分辨流场中的细致结构，体现了其高精度、高分辨率的特性。从图 2中可以发现诱导区厚度明显小于文献^[10]数值模拟结果，分析认为是由于受到化学反应机理影响。激波后诱导区厚度由诱导时间决定，不同的化学反应机理由于考虑的基元反应不同而具有不同的诱导时间^[16-17]。为证实这一分析结果，在3.2节详细讨论了化学动力学模型对燃烧流场数值模拟结果的影响。

3.2 化学动力学模型对诱导区厚度的影响分析

为便于比较，采用与3.1节相同的计算网格与计算条件，仅改变化学动力学模型。模型包含3种：Shang的7组元7反应模型，Moretti^[18]的7组元8反应模型以及Balakrishnan-Williams^[19]的9组元21反应模型，为简化描述，称之为模型1,2,3。图 3给出了采用不同化学反应模型时，流场中驻点线上的参数分布。从图中可以看出，模型1在燃烧阵面后温度升高速率明显比其他两种模型小，反映了该反应模型能量释放速率较小，而模型3在燃烧阵面后温度升高最多，反映了其能量释放总量是三者中最大的。表 2给出了采用不同化学反应模型计算出的驻点线区域的激波脱体距离、诱导区厚度等参数，并与Soetrisno的计算结果进行了对比发现：化学反应模型对激波脱体距离的影响较小，模型1~3的误差都在可以接受的范围之内；化学反应模型直接影响诱导区的厚度，模型1诱导区厚度的误差较大。模型1与模型2,3的区别在于模型1考虑了反应H₂+O₂=OH+OH，这个反应对诱导时间的影响很大，直接减小诱导时间，进而导致诱导区厚度减小。模型3在模型1,2所包含组元基础上，额外考虑了HO₂和H₂O₂两种中间产物，并添加了相应的基元反应，对燃烧过程描述更为细致，但计算量是三者中最大的。

4 结论

以Roe提出的通量差分分裂格式为基础，理论推导了其扩展应用于多组元燃烧流场时对流通量离散的具体形式。采用经典流动问题对数值方法和计算程序进行了测试，比较了不同化学反应模型对燃烧流场结构的影响，并通过数值计算结果与试验数据的对比，得到以下结论：

(1) 通量差分分裂格式Roe方法可以扩展应用于三维多组元燃烧流场的数值模拟与分析；

(2) 扩展后的Roe格式能保持耗散低、分辨率高、激波捕捉能力的特点，可以很好地捕捉和描述多组元化学反应流的流场结构；

(3) 化学动力学模型的选择直接影响激波诱导燃烧流场诱导区厚度大小，但对激波位置影响较小。针对实际问题合理选择化学动力学模型是准确模拟燃烧流场的关键。对于超声速条件下诱导燃烧流场，权衡计算量与计算精度之间的关系，Moretti的7组元8反应模型能够得到较为理想的结果。