多組份流動(dòng)質(zhì)量分?jǐn)?shù)保極值原理算法

唐維軍， 蔣 浪， 程軍波

（1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，計(jì)算物理實(shí)驗(yàn)室，北京 100088；2.中國工程物理研究院研究生部，北京 100088）

多組份流動(dòng)質(zhì)量分?jǐn)?shù)保極值原理算法

唐維軍1， 蔣 浪2， 程軍波1

（1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，計(jì)算物理實(shí)驗(yàn)室，北京 100088；2.中國工程物理研究院研究生部，北京 100088）

對基于質(zhì)量分?jǐn)?shù)的Mie-Gruneisen狀態(tài)方程多流體組份模型提出了新的數(shù)值方法.該模型保持混合流體的質(zhì)量、動(dòng)量、和能量守恒，保持各組份分質(zhì)量守恒，在多流體組份界面處保持壓力和速度一致.該模型是擬守恒型方程系統(tǒng).對該模型系統(tǒng)的離散采用波傳播算法.與直接對模型中所有守恒方程采用相同算法不同的是，在處理分介質(zhì)質(zhì)量守恒方程時(shí)，對波傳播算法進(jìn)行了修正，使之滿足質(zhì)量分?jǐn)?shù)保極值原理.而不作修改的算法則不能保證質(zhì)量分?jǐn)?shù)在［0，1］范圍.數(shù)值實(shí)驗(yàn)驗(yàn)證了該方法有效.

Mie-Gruneisen狀態(tài)方程；質(zhì)量分?jǐn)?shù)保極值原理；波傳播算法

0 引言

可壓縮多介質(zhì)復(fù)雜流動(dòng)有很廣泛的應(yīng)用，包括爆轟，燃燒，慣性約束聚變，天體物理等諸多領(lǐng)域.研究這些復(fù)雜流動(dòng)的目的是要了解多介質(zhì)流動(dòng)的流體動(dòng)力學(xué)機(jī)制.這涉及到多介質(zhì)流動(dòng)的數(shù)學(xué)建模及相關(guān)模型的數(shù)值離散方法.一般而言，如不考慮粘性、表面張力、熱力學(xué)效應(yīng)等物理因素，單純使用歐拉方程描述單介質(zhì)流體無粘流動(dòng)時(shí)，其數(shù)學(xué)模型由質(zhì)量、動(dòng)量、能量的三個(gè)守恒律組成的完全守恒系統(tǒng)構(gòu)成，且該完全守恒系統(tǒng)可用近三十年來發(fā)展成熟的各種守恒格式進(jìn)行數(shù)值離散.但如流動(dòng)涉及多介質(zhì)時(shí)，僅用完全守恒模型來描述卻難以實(shí)現(xiàn).在一些非常特殊的情形，有些作者嘗試用完全守恒系統(tǒng)來描述多介質(zhì)流動(dòng)［1－3］，但涉及混合流動(dòng)時(shí)，相應(yīng)的數(shù)值格式處理則變得相對復(fù)雜；并且在計(jì)算純界面問題時(shí)，界面附近壓力或速度的數(shù)值解可能會(huì)出現(xiàn)非物理振蕩.目前比較公認(rèn)的是采用擬守恒模型來描述多介質(zhì)復(fù)雜流動(dòng).例如對理想氣體狀態(tài)方程一維兩流體組份流動(dòng)，利用對純接觸問題界面壓力和速度數(shù)值解要求保持一致這種思想，Abgrall［4］最早揭示模型可由混合介質(zhì)質(zhì)量、動(dòng)量和能量守恒方程，以及與混合狀態(tài)方程參數(shù)有關(guān)的量的一個(gè)對流方程來描述.并且Abgrall［4］還確認(rèn)該對流方程的離散格式不是獨(dú)立的，而與三個(gè)守恒方程的離散格式有關(guān).這個(gè)事實(shí)顯示多介質(zhì)流動(dòng)的數(shù)學(xué)模型和數(shù)值格式與單介質(zhì)流動(dòng)情形相比要復(fù)雜.隨后Abgrall［4］思想被推廣到剛性氣體狀態(tài)方程［5－7］、van der Waals氣體狀態(tài)方程［8］、密實(shí)介質(zhì)狀態(tài)方程［9］、以及Mie-Gruneisen狀態(tài)方程［10－12］.此外，還有眾多的其它模型來描述多介質(zhì)流動(dòng)，如加入level set方程的擴(kuò)展歐拉方程組模型［13］，Karni的原始變量模型［14］和雜交模型［15］，Allaire等的5方程模型［16］，Saurel和Abgrall的兩相流7方程模型［17］，使用體積分?jǐn)?shù)的4方程模型［18－19］，反擴(kuò)散界面模型等［20－22］，兩相流BGK模型［23］.至于有關(guān)上述模型的數(shù)值格式，散見上述各文及相關(guān)引文，此處不一一列舉.

討論Mie-Gruneisen狀態(tài)方程一維模型的數(shù)值離散格式.該模型基于Abgrall［4］等效方程思想，且分介質(zhì)質(zhì)量保持守恒的模型［12］的數(shù)值離散格式.與基于體積分?jǐn)?shù)對流方程模型［10］相比，本文所用模型［23］不僅能保持混合流體的質(zhì)量、動(dòng)量、和能量守恒，還保持各組份分質(zhì)量守恒，且對純接觸問題保持界面壓力速度一致.本文使用波傳播算法［7，10，24］離散系統(tǒng)，但不是象文［12］對所有守恒方程進(jìn)行相同數(shù)值離散，而是在對分介質(zhì)質(zhì)量守恒方程離散時(shí)，采用Larrouturou［1］對完全守恒系統(tǒng)設(shè)計(jì)的質(zhì)量分?jǐn)?shù)保極值原理算法，對波傳播算法進(jìn)行修改.即對分介質(zhì)質(zhì)量守恒方程的數(shù)值通量進(jìn)行修改，同時(shí)引進(jìn)一個(gè)類似CFL條件的時(shí)間步長約束.其目的是確保守恒系統(tǒng)離散時(shí)，各介質(zhì)質(zhì)量分?jǐn)?shù)始終在［0，1］范圍.如對擬守恒系統(tǒng)直接采用單介質(zhì)流動(dòng)常用的那些數(shù)值格式，如TVD格式，MUSCL格式，PPM格式等，都無法確保質(zhì)量分?jǐn)?shù)在［0，1］范圍.Abgrall［4］在用Roe和MUSCL格式處理理想氣體狀態(tài)方程兩組份完全守恒模型數(shù)值離散時(shí)，借用所謂壓力速度界面一致條件來導(dǎo)出狀態(tài)方程混合參量的某個(gè)特殊離散格式，在這種情形質(zhì)量分?jǐn)?shù)滿足保極值原理.但這種方法難以推廣到本文討論的Mie-Gruneisen狀態(tài)方程情形.特別值得提到的是，文［2］也采用了Larrouturou［1］通量修正格式，但文［2］模型與本文模型不同，且僅限理想氣體狀態(tài)方程，無法推廣到其它狀態(tài)方程.另外，本文對波傳播算法二階格式的通量修改格式，與Larrouturou［1］二階MUSCL的通量修改格式不一樣，后者對質(zhì)量分?jǐn)?shù)的重構(gòu)需加額外的約束.本文算例表明，這種數(shù)值格式的有效性.

1 模型方程

文［12］的模型不僅能保持混合流體的質(zhì)量、動(dòng)量、和能量守恒，還保持各組份分質(zhì)量守恒，保留了Shyue［5］的等效方程，整個(gè)模型是一個(gè)擬守恒系統(tǒng).不妨設(shè)整個(gè)流場有2種組份，其模型系統(tǒng)為

其中ρ，u，p，E分別表示密度（kg·m－3），速度（m·s－1），壓力（Pa），比總能量，Y表示第一個(gè)組份質(zhì)量分?jǐn)?shù).Mie-Gruneisen狀態(tài)方程為

和內(nèi)能.將方程組（1）－（7）寫成擬線性形式

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

矩陣A對應(yīng)的特征值分別為

特征值向量Λ對應(yīng)的右特征向量記為R＝（r1…，r7）．這些特征值和右特征向量滿足Ark＝λkrk，k＝1，…，7．

Jacobia矩陣A對應(yīng)前4個(gè)守恒方程的子矩陣記為AC，而對應(yīng)前4個(gè)守恒方程的守恒通量

2 一維Riemann逼近問題

考慮一維Riemann逼近問題

其中，qL＝qi，qR＝qi＋1．設(shè)矩陣AH是Jacobia矩陣A在Roe平均狀態(tài)qH的局部線性化，即AH＝A（qH）．Roe線性化要滿足如下兩式

上述過程中，平均量（1／Γ）H，（p／Γ）H須同式（14）的選取方式，且滿足如下逼近（Shyue［14］）

則平均量pH，ΓH按下式計(jì)算

要完成AH，還須補(bǔ)充φH，χH，ψH．但注意式（12）方程個(gè)數(shù)與而式（11）方程個(gè)數(shù)并不相等，沒有其它條件可用，故直接仿式（14）Roe平均給出，

線性化矩陣AH的特征值為

其中Roe平均聲速為

這里KH＝

3 數(shù)值格式

基于波傳播算法［7，10，24］來計(jì)算擬守恒系統(tǒng)（9），波傳播算法的標(biāo)準(zhǔn)有限體積方法是Godunov型迎風(fēng)格式.但本文在處理分介質(zhì)質(zhì)量守恒方程，不采用標(biāo)準(zhǔn)波傳播算法，而采用Larrouturou［1］提出的質(zhì)量分?jǐn)?shù)保極值原理算法，對擬守恒系統(tǒng)（9）的第4個(gè)守恒方程的數(shù)值通量進(jìn)行修改.Larrouturou［1］提出的質(zhì)量分?jǐn)?shù)保極值原理算法原用來處理理想氣體狀態(tài)方程兩組份完全守恒系統(tǒng)，其特點(diǎn)是采用數(shù)值格式時(shí)，能保持質(zhì)量分?jǐn)?shù)保極值原理，即完全守恒系統(tǒng)的分介質(zhì)質(zhì)量分?jǐn)?shù)始終在［0，1］范圍.本文將Larrouturou的這個(gè)算法推廣到Mie-Gruneisen方程的擬守恒系統(tǒng)（9），不僅確保純接觸問題界面壓力速度保持一致，這是Larrouturou［1］算法不曾考慮的，而且能確保擬守恒系統(tǒng)（9）的分介質(zhì)質(zhì)量分?jǐn)?shù)滿足保極值原理，這是原來波傳播算法［7，10，24］處理分介質(zhì)守恒方程所不具備的.

3.1 波傳播算法

波傳播算法［7，10］一階格式為

其中

波傳播算法［7，10］二階格式在一階格式基礎(chǔ)上，增加了采用波限制器的二階校正項(xiàng)，即

3.2 分介質(zhì)質(zhì)量守恒方程通量修正格式

上面已提到，計(jì)算擬守恒系統(tǒng)（9）的分介質(zhì)質(zhì)量守恒方程，即系統(tǒng)第4個(gè)方程

這里

考慮到質(zhì)量分?jǐn)?shù)Y＝q（4）／q（1），如直接采用波傳播算法（20）或（21）計(jì)算守恒變量q（1），q（4），則無法確保Y∈［0，1］．故方程（23）的離散采用Larrouturou［1］提出的質(zhì)量分?jǐn)?shù)保極值原理算法進(jìn)行通量修改.

雖然模型（9）整體是擬守恒方程組，波傳播算法離散（9）的守恒方程部分可寫為守恒通量形式.第一個(gè)方程，即混合質(zhì)量守恒方程離散為

其中當(dāng)波傳播算法為一階格式（20）時(shí)，通量分量（F（1））i＋1／2為

當(dāng)波傳播算法為二階格式（21）時(shí)，通量分量（F（1））i＋1／2為

利用Larrouturou質(zhì)量分?jǐn)?shù)保極值原理算法［1］，分介質(zhì)質(zhì)量守恒方程（23）的離散格式為

特別需要提到的是，Larrouturou［1］采用二階MUSCL格式進(jìn)行通量修正時(shí)，其通量分量修正與上式并不一致，即

且必須對質(zhì)量分?jǐn)?shù)重構(gòu)加上如下限制

3.3 時(shí)間步長約束

要滿足質(zhì)量分?jǐn)?shù)保極值原理，波傳播算法的上述數(shù)值格式修改，還需考慮一個(gè)類似CFL條件的時(shí)間步長約束［1］，即

當(dāng)波傳播算法為一階格式（20），流通量分量F（1）為一階精度時(shí)，

當(dāng)波傳播算法為二階格式（21），流通量分量F（1）為二階精度時(shí)，

由此得整個(gè)系統(tǒng)（9）采用波傳播算法和質(zhì)量分?jǐn)?shù)保極值原理的時(shí)間步長為從后面數(shù)值算例實(shí)際情況來看，大部分情況下，ΔtMaxP大于ΔtCFL，因此上述時(shí)間步長限制基本沒有增加計(jì)算量.當(dāng)然，這個(gè)實(shí)際經(jīng)驗(yàn)取決于式（33）中的CFL系數(shù)CCFL取值大小.

4 數(shù)值算例

本文所有算例來自文獻(xiàn)［14］.對擬守恒系統(tǒng)（9）使用兩種方法進(jìn)行數(shù)值離散.第一種方法是直接使用波傳播算法；第二種方法是在波傳播算法基礎(chǔ)上，再使用質(zhì)量分?jǐn)?shù)保極值原理對分介質(zhì)守恒方程數(shù)值通量進(jìn)行修正.由于波傳播算法還分為一階和二階格式，后面為簡明起見，記算法A1和A2分別為直接采用波傳播算法的一階和二階格式，而記算法B1和B2分別為采用波傳播算法一階和二級格式，且同時(shí)采用質(zhì)量分?jǐn)?shù)保極值原理算法進(jìn)行分介質(zhì)質(zhì)量通量修正.將采用算法B2，網(wǎng)格數(shù)為2 000的計(jì)算結(jié)果作為參考解.網(wǎng)格剖分?jǐn)?shù)為100，兩種方法（包括一階和二級格式）的計(jì)算結(jié)果進(jìn)行了對比，結(jié)果顯示，兩種方法計(jì)算得到的物理量差距很小，但質(zhì)量分?jǐn)?shù)是否保極值原理卻差距很大.使用算法B1和B2，網(wǎng)格數(shù)分別為100，200，400的加密計(jì)算顯示了本算法計(jì)算的收斂性.在下面5個(gè)實(shí)際算例計(jì)算中，第四個(gè)算例的CFL系數(shù)CCFL取值為0.75，其它算例的CFL系數(shù)CCFL均取值為0.9.

例1 （單組份問題）氣體TNT爆炸產(chǎn)物Riemann問題，材料為TNT爆炸產(chǎn)物，參考狀態(tài)方程使用如下JWL狀態(tài)方程［10］（34）－（36），材料參數(shù)Γ0，V0，e0，A，B，R1，R2具體見文［14］.

圖1是t＝12 μs時(shí)刻的密度、速度、壓力、聲速以及質(zhì)量分?jǐn)?shù).解的結(jié)構(gòu)包括左向傳播的稀疏波，右向傳播的接觸間斷和激波.算法A1和B1，A2和B2的計(jì)算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近參考解.圖2是兩種方法的一階和二階格式計(jì)算的各時(shí)刻質(zhì)量分?jǐn)?shù)極大值隨時(shí)間變化的歷史曲線.從圖中可以看出，直接采用波傳播算法的二階格式A2，其結(jié)果不能保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.而算法B1、B2及A1則保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.圖3是采用B1和B2對網(wǎng)格數(shù)分別為100，200，400的加密計(jì)算結(jié)果，顯示計(jì)算結(jié)果是收斂的.

例2 （單組份問題）鋁塊撞擊鋁塊，鋁板從右向左撞擊半無限長鋁板.狀態(tài)方程采用如下激波狀態(tài)方程［10］（37）－（39），材料參數(shù)ρ0，c0，s，Γ0，α，e0，p0見文［10］.

圖4是t＝50 μs時(shí)刻的密度、速度、壓力、聲速以及質(zhì)量分?jǐn)?shù).解的結(jié)構(gòu)包括左向和右向傳播的兩個(gè)激波，以及左向傳播的接觸間斷.算法A1和B1，A2和B2的計(jì)算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖5是兩種方法的一階和二階格式計(jì)算的各時(shí)刻質(zhì)量分?jǐn)?shù)極大值隨時(shí)間變化的歷史曲線.從圖中可以看出，直接采用波傳播算法的二階格式A2，其結(jié)果不能保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.而算法B1、B2及A1則保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.圖6是采用B1和B2對網(wǎng)格數(shù)分別為100，200，400的加密計(jì)算結(jié)果，顯示計(jì)算結(jié)果是收斂的.

例3 兩組份碰撞算例，初始為標(biāo)準(zhǔn)大氣壓條件.即壓力為p0＝1 atm，T0＝300 K.銅板以速度u＝1 500 m·s－1右行，與固體炸藥碰撞.銅和固體炸藥的參考狀態(tài)方程都使用如下CC狀態(tài)方程［10］（40）－（43），材料參數(shù)ρ0，cV，T0，A，B，ε1，ε2，Γ0見文［10］，密度ρ＝ρ0．

初始數(shù)據(jù)為

圖7是t＝85 μs時(shí)刻的密度、速度、壓力、熱內(nèi)能以及質(zhì)量分?jǐn)?shù).圖中包含左向和右向傳播的兩個(gè)激波，以及右向傳播的接觸間斷.算法A1和B1，A2和B2的計(jì)算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖8是兩種方法的一階和二階格式計(jì)算的各時(shí)刻質(zhì)量分?jǐn)?shù)極大值隨時(shí)間變化的歷史曲線.從圖中可以看出，直接采用波傳播算法的二階格式A2，其結(jié)果不能保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.而算法B1、B2及A1則保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.圖9是采用B1和B2對網(wǎng)格數(shù)分別為100，200，400的加密計(jì)算結(jié)果，顯示計(jì)算結(jié)果是收斂的.

例4 兩組份問題，材料為銅和氣體爆轟產(chǎn)物，參考狀態(tài)方程銅使用CC狀態(tài)方程（40）－（43），而氣體爆轟產(chǎn)物使用JWL狀態(tài)方程（34）－（36），材料參數(shù)見文［14］.

初始條件為

圖10是t＝73 μs時(shí)刻的密度、速度、壓力、熱內(nèi)能以及質(zhì)量分?jǐn)?shù).解的結(jié)構(gòu)包含左向傳播的稀疏波和右向傳播的激波，以及右向傳播的接觸間斷.算法A1和B1，A2和B2的計(jì)算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖11是兩種方法的一階和二階格式計(jì)算的各時(shí)刻質(zhì)量分?jǐn)?shù)極大值隨時(shí)間變化的歷史曲線.從圖中可以看出，直接采用波傳播算法的一階格式A1和二階格式A2，其結(jié)果均不能保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.而算法B1、B2則保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.圖12是采用B1和B2對網(wǎng)格數(shù)分別為100，200，400的加密計(jì)算結(jié)果，顯示計(jì)算結(jié)果是收斂的.

初始條件為

圖13是t＝120 μs時(shí)刻的密度、速度、壓力、Gruneisen系數(shù)Γ以及質(zhì)量分?jǐn)?shù).圖中包含左向傳播的系數(shù)波和右向傳播的激波和接觸間斷.算法A1和B1，A2和B2的計(jì)算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖14是兩種方法的一階和二階格式計(jì)算的各時(shí)刻質(zhì)量分?jǐn)?shù)極大值隨時(shí)間變化的歷史曲線.從圖中可以看出，直接采用波傳播算法的一階格式A1和二階格式A2，其結(jié)果不能保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.而算法B1、B2則保持質(zhì)量分?jǐn)?shù)在［0，1］范圍.圖15是采用B1和B2對網(wǎng)格數(shù)分別為100，200，400的加密計(jì)算結(jié)果，顯示計(jì)算結(jié)果是收斂的.

5 結(jié)論

對Mie-Gruneisen狀態(tài)方程多流體組份流動(dòng)問題，有別于Shyue［10］采用的基于體積分?jǐn)?shù)擬守恒方程組模型，本文采用基于質(zhì)量分?jǐn)?shù)的Mie-Gruneisen狀態(tài)方程多流體組份模型.對該模型進(jìn)行數(shù)值離散時(shí)，需要設(shè)法保持質(zhì)量分?jǐn)?shù)在［0，1］區(qū)間，否則，所得分介質(zhì)質(zhì)量有可能為負(fù).將Larrouturou［1］研究理想氣體狀態(tài)方程完全守恒系統(tǒng)的多組份模型數(shù)值離散保持質(zhì)量分?jǐn)?shù)保極值原理的格式應(yīng)用到Mie-Gruneisen狀態(tài)方程多流體組份擬守恒系統(tǒng).數(shù)值結(jié)果顯示了修改的波傳播算法格式能保持質(zhì)量分?jǐn)?shù)保極值原理.

［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］柏勁松，陳森華，李平.多介質(zhì)非守恒律歐拉方程的數(shù)值計(jì)算方法［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］柏勁松，陳森華，鐘敏.可壓縮密實(shí)介質(zhì)多流體高精度歐拉算法［J］.高壓物理學(xué)報(bào)，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ǐn)?shù)形式的多介質(zhì)流動(dòng)數(shù)值模擬［J］.南京航空航天大學(xué)學(xué)報(bào)，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

國家自然科學(xué)基金（11272064、11172050）及中國工程物理研究院重點(diǎn)基金（2013A0202011）資助項(xiàng)目

唐維軍（1965－），男，湖南湘潭，研究員，博士，從事計(jì)算數(shù)學(xué)和計(jì)算流體力學(xué)研究，E-mail：tang_weijun＠iapcm.ac.cn