动力学格式在超声速非平衡流数值模拟中的应用

卓长飞，武晓松，封 锋

（南京理工大学 机械工程学院，南京210094）

高超声速飞行器在再入大气层的过程中与大气层相互作用，在其头部周围形成一道强弓形激波，激波后高温气体发生化学反应，影响到飞行器周围流场的气动特性和物面热交换。再入飞行器轨道上的大部分区域，流场处于化学非平衡状态。因此，对于高超钝头体飞行器模拟主要采用有限速率化学非平衡流模拟。国内外在超声速化学非平衡流方面的研究已经取得一定进展，发展了一些高效、高精度的数值方法。

在计算流体力学中常见的求解方法主要分为两大类［1］：Godunov方法和Boltzmann方法。Godunov方法是基于黎曼解构造通量；Boltzmann方法是基于微观粒子分布函数构造通量，这种基于气体分子动力学理论求解可压缩流动的数值格式称为GKS（Gas Kinetic Scheme）。由于气体动力学格式具有健壮性和高效性，在近几十年已经得到广泛的应用，并且在不断发展之中。目前主要有两大类气体动力学格式：KFVS格式（Kinetic Flux－vector Splitting Scheme）和BGK格式（Bhatnagar Gross Krook Scheme）。KFVS格式是基于求解最简单的无碰撞的Boltzmann方程，在采用KFVS格式构造数值通量时需要计算exp指数函数和erf误差函数，因此其计算效率显得略低。Agarwal和Acheson［2］提出了一种由现有的KFVS格式和WPS格式结合而成的高效率高精度的格式Kinetic Wave／Particle Splits Scheme（KWPS格式）。从本质上来说KWPS格式也由气体动力学理论推导而来，因此它保留了KFVS格式的精确性和健壮性，同时在构造数值通量过程中不需要计算exp指数函数和erf误差函数，减小了计算量，保持了WPS格式的高效性。

国内外文献表明气体动力学格式得到广泛的发展和应用［1－8］，但用于化学非平衡流数值模拟的文献较少［4］。为了验证2种动力学格式（KFVS格式和KWPS格式）能否在超声速化学非平衡流场数值模拟中仍然表现出高精度、高分辨率特性，同时验证这2种气体动力学格式能否应用于非结构网格中，本文拟将在结构网格提出的2种气体动力学格式（KFVS格式和KWPS格式），结合物理量线性重构技术和限制器，得到基于非结构网格高精度、高分辨率的动力学通量分裂格式有限体积方法，并推广到超声速化学非平衡流场数值模拟中，为计算超声速化学非平衡流提供一种新方法。

1 数值方法

1.1 控制方程

在笛卡尔坐标系下微分守恒形式的二维轴对称化学非平衡流Euler方程为

式中：Q为守恒变量；F，G分别为轴向和径向无粘通量；S1为轴对称源项；S0为化学反应源项。这里仅列出主要的变量：

式中：ρ为密度；u，v分别为轴向、径向速度；p为压力；E为单位体积的总能；wl（l＝1，…，N）为l组分质量分数，ωl为l组分生成速率，N为总组分数。

1.2 动力学格式

本文采用格心格式的有限体积法，因此首先对微分方程进行体积分，再应用散度定理得到积分形式的控制方程：

式中：Wk，Ak分别为控制体单元第k个边界的通量、面积；Nf为控制体单元的边界个数。

计算控制体单元第k个边界对流通量Wk的KFVS格式表达式为

式中：

计算控制体单元第k个边界对流通量Wk的KWPS格式表达式如下：

式中：

显然，这2种动力学矢通量分裂格式不包含自由参数，而且比常见的FVS类、FDS类和AUSM类格式形式简单，易于编程。

1.3 物理量重构

在计算单元边界面通量时，直接取边界两边单元的物理量，在空间上只有一阶精度，数值耗散较大。为了得到高阶精度，需要进行空间高阶插值。本文利用变量线性重构技术［9］获得空间二阶精度。物理量重构表达式为

式中：U表示守恒变量；i，j分别为左、右单元中心点编号；r为两点位置矢量；守恒变量的梯度通过格林－高斯公式计算：

式中：Vi为i单元的体积；nij为i，j单元交界面的单位法矢量；Aij为i，j单元交界面的面积。

1.4 限制器

在跨声速和超声速计算中，为了抑制激波附近可能出现的非物理振荡，保持数值格式的单调性，需要使用限制器。限制器通过对物理量的梯度进行限制，从而避免在重构过程中出现新的极值。本文采用收敛性较好的Venkatarishnan限制器，具体表达式参见文献［9］，这里不再列出。使用限制器后，物理量重构公式可以写作：

式中：φi，φj分别表示i，j单元内的限制器。

1.5 化学反应模型

化学反应模型是决定成功模拟化学非平衡流的关键之一。在本文的验证算例中，氢气和氧气的化学反应机理采用7组分8步基元反应模型［10］，氮气和氧气的化学反应机理采用高温空气5组分12反应模型［11］。详细反应表达式与参数如表1和表2所示。表中：阿乌累里斯公式表达式为Kf＝ATbexp（－E0R0－1T－1）；指 前 因 子A的 单 位为（cm3／mole）m－1·s－1，其中m为基元反应级数，b为温度指数，E0为反应活化能，R0为气体常数，M表示第三碰撞体。

表1 H2＋O2基元反应模型

表2 N2＋O2基元反应模型

1.6 时间算子分裂算法

化学非平衡流控制方程组分为流动部分和化学反应部分，两者相互耦合。本文采用时间算子分裂算法处理这种耦合过程［12］，即把反应流方程组分为流动和化学反应2部分，如式（8）和式（9）所示，并把求解流动偏微分方程时采用的时间步长进一步细分，作为求解化学反应刚性常微分方程的步长，计算化学反应对流场的贡献。具体做法是先冻结化学反应求解得到流场参数，然后将化学反应看做等容放热或吸热过程，保持内能、速度参数不变，计算各组分的质量变化率，最后迭代求解温度得到流场其他参数。

2 计算结果与分析

2.1 H2＋O2激波诱导燃烧

当按理想化学当量比的混合氢气／氧气或者氢气／空气以超声速或高超声速流过钝头体时，会产生较强的弓形激波，波后温度升高会引起混合气体燃烧，这就对所采用的数值计算方法精度提出了更高要求。Lehr对此进行了大量实验研究，获得激波形状、激波位置的数据。研究人员又进行数值模拟，得到与实验比较吻合的结果。本文选用其中2个工况验证本文数值方法的可靠性。

自由来流条件：来流马赫数5.08，静温291.5K，来流速度2 705.0m／s，静压24 797Pa，氢气和氧气的体积混合比为2∶1。球头半径R＝0.007 5m，参考的实验数据取自文献［13］，参考的数值模拟数据取自文献［14－15］。非结构网格节点总数为12 409，单元边总数为36 578，单元总数为24 170。

图1给出了2种动力学格式计算得到的密度等值线。激波与燃烧波耦合形成爆震波，激波层厚度略有增加，2种格式得到的激波面位置均与实验数据［13］基本吻合。图2～图4分别给出了2种格式计算得到的驻点线（无量纲距离X＊：驻点线上距驻点距离与球头半径之比）上无量纲密度ρ＊（驻点线上密度与来流密度之比）、无量纲压力p＊（驻点线上静压与来流静压之比）与无量纲温度T＊（驻点线上静温与来流静温之比）和主要组分质量分数w分布，所有计算结果与参考文献［14－15］相比基本吻合。其中，本文计算得到的密度和压力与参考文献计算的结果相比，受激波与燃烧波耦合影响而产生的Von Neumann尖峰更为陡峭，说明了本文的数值模拟更加精确，同时说明了动力学格式在超声速化学非平衡流中高精度、高分辨率特性，精确地分辨了流场中的细致结构。

图1 H2＋O2钝体绕流场密度等值线图

图2 H2＋O2钝体绕流场驻点线上无量纲密度分布

图3 H2＋O2钝体绕流场驻点线上无量纲压力和无量纲温度分布

图4 H2＋O2钝体绕流场驻点线上主要组分质量分数分布

2.2 H2＋Air激波诱导燃烧

自由来流条件：来流马赫数6.46，静温286.6K，来流速度2 605.0m／s，静压42 662Pa，氢气、氧气、氮气的体积混合比为2∶1∶3.76，在本算例的波后温度条件下可不计N2与O2之间的反应。球头半径R＝0.007 5m，参考的实验数据取自文献［13］，参考的数值模拟数据取自文献［14－15］。非结构网格节点总数为9 663，单元边总数为28 440，单元总数为18 778。

图5给出了2种动力学格式计算得到的密度等值线。由于N2的稀释作用，燃烧释放出的热量小于自由来流动能，此时虽然激波与燃烧波重合，但整个流场中不存在较强的爆震波，2种格式得到的激波面位置均与实验数据［13］基本吻合。图6、图7和图8分别给出了2种格式计算得到的驻点线上无量纲密度ρ＊、无量纲压力p＊与无量纲温度T＊和主要组分质量分数w分布。由于N2的稀释作用，压力不再出现Von Neumann尖峰，与参考文献［14－15］的计算结果相比，本文计算得到的Von Neumann尖峰仍然更为陡峭，说明了动力学格式在该算例中精确描述了激波与燃烧波耦合产生的陡峭的Von Neumann尖峰。需要说明的是，国内多篇文献采用不同数值方法模拟球头激波燃烧算例均采用文献［14－15］作为参考，计算得到的Von Neumann尖峰均比参考文献陡峭。这可能是参考文献年代比较久远，限于当时数值模拟的水平，数值模拟结果没能较好地描述Von Neumann尖峰。

图5 H2＋Air钝体绕流场密度等值线图

图6 H2＋Air钝体绕流场驻点线上无量纲密度分布

图7 H2＋Air钝体绕流场驻点线上无量纲压力和无量纲温度分布

图8 H2＋Air钝体绕流场驻点线上主要组分质量分数分布

2.3 N2＋O2球锥高超声速非平衡流场

KIM［16］推导了AUSMPW＋格式并用于数值模拟高超声速空气通过球锥的非平衡流场，与实验结果对比表明AUSMPW＋格式在高超非平衡流数值模拟中表现效果较好。本文采用该算例作为验证数值方法的可靠性算例之一，化学动力学模型采用高温空气5组分12反应模型。

自由来流条件：来流速度3 630m／s，静温293K，静压2 400Pa。球锥半径R＝0.007m，参考的实验数据取自文献［16］。非结构网格节点总数为16 145，单元边总数为48 001，单元总数为31 856。文献［16］的数值模拟结果中没有给出驻点线上参数分布，因此本文给出驻点线上主要参数分布。

图9中给出了计算得到的密度、压力等值线，其中每部分的上半部分为KFVS格式计算结果，下半部分为KWPS格式计算结果。2种格式计算得到的激波面位置均与文献［16］给出的实验数据吻合。此外，本算例的来流马赫数约为11.0，与文献［16］采用AUSMPW＋格式计算结果相同，并且不会出现文献中提及的在高马赫数下计算产生的carbuncle现象。

图10和图11分别给出了驻点线上2种格式计算得到的驻点线上的无量纲温度和主要组分分布。高速空气通过球锥产生激波，波后温度迅速升高，N2和O2开始吸收热量并发生离解和置换反应，由于N2与O2反应过程吸收热量，因此沿驻点线上温度逐渐降低，N2和O2组分质量分数也均降低，NO组分质量分数升高。在该算例温度下，N2离解较少，N2质量分数降低的主要原因是发生置换反应生成NO。而O2质量分数降低的主要原因一方面是O2离解，另一方面是生成NO。结合图10和图11可以看出，当驻点线轴线温度低于N2离解温度时（本算例条件下离解起始温度为4 000K左右），N2不再离解，仅发生置换反应，但O2的离解和置换反应仍会进行。因此，在趋近于壁面处温度和组分质量分数变化均减缓。

图9 空气高超声速钝体绕流场密度和压力的等值线图

图10 空气高超声速钝体绕流场驻点线上无量纲温度分布

图11 空气高超声速钝体绕流场驻点线上主要组分质量分数分布

从以上3个算例看出，KFVS和KWPS格式计算的结果基本重合，并且均与参考文献数据基本吻合，说明这2种动力学格式在超声速化学非平衡流场数值模拟中均能较好地分辨出激波和燃烧波等物理现象，流场结构清晰，在强间断处基本无振荡，没有出现非物理解，具有较高精度和分辨率。

3 结论

本文将2种基于气体分子动力学理论的动力学格式（KFVS格式和KWPS格式）用于基于非结构网格超声速化学非平衡流场数值模拟中，并且成功模拟了几个经典的超声速化学非平衡流算例。针对本文数值模拟的算例，可以得到以下结论：

①将原来在结构网格中提出的2种动力学通量格式（KFVS格式和KWPS格式）推广应用于非结构网格中，数值结果表明该方法是可行的。

②将2种动力学格式（KFVS格式和KWPS格式）应用到超声速化学非平衡流场的数值模拟中，数值模拟结果与前人实验数据和数值模拟结果进行了比较，说明了本文将2种动力学格式推广到超声速化学非平衡流场数值模拟中是可行的，并且2种动力学格式在超声速化学非平衡流场数值模拟中均保持高精度、高分辨率特性，为计算超声速化学非平衡流提供一种新方法。

③2种动力学格式不含自由参数，形式简单，同时能够正确地分辨流场中的激波、燃烧波和其他物理现象，流场结构清晰，具有较高精度和高分辨率；2种动力学格式在超声速化学非平衡流场中计算精度和分辨率基本一致。

［1］LI S H，XU K.Entropy analysis of kinetic flux vector splitting schemes for the compressible Euler equations［J］.Zeitschrift f¨ur Angewandte Mathematik und Physik，2001，52：62－78.

［2］ACHESON K E，AGARWAL R K.A kinetic theory－based wave／particle flux－splitting scheme for the Euler equations.AIAA95－2178［R］.1995.

［3］REKSOPRODJO H S R，AGARWAL R.Implicit kinetic schemes for Euler equations，AIAA2001－2629［R］.2001.

［4］EPPARD W M，GROSSMAN B.An upwind kinetic flux－vector splitting method for flows chemical and thermal non－equilibrium，AIAA93－0894［R］.1993.

［5］RAVICHANDRAN K S.Higher order KFVS algorithms using compact upwind difference operators［J］.Journal of Computational Physics，1997，130：161－173.

［6］MANDAL J C，JAIN K.A new implicit formulation of KFVSscheme for Euler equations，AIAA2006－3 709［R］.2006.

［7］LEI T.Progress in gas－kinetic upwind schemes for the solution of Euler／Navier－Stokes equations－I：overview［J］.Computer ＆Fluids，2012，56：39－48.

［8］汤华中，邬华谟.高分辨率KFVS有限体积方法及其CFD应用［J］.计算数学，1999，21（3）：375－384.TANG Hua－zhong，WU Hua－mo.High resolution KFVS finite volume methods and its appliciton in CFD［J］.Mathematica Numerica Sinica，1999，21（3）：375－384.（in Chinese）

［9］VENKATARISHNAN V.On the accuracy of limiters and convergence to steady state solutions，AIAA93－0880［R］.1993.

［10］刘君，周松柏，徐春光.超声速流动中燃烧现象的数值模拟方法及应用［M］.长沙：国防科技大学出版社，2008：288－289.LIU Jun，ZHOU Song－bai，XU Chun－guang.The application and method of combustion phenomena in supersonic flow［M］.Changsha：National University of Defense Technology Press，2008：288－289.（in Chinese）

［11］PARK C，YOON S.Caluclation of real－gas effects on blunt－body trim angles［J］.AIAA Joural，1992，30（4）：999－1 007.

［12］胡湘渝，张德良，姜宗林.气相爆轰基元反应模型数值模拟［J］.空气动力学报，2003，21（1）：59－66.HU Xiang－yu，ZHANG De－liang，JIANG Zong－lin.Numerical simulation of gaseous detonation with detailed chemical reaction model［J］.Acta Aerodynamica Sinica，2003，21（1）：59－66.（in Chinese）

［13］LEHR H F.Experiments on shock－induced combustion［J］.Aeronautica Acta，1972，17：589－597.

［14］SOETRISNO M，IMLAY S T.Simulation of the flow field of a ram aceelerator，AIAA91－1915［R］.1991.

［15］YUNGSTER S，EBERHARDT S，BRUCKNER A P.Numerical simulation of hypervelocity projectiles in detonable gases［J］.AIAA Journal，1991，29（2）：187－192.

［16］KIM K H.Accurate computation of hypersonic flows using AUS－MPW ＋ scheme and shock－induced grid technique，AIAA98－2442［R］.1998.