多组份流动质量分数保极值原理算法

唐维军， 蒋 浪， 程军波

（1.北京应用物理与计算数学研究所，计算物理实验室，北京 100088；2.中国工程物理研究院研究生部，北京 100088）

多组份流动质量分数保极值原理算法

唐维军1， 蒋 浪2， 程军波1

（1.北京应用物理与计算数学研究所，计算物理实验室，北京 100088；2.中国工程物理研究院研究生部，北京 100088）

对基于质量分数的Mie-Gruneisen状态方程多流体组份模型提出了新的数值方法.该模型保持混合流体的质量、动量、和能量守恒，保持各组份分质量守恒，在多流体组份界面处保持压力和速度一致.该模型是拟守恒型方程系统.对该模型系统的离散采用波传播算法.与直接对模型中所有守恒方程采用相同算法不同的是，在处理分介质质量守恒方程时，对波传播算法进行了修正，使之满足质量分数保极值原理.而不作修改的算法则不能保证质量分数在［0，1］范围.数值实验验证了该方法有效.

Mie-Gruneisen状态方程；质量分数保极值原理；波传播算法

0 引言

可压缩多介质复杂流动有很广泛的应用，包括爆轰，燃烧，惯性约束聚变，天体物理等诸多领域.研究这些复杂流动的目的是要了解多介质流动的流体动力学机制.这涉及到多介质流动的数学建模及相关模型的数值离散方法.一般而言，如不考虑粘性、表面张力、热力学效应等物理因素，单纯使用欧拉方程描述单介质流体无粘流动时，其数学模型由质量、动量、能量的三个守恒律组成的完全守恒系统构成，且该完全守恒系统可用近三十年来发展成熟的各种守恒格式进行数值离散.但如流动涉及多介质时，仅用完全守恒模型来描述却难以实现.在一些非常特殊的情形，有些作者尝试用完全守恒系统来描述多介质流动［1－3］，但涉及混合流动时，相应的数值格式处理则变得相对复杂；并且在计算纯界面问题时，界面附近压力或速度的数值解可能会出现非物理振荡.目前比较公认的是采用拟守恒模型来描述多介质复杂流动.例如对理想气体状态方程一维两流体组份流动，利用对纯接触问题界面压力和速度数值解要求保持一致这种思想，Abgrall［4］最早揭示模型可由混合介质质量、动量和能量守恒方程，以及与混合状态方程参数有关的量的一个对流方程来描述.并且Abgrall［4］还确认该对流方程的离散格式不是独立的，而与三个守恒方程的离散格式有关.这个事实显示多介质流动的数学模型和数值格式与单介质流动情形相比要复杂.随后Abgrall［4］思想被推广到刚性气体状态方程［5－7］、van der Waals气体状态方程［8］、密实介质状态方程［9］、以及Mie-Gruneisen状态方程［10－12］.此外，还有众多的其它模型来描述多介质流动，如加入level set方程的扩展欧拉方程组模型［13］，Karni的原始变量模型［14］和杂交模型［15］，Allaire等的5方程模型［16］，Saurel和Abgrall的两相流7方程模型［17］，使用体积分数的4方程模型［18－19］，反扩散界面模型等［20－22］，两相流BGK模型［23］.至于有关上述模型的数值格式，散见上述各文及相关引文，此处不一一列举.

讨论Mie-Gruneisen状态方程一维模型的数值离散格式.该模型基于Abgrall［4］等效方程思想，且分介质质量保持守恒的模型［12］的数值离散格式.与基于体积分数对流方程模型［10］相比，本文所用模型［23］不仅能保持混合流体的质量、动量、和能量守恒，还保持各组份分质量守恒，且对纯接触问题保持界面压力速度一致.本文使用波传播算法［7，10，24］离散系统，但不是象文［12］对所有守恒方程进行相同数值离散，而是在对分介质质量守恒方程离散时，采用Larrouturou［1］对完全守恒系统设计的质量分数保极值原理算法，对波传播算法进行修改.即对分介质质量守恒方程的数值通量进行修改，同时引进一个类似CFL条件的时间步长约束.其目的是确保守恒系统离散时，各介质质量分数始终在［0，1］范围.如对拟守恒系统直接采用单介质流动常用的那些数值格式，如TVD格式，MUSCL格式，PPM格式等，都无法确保质量分数在［0，1］范围.Abgrall［4］在用Roe和MUSCL格式处理理想气体状态方程两组份完全守恒模型数值离散时，借用所谓压力速度界面一致条件来导出状态方程混合参量的某个特殊离散格式，在这种情形质量分数满足保极值原理.但这种方法难以推广到本文讨论的Mie-Gruneisen状态方程情形.特别值得提到的是，文［2］也采用了Larrouturou［1］通量修正格式，但文［2］模型与本文模型不同，且仅限理想气体状态方程，无法推广到其它状态方程.另外，本文对波传播算法二阶格式的通量修改格式，与Larrouturou［1］二阶MUSCL的通量修改格式不一样，后者对质量分数的重构需加额外的约束.本文算例表明，这种数值格式的有效性.

1 模型方程

文［12］的模型不仅能保持混合流体的质量、动量、和能量守恒，还保持各组份分质量守恒，保留了Shyue［5］的等效方程，整个模型是一个拟守恒系统.不妨设整个流场有2种组份，其模型系统为

其中ρ，u，p，E分别表示密度（kg·m－3），速度（m·s－1），压力（Pa），比总能量，Y表示第一个组份质量分数.Mie-Gruneisen状态方程为

和内能.将方程组（1）－（7）写成拟线性形式

其中q＝（ρ，ρu，ρE，ρY，1／Γ，pref／Γ，ρeref）T，Jacobia矩阵A

矩阵A对应的特征值分别为

特征值向量Λ对应的右特征向量记为R＝（r1…，r7）．这些特征值和右特征向量满足Ark＝λkrk，k＝1，…，7．

Jacobia矩阵A对应前4个守恒方程的子矩阵记为AC，而对应前4个守恒方程的守恒通量

2 一维Riemann逼近问题

考虑一维Riemann逼近问题

其中，qL＝qi，qR＝qi＋1．设矩阵AH是Jacobia矩阵A在Roe平均状态qH的局部线性化，即AH＝A（qH）．Roe线性化要满足如下两式

上述过程中，平均量（1／Γ）H，（p／Γ）H须同式（14）的选取方式，且满足如下逼近（Shyue［14］）

则平均量pH，ΓH按下式计算

要完成AH，还须补充φH，χH，ψH．但注意式（12）方程个数与而式（11）方程个数并不相等，没有其它条件可用，故直接仿式（14）Roe平均给出，

线性化矩阵AH的特征值为

其中Roe平均声速为

这里KH＝

3 数值格式

基于波传播算法［7，10，24］来计算拟守恒系统（9），波传播算法的标准有限体积方法是Godunov型迎风格式.但本文在处理分介质质量守恒方程，不采用标准波传播算法，而采用Larrouturou［1］提出的质量分数保极值原理算法，对拟守恒系统（9）的第4个守恒方程的数值通量进行修改.Larrouturou［1］提出的质量分数保极值原理算法原用来处理理想气体状态方程两组份完全守恒系统，其特点是采用数值格式时，能保持质量分数保极值原理，即完全守恒系统的分介质质量分数始终在［0，1］范围.本文将Larrouturou的这个算法推广到Mie-Gruneisen方程的拟守恒系统（9），不仅确保纯接触问题界面压力速度保持一致，这是Larrouturou［1］算法不曾考虑的，而且能确保拟守恒系统（9）的分介质质量分数满足保极值原理，这是原来波传播算法［7，10，24］处理分介质守恒方程所不具备的.

3.1 波传播算法

波传播算法［7，10］一阶格式为

其中

波传播算法［7，10］二阶格式在一阶格式基础上，增加了采用波限制器的二阶校正项，即

3.2 分介质质量守恒方程通量修正格式

上面已提到，计算拟守恒系统（9）的分介质质量守恒方程，即系统第4个方程

这里

考虑到质量分数Y＝q（4）／q（1），如直接采用波传播算法（20）或（21）计算守恒变量q（1），q（4），则无法确保Y∈［0，1］．故方程（23）的离散采用Larrouturou［1］提出的质量分数保极值原理算法进行通量修改.

虽然模型（9）整体是拟守恒方程组，波传播算法离散（9）的守恒方程部分可写为守恒通量形式.第一个方程，即混合质量守恒方程离散为

其中当波传播算法为一阶格式（20）时，通量分量（F（1））i＋1／2为

当波传播算法为二阶格式（21）时，通量分量（F（1））i＋1／2为

利用Larrouturou质量分数保极值原理算法［1］，分介质质量守恒方程（23）的离散格式为

特别需要提到的是，Larrouturou［1］采用二阶MUSCL格式进行通量修正时，其通量分量修正与上式并不一致，即

且必须对质量分数重构加上如下限制

3.3 时间步长约束

要满足质量分数保极值原理，波传播算法的上述数值格式修改，还需考虑一个类似CFL条件的时间步长约束［1］，即

当波传播算法为一阶格式（20），流通量分量F（1）为一阶精度时，

当波传播算法为二阶格式（21），流通量分量F（1）为二阶精度时，

由此得整个系统（9）采用波传播算法和质量分数保极值原理的时间步长为从后面数值算例实际情况来看，大部分情况下，ΔtMaxP大于ΔtCFL，因此上述时间步长限制基本没有增加计算量.当然，这个实际经验取决于式（33）中的CFL系数CCFL取值大小.

4 数值算例

本文所有算例来自文献［14］.对拟守恒系统（9）使用两种方法进行数值离散.第一种方法是直接使用波传播算法；第二种方法是在波传播算法基础上，再使用质量分数保极值原理对分介质守恒方程数值通量进行修正.由于波传播算法还分为一阶和二阶格式，后面为简明起见，记算法A1和A2分别为直接采用波传播算法的一阶和二阶格式，而记算法B1和B2分别为采用波传播算法一阶和二级格式，且同时采用质量分数保极值原理算法进行分介质质量通量修正.将采用算法B2，网格数为2 000的计算结果作为参考解.网格剖分数为100，两种方法（包括一阶和二级格式）的计算结果进行了对比，结果显示，两种方法计算得到的物理量差距很小，但质量分数是否保极值原理却差距很大.使用算法B1和B2，网格数分别为100，200，400的加密计算显示了本算法计算的收敛性.在下面5个实际算例计算中，第四个算例的CFL系数CCFL取值为0.75，其它算例的CFL系数CCFL均取值为0.9.

例1 （单组份问题）气体TNT爆炸产物Riemann问题，材料为TNT爆炸产物，参考状态方程使用如下JWL状态方程［10］（34）－（36），材料参数Γ0，V0，e0，A，B，R1，R2具体见文［14］.

图1是t＝12 μs时刻的密度、速度、压力、声速以及质量分数.解的结构包括左向传播的稀疏波，右向传播的接触间断和激波.算法A1和B1，A2和B2的计算结果显示差距甚为细微，但二阶格式比一阶格式要更靠近参考解.图2是两种方法的一阶和二阶格式计算的各时刻质量分数极大值随时间变化的历史曲线.从图中可以看出，直接采用波传播算法的二阶格式A2，其结果不能保持质量分数在［0，1］范围.而算法B1、B2及A1则保持质量分数在［0，1］范围.图3是采用B1和B2对网格数分别为100，200，400的加密计算结果，显示计算结果是收敛的.

例2 （单组份问题）铝块撞击铝块，铝板从右向左撞击半无限长铝板.状态方程采用如下激波状态方程［10］（37）－（39），材料参数ρ0，c0，s，Γ0，α，e0，p0见文［10］.

图4是t＝50 μs时刻的密度、速度、压力、声速以及质量分数.解的结构包括左向和右向传播的两个激波，以及左向传播的接触间断.算法A1和B1，A2和B2的计算结果显示差距甚为细微，但二阶格式比一阶格式要更靠近细网格结果.图5是两种方法的一阶和二阶格式计算的各时刻质量分数极大值随时间变化的历史曲线.从图中可以看出，直接采用波传播算法的二阶格式A2，其结果不能保持质量分数在［0，1］范围.而算法B1、B2及A1则保持质量分数在［0，1］范围.图6是采用B1和B2对网格数分别为100，200，400的加密计算结果，显示计算结果是收敛的.

例3 两组份碰撞算例，初始为标准大气压条件.即压力为p0＝1 atm，T0＝300 K.铜板以速度u＝1 500 m·s－1右行，与固体炸药碰撞.铜和固体炸药的参考状态方程都使用如下CC状态方程［10］（40）－（43），材料参数ρ0，cV，T0，A，B，ε1，ε2，Γ0见文［10］，密度ρ＝ρ0．

初始数据为

图7是t＝85 μs时刻的密度、速度、压力、热内能以及质量分数.图中包含左向和右向传播的两个激波，以及右向传播的接触间断.算法A1和B1，A2和B2的计算结果显示差距甚为细微，但二阶格式比一阶格式要更靠近细网格结果.图8是两种方法的一阶和二阶格式计算的各时刻质量分数极大值随时间变化的历史曲线.从图中可以看出，直接采用波传播算法的二阶格式A2，其结果不能保持质量分数在［0，1］范围.而算法B1、B2及A1则保持质量分数在［0，1］范围.图9是采用B1和B2对网格数分别为100，200，400的加密计算结果，显示计算结果是收敛的.

例4 两组份问题，材料为铜和气体爆轰产物，参考状态方程铜使用CC状态方程（40）－（43），而气体爆轰产物使用JWL状态方程（34）－（36），材料参数见文［14］.

初始条件为

图10是t＝73 μs时刻的密度、速度、压力、热内能以及质量分数.解的结构包含左向传播的稀疏波和右向传播的激波，以及右向传播的接触间断.算法A1和B1，A2和B2的计算结果显示差距甚为细微，但二阶格式比一阶格式要更靠近细网格结果.图11是两种方法的一阶和二阶格式计算的各时刻质量分数极大值随时间变化的历史曲线.从图中可以看出，直接采用波传播算法的一阶格式A1和二阶格式A2，其结果均不能保持质量分数在［0，1］范围.而算法B1、B2则保持质量分数在［0，1］范围.图12是采用B1和B2对网格数分别为100，200，400的加密计算结果，显示计算结果是收敛的.

初始条件为

图13是t＝120 μs时刻的密度、速度、压力、Gruneisen系数Γ以及质量分数.图中包含左向传播的系数波和右向传播的激波和接触间断.算法A1和B1，A2和B2的计算结果显示差距甚为细微，但二阶格式比一阶格式要更靠近细网格结果.图14是两种方法的一阶和二阶格式计算的各时刻质量分数极大值随时间变化的历史曲线.从图中可以看出，直接采用波传播算法的一阶格式A1和二阶格式A2，其结果不能保持质量分数在［0，1］范围.而算法B1、B2则保持质量分数在［0，1］范围.图15是采用B1和B2对网格数分别为100，200，400的加密计算结果，显示计算结果是收敛的.

5 结论

对Mie-Gruneisen状态方程多流体组份流动问题，有别于Shyue［10］采用的基于体积分数拟守恒方程组模型，本文采用基于质量分数的Mie-Gruneisen状态方程多流体组份模型.对该模型进行数值离散时，需要设法保持质量分数在［0，1］区间，否则，所得分介质质量有可能为负.将Larrouturou［1］研究理想气体状态方程完全守恒系统的多组份模型数值离散保持质量分数保极值原理的格式应用到Mie-Gruneisen状态方程多流体组份拟守恒系统.数值结果显示了修改的波传播算法格式能保持质量分数保极值原理.

［1］Larouturou B.How to preserve the mass fraction positive when computing compressible multi-component flows［J］.J Comput Phys，1991，95：59－84.

［2］Ton V.Improved shock-capturing methods for multicomponet and reacting flows［J］.J Comput Phys，1996，128：237－253.

［3］Wang S，Anderson M，et al.A thermaodynamically consistent and fully conservative treatement of contact discontinuities for compressible multicomponet flows［J］.J Comput Phys，2004，195：528－559.

［4］Abgrall R.How to prevent pressure oscillations in multicomponent flows：A quasi conservative approach［J］.J Comput Phys，1996，125：150－160.

［5］Saurel R，Abgrall R.A simple method for compressible multifluid flows［J］.SIAM J Sci Comp，1999，21（3）：1115－1145.

［6］Shyue K-M.An efficient shock-capturing algorithm for compressible multicomponent problems［J］.J Comput Phys，1998，1：208－242.

［7］柏劲松，陈森华，李平.多介质非守恒律欧拉方程的数值计算方法［J］.爆炸与冲击，2001，21（4）：265－270.

［8］Shyue K-M.A fluid-mixture type algorithm for compressible multicomponent flow with Van der Waals equation of state［J］.J Comput Phys，1999，1：43－88.

［9］柏劲松，陈森华，钟敏.可压缩密实介质多流体高精度欧拉算法［J］.高压物理学报，2002，16（3）：204－212.

［10］Shyue K-M.A fluid-mixture type algorithm for compressible multicomponent flow with Mie-Gruneisen equation of state［J］.J Comput Phys，2001，171（2）：678－707.

［11］Ward G，Pullin D.A hybrid，center-difference，limiter method for simulations of compressible multicomponent flows with Mie-Grüneisen equation of state［J］.J Comput Phys，2010，229：2999－3018.

［12］Wu Zongduo，Zong Zhi.Numerical calculation of multi-component conservative Euler equations under Mie-Gruneisen equation of state［J］.Chinese J Comput Phys，2011，28（6）：113－119.

［13］Abgrall R，Karni S.Computations of compressible multifluids［J］.J Couput Phys，2001，169：594－623.

［14］Karni S.Multi-component flow calculations by a consistent primitive algorithm［J］.J Comput Phys，1994，112：31－43.

［15］Karni S.Hybrid multifluid algorithms［J］.SIAM J Sci Comput，1996，17：1019－1039.

［16］Allaire G，Clerc S，Kokh S.A five-equation model for the simulation of interfaces between compressible fluids［J］.J Comput Phys，2002，181：577－616.

［17］Saurel R，Abgrall R.Some models and methods for compressible multifluid and multiphase flows［J］.J Comput Phys，1999，150：425－467.

［18］Wang Chenxing，Tang Weijun，et al.Numerical computations of the refraction of a shock wave at interface by multi-component PPM algorimthm［J］.Chinese J Comput Phys，2004，21（6）：532－537.

［19］张雪莹，赵宁.基于体积分数形式的多介质流动数值模拟［J］.南京航空航天大学学报，2005，37（4）：436－440.

［20］Favrie N，Gavrilyuk S.Diffuse interface model for compressible fluid-compressible elastic-plastic solid interaction［J］.J Comput Phys，2012，231：2695－2723.

［21］Kokh S，Lagoutière F.An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model［J］.J Comput Phys，2010，229：2773－2809.

［22］Shyue K-M.An anti-diffusion based Eulerian interface-sharpening algorithm for compressible two-phase flow with cavitation［C］.Proceedings of the 8th International Symposium on Cavitation CAV2012-Abstract No.198，August 14－16，2012，Singapore.

［23］Pan L，Zhao G，Tian B，Wang S.A gas kinetic for the Baer-Nuziato two-phase flow model［J］.J Comput Phys，2012，231：7518－7536.

［24］LeVeque R.Wave propagation algorithms for multi-dimensional hyperbolic systems［J］.J Comput Phys，1997，131：327－353.

Algorithm Preserving Mass Fraction Maximum Principle for Multi-component Flows

TANG Weijun1，JIANG Lang2，CHENG Junbo1
（1.Institute of Applied Physics and Computational Mathematics，Laboratory of Computational Physics，Beijing 100088，China；2.Graduate School of CAEP，P．O．Box 2101，Beijing 100088，China）

We propose a new method for compressible multi-component flows with Mie-Gruneisen equation of state based on mass fraction.The model preserves conservation law of mass，momentum and total energy for mixture flows.It also preserves conservation of mass of all single components.Moreover，it prevents pressure and velocity from jumping across interface that separate regions of different fluid components.Wave propagation method is used to discretize this quasi-conservation system.Modification of numerical method is adopted for conservative equation of mass fraction.This preserves the maximum principle of mass fraction.The wave propagation method which is not modified for conservation equations of flow components mass，cannot preserve the mass fraction in the interval［0，1］.Numerical results confirm validity of the method.

Mie-Gruneisen equation of state；maximum principle of mass fraction；wave propagation method

date： 2013－06－17；Revised date： 2013－10－10

O241.82

A

1001-246X（2014）03-0292-15

2013－06－17；

2013－10－10

国家自然科学基金（11272064、11172050）及中国工程物理研究院重点基金（2013A0202011）资助项目

唐维军（1965－），男，湖南湘潭，研究员，博士，从事计算数学和计算流体力学研究，E-mail：tang_weijun＠iapcm.ac.cn