吴 伟，许厚谦，王 亮，薛 锐

（南京理工大学 能源与动力工程学院，南京210094）

传统的数值方法主要包括有限元、有限差分、有限体积等，上述方法原理各不相同，但都需要采用网格划分求解域，因而其应用受到了网格的制约，如结构网格处理相对复杂外形的求解区域能力有限；非结构网格虽然对于复杂外形的处理灵活方便，但对于变形、动边界及不连续等复杂情形，生成高质量的网格可能花费与计算相当的时间；而笛卡尔网格虽然方便、快速，但是对于边界的处理却较为复杂、低效。为了摆脱网格的束缚，Batina［1］首先将无网格方法引入到计算流体力学领域，采用一系列点云离散求解域。无网格方法无需将节点连成网格，布点快捷，无需人为干涉，对于复杂外形更具灵活性。近十几年，无网格方法受到国内外学者的关注、研究，其内容主要包含自动布点生成点云，发展适用无网格方法的高激波分辨率、高精度数值格式，处理非定常问题的无网格方法以及混合网格方法等［2－5］，然而，现有的无网格方法多见应用于空气动力学领域，对于航空航天、兵器等领域中广泛存在的复杂化学反应流场的模拟还显得无能为力，对上述问题展开深入研究具有重要的意义和应用价值。

本文基于PC机集群，研究了耦合化学反应计算的并行无网格方法，用于包含复杂激波、燃烧波的非平衡化学反应流的数值模拟。其基本变量解函数采用线性函数描述，空间导数采用最小二乘拟合逼近，将多组分HLLC格式引入无网格方法用于数值通量的计算，时间项采用4阶Runge－Kutta法显式推进，反应源项采用有限速率反应模型处理，进程间的数据通信采用MPI。采用上述无网格方法对激波诱导燃烧流场以及高速运动弹丸诱导斜爆轰流场进行了数值模拟。

1 控制方程

本文采用含化学反应源项的二维多组分Euler方程，具体形式如下：

守恒变量U，对流通量F、G，化学反应源项W，轴对称源项S的定义如下：

式中：t为流动时间；ρs为组分s的质量密度；N为组分总数；ρ为混合气体的质量密度；u，v分别为x，y方向速度分量；p为混合气体总压；ωs为化学反应引起组分s的质量生成率；ρE为单位体积总能。

式中：hs为组分s的焓值。各组分的焓值是温度的函数，可按下多项式拟合获得：

式中：Ru为通用气体常数；Ms为组分s的摩尔质量；T为混合气体温度；系数a1～a6可由JANAF表［6］查得。反应中认为混合气体及各组分气体满足理想气体状态方程：

式中：M为混合气体的平均摩尔质量。

2 数值方法

2.1 对流通量计算

本文假设流动解函数满足线性关系，通过最小二乘拟合逼近空间导数。由此可得任意中心点i的对流通量为

形函数bij、cij以及空间导数具体求解过程参见文献［5］，其中，中心点i及其卫星点j的中点的数值通量Wij采用多组分HLLC格式计算，具体如下所示：

式中：下标K根据式（7）取i或j，Si、Sj分别为中心点i、中心点j的波速，SM为接触间断面速度，lx、ly分别为单位法向量在x、y方向的分量，unK为中心点i或j处流体速度的法向分量，p＊按下式计算：

具体计算参见文献［4－5］，在此不再赘述。

2.2 化学反应模型

化学反应源项采用有限速率反应模型计算，反应体系中的任意反应可以表示为

式中：ν′sm、ν″sm为组分s在反应m中的化学反应当量系数，χs为组分s的化学符号，正向反应速率常数Kfm可由Arrhenius公式计算获得：

式中：Am为指前因子，bm为温度因子，Em为活化能。逆向反应速率常数Kbm可由反应m的平衡常数或Arrhenius公式计算获得。任意组分s的化学反应源项可由下式计算：

式中：NR为化学反应总数。

2.3 时间离散

对于任意中心点i，式（1）可写成如下半离散形式：

式中：R为点i的残差矢量。对式（14）可采用4阶Runge－Kutta法进行显式时间推进，形式如下：

式中：l＝1，2，3，4，φl分别为1／4、1／3、1／2、1。为了保证计算的稳定，采用各进程中最小时间步长作为全局时间步长Δt。流动与化学反应解耦计算，为克服“刚性”问题，将化学反应时间步长作进一步细分处理。

3 并行处理

本文并行计算在多台通过D－Link集线器互连的Dell台式机上进行，进程之间的数据通信均采用具有良好可移植性、可扩展性的MPI。并行计算中负载平衡是影响计算效率的重要因素。合理的分配策略是保证并行效率的关键。为了保证各进程分配均匀，本文采用相对简单的一维划分，将离散点进行分区编号，根据进程数平均分配，各进程保存各自分区的离散点，并生成相应的分区通信边界，建立向相邻进程发送、接受数据离散点链表。计算每推进1个时间步长，各进程向相邻进程发送并从相邻进程接受通信边界节点数据。

4 算例

共选取2个算例验证本文算法的有效性，分别为激波诱导燃烧和高速飞行弹丸诱导爆轰实验，算例1是半径R为7.5mm的球－柱体钝头弹射入等化学当量比的H2／Air预混气体，算例2为尖头弹（简化结构如图1所示）射入等化学当量比的CH4／Air预混气体，压力p0、温度T0、速度v0如表1所示。氢气／空气采用7组分8反应模型，正／逆反应速率常数计算系数取自文献［7］，甲烷／空气采用13组分19反应模型，正向反应速率常数计算系数同文献［8］，逆向反应速率由平衡常数计算。本文计算离散点通过非结构网格无网格化方法获得。

图1 驻定斜爆轰波实验弹丸结构

表1 实验条件

算例1采用非均匀分布的121 226个节点离散流场，其中驻点流线上布点140个。图2、图3分别为采用4进程计算得到的等温线和H2O质量分数分布，图2中“■”为实验照片中激波阵面的位置，可见自由来流速度大于临界爆轰速度，激波阵面与燃烧阵面耦合。图4、图5为沿驻点流线无量纲压力和温度的分布曲线，可见激波后温度、压力迅速上升并达到平衡状态，同文献［9］采用网格方法得到的结果吻合较好。此外本文算法成功捕捉到压力的VonNeumann尖峰。

图2 等温线分布

图3 H2O质量分数分布

图4 沿驻点流线无量纲压力分布

图5 沿驻点流线无量纲温度分布

算例2采用非均匀分布的207 532个节点离散，对弹丸表面附近的节点进行了加密处理。图6为采用3进程计算得到的密度云图，作为比较，同样给出了实验阴影照片。

图7上下分别为无网格方法与非结构网格［8］方法计算得到的温度云图，可见弹丸头部激波较弱，不足以点燃可燃混合气体，仅形成驻定的斜激波，计算得到激波角为35.4°，同实验获得的37.8°以及非结构网格计算得到的36.5°吻合较好。图8为CH4以及CO2质量分数分布云图，可见仅在弹丸中部凹槽内以及弹尾发生了化学反应。

图6 密度云图与实验阴影照片对比

图7 温度云图对比

图8 CH4和CO2质量分数分布

5 效率分析

本文分别使用2～8个进程对上述2个算例进行了并行计算。由于并行程序同串行程序存在差异，将两者进行比较不够合理，因而本文采用2进程计算耗时为参考量计算加速比δ、并行效率η，其部分结果见表2。可见随着进程个数NP的增加，并行效率逐渐下降，这是由于NP的增加导致各进程分配有效计算节点的减少，使得计算与通信的耗时比减小。此外虽然算例2离散点个数多于算例1，但是需要计算化学反应的节点却远少于算例1，其计算时间与通信时间相比所占比重较低，因此算例2的并行效率均低于算例1。

表2 并行加速比及并行效率

6 结束语

本文基于MPI并行坏境，发展了非平衡化学反应流数值模拟的并行无网格算法，并对激波诱导燃烧和高速飞行弹丸诱导爆轰流场进行了数值模拟，其结果同实验以及网格方法结果吻合较好，并且获得了较为理想的并行效率，表明了计算形式简单的多组分HLLC格式具有较高的激波分辨率，验证了本文并行无网格算法的正确性。本文研究不仅拓展了无网格方法处理复杂流场的能力，丰富了非平衡反应流模拟的数值方法，而且为航空航天、兵器等领域中众多化学反应流的大规模数值计算提供了有效的工具。

［1］BATINA J T.A gridless Euler／Navier－Stokes solution algorithm for complex－aircraft applications，AIAA 93－0333［R］.1993.

［2］孙迎丹，王刚，叶正寅.AUSM＋－up格式在无网格算法中的推广［J］.空气动力学学报，2005，23（4）：511－515.SUN Ying－dan，WANG Gang，YE Zheng－yin.Extending the AUSM＋－up scheme into the gridless method［J］.Acta Aerodynamica Sinica，2005，23（4）：511－515.（in Chinese）

［3］WANG H，CHEN H Q，PERIAUX J.A study of gridless method with dynamic clouds of points for solving unsteady CFD problems in aerodynamics［J］.Int J Meth Fluids，2010，64：98－118.

［4］周星，许厚谦.无网格算法在膛口流场数值模拟中的应用［J］.弹道学报，2011，23（1）：23－26.ZHOU Xing，XU Hou－qian.Application of gridless method in numerical simulation of muzzle flow field［J］.Journal of ballistics，2011，23（1）：23－26.（in Chinese）

［5］ZHOU X，XU H Q.Gridless method for unsteady flows involving moving discrete points and its applications［J］.Engineering Applications of Computational Fluid Mechanics，2010，4（2）：276－286.

［6］STULL D R，PROPHET H.JANAF thermochemical tables［J］.Journal of Physical and Chemical Reference Data，1975，4（1）：1－176.

［7］EVANS J S，SCHEXNAYDER C J.Influence of chemical kinetics and unmixedness on burning in supersonic hydrogen flames［J］.AIAA，1980，18（2）：188－193.

［8］代淑兰，许厚谦.驻定斜爆轰波并行数值模拟［J］.爆炸与冲击，2006，26（6）：572－576.DAI Shu－lan，XU Hou－qian.Numerical simulation of standing oblique detonation wave using parallel computation method［J］.Explosion and Shock Waves，2006，26（6）：572－576.（in Chinese）

［9］孙明波，梁剑寒，王振国.非平衡流解耦方法及其计算激波诱导燃烧的应用验证［J］.航空动力学报，2008，23（11）：2 055－2 061.SUN Ming－bo，LIANG Jian－han，WANG Zhen－guo.Validation of an uncoupled solver of non－equilibrium flow for shock－induced combustion［J］.Journal of Aerospace Power，2008，23（11）：2 055－2 061.（in Chinese）