多組份流動質量分數保極值原理算法

唐維軍， 蔣 浪， 程軍波

（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