真實虛擬流方法在多介質可壓縮流動模擬中的應用＊

陳 宏，朱衛兵，張小彬，孫潤鵬，郭金鑫

（1.哈爾濱工程大學航天與建筑工程學院，黑龍江 哈爾濱 150001；2.中國航天科工集團三十一研究所高超聲速沖壓發動機技術重點實驗室，北京 100074）

多介質可壓縮流應用范圍較廣泛，如水下爆炸、航空發動機及超音速沖壓發動機燃料的混合和燃燒效率的提高、慣性約束熱核聚變的控制等［1－2］，而多介質可壓縮流中激波及大密度比介質界面附近的計算精度問題一直是計算流體力學研究中的難點和熱點。由于在界面兩側流體的狀態方程不同，流體密度在界面處出現階躍，采用傳統的單相流方式處理此類問題時，界面處會出現非物理振蕩。當界面變化劇烈時，會使計算無法持續進行。因此，界面附近流場的處理就顯得尤為重要。

為了解決多介質流計算中產生的振蕩現象，R.P.Fedkiw等［3］提出了采用GFM處理多介質流動問題。該方法在抑制界面附近非物理振蕩方面取得了較好效果，但是在解決強間斷和大密度比問題（如氣－水界面問題）時會引起虛假的物理解。R.Caiden等［4］用GFM處理氣－水界面問題，通過對水介質中的速度進行插值和對空氣側的壓力進行插值定義了ghost流體的狀態。該方法盡管處理氣－水界面表現優異，但不適合用于求解氣－氣界面。T.G.Liu等［5］詳細分析了原始GFM的缺陷，認為其處理波系與界面作用時不夠精確，在此基礎上他們提出了修正虛擬流法（modified ghost fluid method，MGFM），利用特征關系和激波跳躍條件求解當地Riemann解問題，進而更新對應的ghost點。C.M.Wang等［6］對MGFM進行了改進，提出了新的處理方式，即RGFM。該方法對靠近界面的真實流體狀態采用Riemann解進行了修正，從而減小了計算誤差。在此基礎上，C.M.Wang等［6］初步計算了氣－水界面變化問題。陳榮三等［7］通過在界面處構造Riemann問題對原始GFM進行了改進，對一維強激波與氣－氣、氣－液界面的相互作用問題以及射流問題進行了計算。本文中，將level set方法［8］、HLLC方法［9］結合GFM及RGFM運用于氣－氣、氣－水界面兩相流的數值計算，并比較用2種ghost方法得到的計算結果。在氣－氣界面兩相流中，用RGFM捕捉到的波系和界面分布更準確，同時采用RGFM處理氣－水界面復雜變化時能得到較好的結果，而采用GFM處理氣－水界面問題時遇到了困難，程序無法進行下去。

1 計算方法

1.1 基本方程

守恒方程為二維Euler方程，表示如下：

式中：x和y分別為橫向和縱向坐標，ρ和p分別為密度和壓力，u和v為分別為x和y方向上的速度，體積總能量

由于HLLC格式在HLL（Harten－Lax－Van Leer）格式的基礎上考慮了接觸間斷和剪切波問題，因此，本文中采用HLLC格式離散數值通量

式中：Sl和Sr分別為左、右波速，Sm為中間波速，詳見文獻［9］。中間狀態變量可由Rankine－Hugonoit關系式得到：

式中：k＝l，r；ρ＊k為左右2束波之間的物質密度；q≡unx＋vny。

中間波速和最大、最小波速分別估計為：

式中：c為聲速，u＊、v＊為Roe平均變量。

1.2 level set方法

level set方程為基本思想為：構造函數φ（x，t），使得在任意時刻，運動界面Γ（t）恰是φ（x，t）的零等值面，一般可取φ（x，t）為點x到界面Γ（0）的符號距離。φ（x，t）以適當的速度移動，在任意時刻，只要求出φ（x，t）的值，物質界面的位置也就確定了。

level set方法通過距離函數φ確定物質界面的位置。x到Γ（0）的符號距離

level set方程的數值求解：方程（12）空間項采用五階 WENO格式［10］求解，時間項采用三階TVDRunge－Kutta方法［11］求解。level set重構方程表示為

1.3 真實虛擬流法（RGFM）

與原始的GFM［3］相比，真實虛擬流法（RGFM）允許修改界面處真實流體的壓力、速度和密度，使計算結果更精確。通過求解Riemann解得出界面處的參數，并將其賦予虛擬網格節點，如圖1所示。假設物質界面位于節點i和i＋1之間，首先在界面處定義Riemann問題，界面左側狀態Ul＝Ui，右側狀態Ur＝Ui＋1，然后在此基礎上求解Riemann解（本文在此處采用精確Riemann解［9］），可以得到界面處的速度ui和壓力pi以及界面左側的密度ρi，l和右側的密度ρi，r，將所求得的Riemann解重新賦值（ui、pi和ρi，l）給ghost點（i＋1和i＋2），同時將（ui、pi和ρi，l）賦給真實流體中緊鄰界面的節點i。

圖1 真實虛擬流方法Fig.1Real ghost fluid method

2 數值模擬

2.1 一維算例—Sod激波管問題

初始時刻壓力比p2／p1＝10的Sod激波管問題初始條件為：當x＜0.6時，ρ2＝1.000，p2＝10.0，u2＝0.0，γ＝1.4；當x＞0.6時，ρ1＝0.125，p1＝1.0，u1＝0.0，γ＝1.4。計算結果如圖2所示，計算時間為0.035s，圖中給出了用GFM和RGFM得到的密度、速度和壓力的數值解。從圖2中可以看出，與解析解相比，RGFM捕捉到的激波、膨脹波及界面的分布明顯好于GFM的計算結果。

壓力比p2／p1＝1 000時，由于初始時刻流場中存在較大的間斷（壓力和密度），隨著時間的延長，在流場中產生了向右運動的強激波和接觸間斷，同時還產生了向左運動的膨脹波。由于初期接觸間斷的速度較高，其緊隨激波向右運動。圖3為t＝0.005s時計算結果與解析解的比較。由圖3可以看出，在這種極端的情況下采用RGFM得到的計算結果與解析解仍然符合較好，而采用GFM得到的結果則在界面處出現了振蕩現象，得到了非物理解。

圖2 p2／p1＝10時的Sod激波管問題Fig.2Sod’s shock tube problem at p2／p1＝10

圖3 p2／p1＝1 000時的Sod激波管問題Fig.3Sod’s shock tube problem at p2／p1＝1 000

2.2 二維算例

2.2.1 激波打柱形氦氣泡

激波與氣泡之間的相互作用應用背景廣泛，從天體物理到人體組織的空泡損傷（如對腎結石的震波碎石等）。J.L.Hass等［13］實驗研究了激波與圓柱形氦氣泡的相互作用。本算例即以該實驗為基礎，研究了激波由重氣體進入輕氣體的物理過程，本算例也在該實驗的基礎上對馬赫數為1.22的激波在空氣中與氦氣泡的相互作用過程進行了數值模擬。由于計算域對稱，只計算了上半區域。圖4為計算域示意圖，對幾何尺寸進行了量綱一化，以下各算例與此相同。氣泡內區域為Ⅰ，氣泡外激波前區為Ⅱ，激波后區域為Ⅲ，入射激波速度為1.22Ma，從右向左運動。計算網格為1 600×200，其量綱一初始條件［3］為：氣泡內，ρ＝0.138 0，p＝1.0000，u＝0.000，v＝0.0，γ＝1.67；激波前，ρ＝1.0000，p＝1.000，u＝0.000，v＝0.0，γ＝1.4；當x＞225時，ρ＝1.376 4，p＝1.569 8，u＝－0.394，v＝0.0，γ＝1.4。

圖5給出了對該物理過程采用RGFM得到的計算結果和實驗圖像［13］，為了能清晰地觀察到波系的變化，采用紋影技術［14］進行了后處理。入射激波在空氣中從右向左運動，首先打到氣泡右側界面上，右界面受此作用開始向左運動，同時在右界面上形成Richtmyer－Meshkov（RM）不穩定性。入射激波打到界面上產生的反射波為膨脹波，而折射波為激波；由于激波在氦氣泡中的傳播速度比在空氣中的快，所以可以看到折射激波在氣泡中形成弧形激波走在了入射激波的前面，這樣就在氣泡上下界面形成了速度差，使得氣泡上下界面開始出現Kelvin－Helmholtz（KH）不穩定性，如圖5（a）所示。

當氦氣泡中的弧形折射激波遇到氣泡左界面時，同樣會發生反射和折射，反射波和折射波均為激波，并且在左界面也產生了RM不穩定性。由圖5（b）可以看出，此時反射波較弱，在實驗紋影圖中已經很難觀察到；而折射波以圓弧形從氦氣泡中進入空氣，在圖中可以看到圓弧形激波遇到激波管上下壁面發生反射，且該反射波再次與氣泡的上下界面作用，形成新的RM和KH不穩定性；此時氣泡在波系的作用下呈半圓形。

由圖5（c）和（d）可以看出，氣泡右側界面受氣流和不穩定性影響向內塌陷，并逐漸形成丁狀結構。隨著丁狀結構的發展，最終導致氣泡斷裂，并在斷裂處形成射流。

圖4 激波管內激波打氣泡示意圖Fig.4Sketch of the dimensionless computational domain for shock－bubble interaction in a shock tube

圖5 采用RGFM得到的計算結果與實驗圖像［13］的比較Fig.5Comparison between computational results by RFGM and experimental images［13］

圖6為激波打氣泡中后期氣泡斷裂后的形態圖像，圖中同時給出了對應時刻采用GFM得到的計算結果。從圖6中可以看出，采用RGFM得到的計算結果與實驗結果符合更好，而采用GFM得到的計算結果誤差稍大。目前關于激波與氦氣泡相互作用問題的數值研究中，大部分文獻只給出了氣泡變化的初期形態，很少給出氣泡大變形甚至斷裂后的形態。本文算法給出了整個氣泡變化過程的圖像，而且均與實驗結果［13］符合較好，驗證了本文算法的適用性。

圖6 用不同方法得到的激波打氦氣泡中后期氣泡斷裂后的形態Fig.6Shapes of helium bubble after fracture in the midanaphase of the interaction between shock wave and helium bubble by different methods

2.2.2 水下激波打氣泡

本算例主要分析水下一入射激波掃過圓柱形氣泡后的物理變化過程，計算域如圖7所示，各邊界均為自由邊界。初始條件為：氣泡內，ρ＝1.00，p＝1.0，u＝0.0000，v＝1.0，γ＝1.4；激波后，ρ＝1 220.75，p＝9 120.0，u＝1.284 1，v＝0.0，γ＝7.0；其余為，ρ＝1 000.0，p＝1.00，u＝0.0000，v＝0.0，γ＝7.0。

采用GFM處理氣－水界面時，程序不能順利計算下去，因此本部分只給出了采用RGFM得到的計算結果。圖8給出了網格數為840×720的情況下不同時刻氣泡界面的位置，從圖中可以看出氣泡從變形到破碎的變化過程。

圖7 計算域Fig.7Computational domain

圖8 不同時刻的界面分布Fig.8Distribution of the interface at different times

圖9給出了不同時刻流場密度變化的紋影圖。激波在水中從左向右運動掃過柱形氣泡物質界面，在空氣側產生一道折射激波，而在水中產生反射膨脹波，如圖9（a）所示，由于界面兩側聲阻不同，且差異較大，氣泡內折射激波速度滯后于入射激波，氣泡受到波系作用而產生變形，呈現出“腎”形。隨著時間的延長，氣泡內折射激波在氣泡右側界再次發生反射和折射，同時在氣泡凹陷處形成水射流，如圖9（b）所示。圖9（c）～（d）中水射流已經將氣泡截斷，并產生了沖擊波。受該沖擊波的影響，斷裂的小氣泡中形成了二次水射流現象。該算例驗證了本文算法處理多介質流中激波與大密度比界面相互作用問題的適用性。

圖9 不同時刻的密度紋影Fig.9Density schlieren images at different times

3 結 論

（1）將RGFM引入到多介質可壓縮流的數值模擬中，研究激波與大密度比介質界面的相互作用問題，實現了界面不穩定的精確模擬。研究和開發了適用于可壓縮多介質流體界面（氣－氣界面、氣－水界面）不穩定性發展演化的計算方法和程序。

（2）一維和二維算例表明，RGFM在處理氣－氣界面時優于原始GFM，能得到更精確的計算結果。

（3）在水下激波打氣泡的算例中，采用RGFM成功的捕捉到了大密度比（1 000∶1）界面從變形到破碎的變化過程，證明該方法可以處理大密度比多介質流動中的界面大變形問題。

［1］Brouillette M.The Richtmyer－Meshkov instability［J］.Annual Review of Fluid Mechanics，2002，34：445－468.

［2］Waitz I A，Marble F E，Zukoski E E.Investigation of a contoured wall injector for hyper－velocity mixing augmentation［J］.AIAA Journal，1993，31（6）：1014－1021.

［3］Fedkiw R P，Aslam T，Merriman B，et al.A non－oscillatory Eulerian approach to interfaces in multimaterial flows：The ghost fluid method［J］.Journal of Computational Physics，1999，152（2）：457－492.

［4］Caiden R，Fedkiw R P，Anderson C.A numerical method for two－phase flow consisting of separate compressible and incompressible regions［J］.Journal of Computational Physics，2001，166（1）：1－27.

［5］Liu T G，Khoo B C，Yeo K S.Ghost fluid method for strong shock impacting on material interface［J］.Journal of Computational Physics，2003，190（2）：651－681.

［6］Wang C W，Liu T G，Khoo B C.A real ghost fluid method for the simulation of multimedium compressible flow［J］.SIAM Journal on Scientific Computing，2006，28（1）：278－302.

［7］陳榮三，盛萬成.用改進的ghost fluid方法模擬強激波與界面的相互作用［J］.爆炸與沖擊，2008，28（2）：144－148.

Chen Rong－san，Sheng Wan－cheng.A modified ghost fluid method for strong shock wave impacting on material interface［J］.Explosion and Shock Waves，2008，28（2）：144－148.

［8］Osher S，Sethian J A.Fronts propagating with curvature－dependent speed：Algorithms based on Hamilton－Jacobi formulations［J］.Journal of Computational Physics，1988，79（1）：12－49.

［9］Toro E F.Riemann solvers and numerical methods for fluid dynamics：A pratical introduction［M］.2nd ed.Springer，1999：115－162.

［10］Jiang G S，Shu C W.Efficient implementation of weighted ENO schemes［J］.Journal of Computational Physics，1996，126（1）：202－228.

［11］Chen S，Merriman B，Osher S，et al.A simple level set method for solving Stefan problem［J］.Journal of Computational Physics，1997，135（1）：8－29.

［12］Tang H S，Huang D.A second－order accurate capturing scheme for 1Dinviscid flows of gas and water with vacuum zones［J］.Journal of Computational Physics，1996，128（2）：301－318.

［13］Haas J F，Sturtevant B.Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities［J］.Journal of Fluid Mechanics，1987，181（1）：41－76.

［14］Quirk J J，Karni S.On the dynamics of a shock－bubble interaction［J］.Journal of Fluid Mechanics，1996，318：129－163.