唐維軍， 蔣 浪， 程軍波

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

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

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

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

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

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

0 引言

可壓縮多介質(zhì)復(fù)雜流動有很廣泛的應(yīng)用，包括爆轟，燃燒，慣性約束聚變，天體物理等諸多領(lǐng)域.研究這些復(fù)雜流動的目的是要了解多介質(zhì)流動的流體動力學(xué)機(jī)制.這涉及到多介質(zhì)流動的數(shù)學(xué)建模及相關(guān)模型的數(shù)值離散方法.一般而言，如不考慮粘性、表面張力、熱力學(xué)效應(yīng)等物理因素，單純使用歐拉方程描述單介質(zhì)流體無粘流動時，其數(shù)學(xué)模型由質(zhì)量、動量、能量的三個守恒律組成的完全守恒系統(tǒng)構(gòu)成，且該完全守恒系統(tǒng)可用近三十年來發(fā)展成熟的各種守恒格式進(jìn)行數(shù)值離散.但如流動涉及多介質(zhì)時，僅用完全守恒模型來描述卻難以實現(xiàn).在一些非常特殊的情形，有些作者嘗試用完全守恒系統(tǒng)來描述多介質(zhì)流動［1－3］，但涉及混合流動時，相應(yīng)的數(shù)值格式處理則變得相對復(fù)雜；并且在計算純界面問題時，界面附近壓力或速度的數(shù)值解可能會出現(xiàn)非物理振蕩.目前比較公認(rèn)的是采用擬守恒模型來描述多介質(zhì)復(fù)雜流動.例如對理想氣體狀態(tài)方程一維兩流體組份流動，利用對純接觸問題界面壓力和速度數(shù)值解要求保持一致這種思想，Abgrall［4］最早揭示模型可由混合介質(zhì)質(zhì)量、動量和能量守恒方程，以及與混合狀態(tài)方程參數(shù)有關(guān)的量的一個對流方程來描述.并且Abgrall［4］還確認(rèn)該對流方程的離散格式不是獨立的，而與三個守恒方程的離散格式有關(guān).這個事實顯示多介質(zhì)流動的數(shù)學(xué)模型和數(shù)值格式與單介質(zhì)流動情形相比要復(fù)雜.隨后Abgrall［4］思想被推廣到剛性氣體狀態(tài)方程［5－7］、van der Waals氣體狀態(tài)方程［8］、密實介質(zhì)狀態(tài)方程［9］、以及Mie-Gruneisen狀態(tài)方程［10－12］.此外，還有眾多的其它模型來描述多介質(zhì)流動，如加入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ì)量、動量、和能量守恒，還保持各組份分質(zhì)量守恒，且對純接觸問題保持界面壓力速度一致.本文使用波傳播算法［7，10，24］離散系統(tǒng)，但不是象文［12］對所有守恒方程進(jìn)行相同數(shù)值離散，而是在對分介質(zhì)質(zhì)量守恒方程離散時，采用Larrouturou［1］對完全守恒系統(tǒng)設(shè)計的質(zhì)量分?jǐn)?shù)保極值原理算法，對波傳播算法進(jìn)行修改.即對分介質(zhì)質(zhì)量守恒方程的數(shù)值通量進(jìn)行修改，同時引進(jìn)一個類似CFL條件的時間步長約束.其目的是確保守恒系統(tǒng)離散時，各介質(zhì)質(zhì)量分?jǐn)?shù)始終在［0，1］范圍.如對擬守恒系統(tǒng)直接采用單介質(zhì)流動常用的那些數(shù)值格式，如TVD格式，MUSCL格式，PPM格式等，都無法確保質(zhì)量分?jǐn)?shù)在［0，1］范圍.Abgrall［4］在用Roe和MUSCL格式處理理想氣體狀態(tài)方程兩組份完全守恒模型數(shù)值離散時，借用所謂壓力速度界面一致條件來導(dǎo)出狀態(tài)方程混合參量的某個特殊離散格式，在這種情形質(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ì)量、動量、和能量守恒，還保持各組份分質(zhì)量守恒，保留了Shyue［5］的等效方程，整個模型是一個擬守恒系統(tǒng).不妨設(shè)整個流場有2種組份，其模型系統(tǒng)為

其中ρ，u，p，E分別表示密度（kg·m－3），速度（m·s－1），壓力（Pa），比總能量，Y表示第一個組份質(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個守恒方程的子矩陣記為AC，而對應(yīng)前4個守恒方程的守恒通量

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按下式計算

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

線性化矩陣AH的特征值為

其中Roe平均聲速為

這里KH＝

3 數(shù)值格式

基于波傳播算法［7，10，24］來計算擬守恒系統(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個守恒方程的數(shù)值通量進(jìn)行修改.Larrouturou［1］提出的質(zhì)量分?jǐn)?shù)保極值原理算法原用來處理理想氣體狀態(tài)方程兩組份完全守恒系統(tǒng)，其特點是采用數(shù)值格式時，能保持質(zhì)量分?jǐn)?shù)保極值原理，即完全守恒系統(tǒng)的分介質(zhì)質(zhì)量分?jǐn)?shù)始終在［0，1］范圍.本文將Larrouturou的這個算法推廣到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ǔ)上，增加了采用波限制器的二階校正項，即

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

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

這里

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

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

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

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

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

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

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

3.3 時間步長約束

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

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

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

由此得整個系統(tǒng)（9）采用波傳播算法和質(zhì)量分?jǐn)?shù)保極值原理的時間步長為從后面數(shù)值算例實際情況來看，大部分情況下，ΔtMaxP大于ΔtCFL，因此上述時間步長限制基本沒有增加計算量.當(dāng)然，這個實際經(jīng)驗取決于式（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分別為采用波傳播算法一階和二級格式，且同時采用質(zhì)量分?jǐn)?shù)保極值原理算法進(jìn)行分介質(zhì)質(zhì)量通量修正.將采用算法B2，網(wǎng)格數(shù)為2 000的計算結(jié)果作為參考解.網(wǎng)格剖分?jǐn)?shù)為100，兩種方法（包括一階和二級格式）的計算結(jié)果進(jìn)行了對比，結(jié)果顯示，兩種方法計算得到的物理量差距很小，但質(zhì)量分?jǐn)?shù)是否保極值原理卻差距很大.使用算法B1和B2，網(wǎng)格數(shù)分別為100，200，400的加密計算顯示了本算法計算的收斂性.在下面5個實際算例計算中，第四個算例的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時刻的密度、速度、壓力、聲速以及質(zhì)量分?jǐn)?shù).解的結(jié)構(gòu)包括左向傳播的稀疏波，右向傳播的接觸間斷和激波.算法A1和B1，A2和B2的計算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近參考解.圖2是兩種方法的一階和二階格式計算的各時刻質(zhì)量分?jǐn)?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的加密計算結(jié)果，顯示計算結(jié)果是收斂的.

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

圖4是t＝50 μs時刻的密度、速度、壓力、聲速以及質(zhì)量分?jǐn)?shù).解的結(jié)構(gòu)包括左向和右向傳播的兩個激波，以及左向傳播的接觸間斷.算法A1和B1，A2和B2的計算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖5是兩種方法的一階和二階格式計算的各時刻質(zhì)量分?jǐn)?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的加密計算結(jié)果，顯示計算結(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時刻的密度、速度、壓力、熱內(nèi)能以及質(zhì)量分?jǐn)?shù).圖中包含左向和右向傳播的兩個激波，以及右向傳播的接觸間斷.算法A1和B1，A2和B2的計算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖8是兩種方法的一階和二階格式計算的各時刻質(zhì)量分?jǐn)?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的加密計算結(jié)果，顯示計算結(jié)果是收斂的.

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

初始條件為

圖10是t＝73 μs時刻的密度、速度、壓力、熱內(nèi)能以及質(zhì)量分?jǐn)?shù).解的結(jié)構(gòu)包含左向傳播的稀疏波和右向傳播的激波，以及右向傳播的接觸間斷.算法A1和B1，A2和B2的計算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖11是兩種方法的一階和二階格式計算的各時刻質(zhì)量分?jǐn)?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的加密計算結(jié)果，顯示計算結(jié)果是收斂的.

初始條件為

圖13是t＝120 μs時刻的密度、速度、壓力、Gruneisen系數(shù)Γ以及質(zhì)量分?jǐn)?shù).圖中包含左向傳播的系數(shù)波和右向傳播的激波和接觸間斷.算法A1和B1，A2和B2的計算結(jié)果顯示差距甚為細(xì)微，但二階格式比一階格式要更靠近細(xì)網(wǎng)格結(jié)果.圖14是兩種方法的一階和二階格式計算的各時刻質(zhì)量分?jǐn)?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的加密計算結(jié)果，顯示計算結(jié)果是收斂的.

5 結(jié)論

對Mie-Gruneisen狀態(tài)方程多流體組份流動問題，有別于Shyue［10］采用的基于體積分?jǐn)?shù)擬守恒方程組模型，本文采用基于質(zhì)量分?jǐn)?shù)的Mie-Gruneisen狀態(tài)方程多流體組份模型.對該模型進(jìn)行數(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］.爆炸與沖擊，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］柏勁松，陳森華，鐘敏.可壓縮密實介質(zhì)多流體高精度歐拉算法［J］.高壓物理學(xué)報，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ì)流動數(shù)值模擬［J］.南京航空航天大學(xué)學(xué)報，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）及中國工程物理研究院重點基金（2013A0202011）資助項目

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