水/氣多介質問題的界面處理方法*

徐 爽,趙 寧,王春武,王東紅

(1.南京 航空航天大 學航空宇航 學院,江蘇 南 京 210016; 2.南京航空航天大學理學院,江蘇 南京210016)

水/氣多介質問題的界面處理方法*

徐 爽1,趙 寧1,王春武2,王東紅2

(1.南京 航空航天大 學航空宇航 學院,江蘇 南 京 210016; 2.南京航空航天大學理學院,江蘇 南京210016)

針對不可壓縮可壓縮水/氣多介質問題,提出一種新的界面處理方法。在可壓縮水/氣界面處構造Riemann問題,在水中設音速趨于無窮大,求解 Riemann問題得到不可壓縮可壓縮水/氣界面處流體的準確流動狀態;然后以此狀態結合 GFM(ghost fluid method)方法分別為2種流體定義界面邊界條件,將兩相流 問題轉化為單相流問題計算,通過求解level set方程來跟蹤界面的位置。對各種不同的界面邊界條件定義方法進行了比較,數值模擬結果表明算法能準確地捕捉各類間斷的位置,證明了算法的有效性和穩健性。

爆炸力學;Riemann問題;GFM 方法;level set;水/氣界面 ;多介質流 動

在多介質問題的數值模擬中,由于水/氣物質屬性的巨大差異,使得水/氣多介質問題成為多介質問題的難點之一。針對具有大密度比的水/氣多 介質 問 題,傳 統 方 法 主 要 有 3 類[1]:(1)水/氣 皆 被 作 為 可壓縮流體處理,其中水的狀態方程為 Tait方程;(2)水/氣皆作為不可壓縮流體,這類方法往往用來處理低速水/氣多介質運動;(3)水作為不可壓流體,氣作為可壓縮流體,相應的分別采用不可壓方程和可壓縮歐拉方程作為水/氣的流場控制方程。對于第1類方法,由于界面兩邊的物質屬性差異較大,過大的密度比和狀態方程差異,容易在界面處產生非物理震蕩,并且由于水的剛性很強,聲速較大會導致計算時間步長過小,降低計算效率。第2類方法在處理氣體中有激波存在以及高速問題時會遇到困難,通常用來模擬低速多介質運動,如水中氣泡上升,靜止水滴下落等。采用第3類方法時,充分考慮到水/氣物理屬性差異,可以處理氣體中有激波等高速運動問題,問題的關鍵在于如何給定合適的水/氣界面邊界條件,使得界面邊界條件能更合理的反映出水/氣界面處真實的流動狀態。

R.P.Fedkiw 等[2]提出了用 GFM 方法來處理多介質問題界面邊界條件,最初只是用來處理氣氣多介質問題,有效的 抑制了 在界面 處產生 的非物 理震蕩 。T.G.Liu等[3-4]在 其 基 礎 上 發 展 了 MGFM 方 法(modified ghost fluid method),通過在界面處定義 Riemann問題,并利用 Riemann問題的解定義了界面處虛擬流體點的速度和壓力值,通過等熵修正定義了真實流體點和虛擬流體點的密度值,MGFM 方法可以有效的處理各種氣氣多介質問題,但是在處理類似激波阻尼的臨界問題時會遇到困難;C.W. Wang等[5-6]提出了 RGFM 方法(real ghost fluid method),利用 Riemann問題的解不 僅 定 義 了 虛 擬 流體點的值,而且對真實流體中點的密度,速度和壓力都進行了更新,使得在界面處的流體狀態滿足接觸間斷的性質,得到了更加準確的界面邊界條件。利用上述幾種GFM 方法在處理水/氣多介質問題時,把水作為可壓縮流體處理,狀態方程采用 Tait方程。而當將水作為不可壓流體,氣體作為可壓縮流體處理時,在界面處需要采用更加合理的界面邊界條件,才能滿足由于界面兩邊不同的控制方程和物質屬性差異對界面處真實流動 狀 態 的 影響。針 對 水/氣多 介 質 問題,R.Caiden 等[1]提 出 了 new GFM 方法,界面處的速度采用界面附近水的速度值,而將界面附近氣體的壓力直接作為界面處的壓力。數值結果表明,該方法在一定程度上可以有效地捕捉物理現象。

new GFM 方法給出的界面邊界條件過于簡單,并不能充分反應出界面處真實的流動狀態,特別是當界面兩邊水/氣的物質屬性差異很大時。由于流體的流動性質很大程度上是物質屬性相關的,因此界面處的流動狀態通常非線性的依賴于界面兩邊流體的物質屬性。在文獻[7-9]中有以下結論:在有限區域內,在馬赫數趨于零的情況下,無黏不可壓縮方程是歐拉方程的收斂極限。在這個結論的基礎上,考慮在界面處定義新的水/氣 Riemann問題,其中可以把水視為聲速趨于無窮大的可壓縮流體,并求解Riemann問題的解,從而得到界面處的真實流體狀態。在本文中,利用得到的 Riemann問題的解分別采用 MGFM 與RGFM 這2種方法定義水/氣界面邊界條件。通過一維數值模擬實驗,以及與利用new GFM 方法的計算結果比較分析,2種定義界面邊界條件的方法均能有效的捕捉到激波與界面的位置,同時分析3種方法的計算結果的差異及其原因。

1 控制方程

1.1 可壓縮流體控制方程

一維無黏可壓縮流體控制方程為歐拉方程:

式 中 :ρ為 密 度 ,u為 速 度 ,p為 壓 力 ,為 體 積 總 能 ,e為 質 量 內 能 。 理 想 氣 體 狀 態 方 程 為,其 中

1.2 不可壓縮流體控制方程

采 用 投 影 法 求 解 不 可 壓 控 制 方 程,對 式(2)關 于x方 向 求 偏 導 得 到,由 于 密 度 為 常 數,則有。 表 明p為 常 數,那 么 在 不 可 壓 流 體 區 域 內 的 壓 力 滿 足 線 性 關 系 式,即 有,其 中xpr和pl分 別 為 不 可 壓 流 體 區 域 左 右 兩 個 端 點 上 的 的 壓 力 值,l為 不 可 壓 區 域 的 長 度 。

1.3 Level set方 程

通過求解level set方程隱式追蹤界面:

Level set函 數初始 化為距 離符號 函數,空間方向離散采用五階 WENO 格式[10],時間方向離散則采用三階 TVD Runge-Kutta方法。

2 界面邊界條件

在 new GFM 方法[1]中,只是簡單地直接取界面附近真實氣體的壓力作為界面處的壓力 值,界面的速度取臨近界面處水的速度值。事實上界面處的流動狀態都是與物質相關的,界面處的流動狀態非線性的依賴于界面兩邊介質的物質屬性,由于水/氣的物質屬性差異很大,應考慮到其對界面處真實流動狀態的影響。在 MGFM 和 RGFM 方法中,水作為可壓縮流體處理,采用剛性較強的 Tait方程作為其狀態方程,通過在界面處定義氣氣 Riemann問題,求解給出了界面處的流體狀態,體現了界面兩邊水/氣的物質屬性對于界面處流動狀態的影響。但是考慮到水作為不可壓流體時,具有密度不變等性質,不能求解氣氣 Riemann問題,上述方法則無法成立。當將水作為不可壓流體處理時,在文獻[9]中,有以下結論:隨著流體馬赫數趨向于零,可壓縮歐拉方程組的解會收斂到無黏不可壓流體方程組的解。另外,R.Agemi[7]證明了在有限區域內,當馬赫數趨向于零時,無 黏可壓縮流體會 收 斂 為相應的無黏 不 可壓流體。基于以上結論,考慮將不可壓流體作為聲速無窮大的可壓縮流體處理,在這一基礎上,定義新的針對水/氣的 Riemann問題,并求解出 Riemann問題的解,利用其定義更合理的水/氣界面邊界條件。

2.1 水/氣 Riemann問題

在 水/氣 界 面 處 ,定 義 Riemann問 題 ,初 始 條 件ρl,ul,pl和ρr,ur,pr分 別 為 界 面 左 右 兩 側 的 密 度 ,速度和壓力。為不失一般性,假設左側為氣體,右側為水。根據以上結論:無黏不可壓流體是可壓縮流體在聲速趨向于無窮大時的收斂極限。在求解水/氣 Riemann問題時可以將水作為聲速趨向于無窮大的可壓流體處理。可壓縮流體中的 Riemann問題可以用雙激波近似方法求解,雙激波近似求解 Riemann問題有以下關系式:

2.2 定義界面邊界條件

在處理多介質問題時,GFM 通過定義虛擬流體點來給定合適的界面邊界條件。首先,在水/氣界面處建立求解水/氣 Riemann問題,利用得到的 Riemann問題的解,定義界面處真實流體點和虛擬流體點的流體狀態值,從而給定合適的界面邊界條件。在得到的水/氣 Riemann問題解的基礎上,本文中采用兩種方式定義界面邊界條件,分別為 MGFM 方法和RGFM 方法。

如圖1所示,界面在網格點i與i+1之間,氣體在界面左側,水在界面的右側,取網格點和2的流場狀態分別作為 Riemann 問 題的 初 始 左 右 狀 態值,求 得 相 應 的 解。接 觸 間 斷 左 側 狀 態 為,右 側 狀 態 為。 在 MGFM 方 法 中[3],定 義 氣 體 界 面 邊 界 條 件 ,利 用 得 到 的 Rie-mann問題的解左側狀態 定義網格點i+1上的虛擬流體點的速度和壓力,網格點i+2和i+3處的速度和壓力取當地真實流體點的值,在求得界面左側熵值的基礎上利用等熵修正得到網格點i上真實流體 以 及 網 格 點i+1、i+2、i+3 上 出 虛 擬 流 體 的 密 度 值 ,上 述 則 給 定 了 氣 體 的 界 面 邊 界 條 件 。 可 以 采 用類似方法定義水的界面邊界條件,值得注意的是,Riemann問題的解中界面右側狀態中的密度也是保持不變的,滿足水的密度不變的性質,虛擬點上的密度都為水的密度。

圖1 利用 MGFM方法對氣體定義界面邊界條件Fig.1 The definition of interface boundary condition for gas by MGFM method

如 圖 2 所 示 ,在 RGFM 方 法 中[6],對 于 氣 體 ,得 到 的 Riemann 問 題 解 的 左 側 狀 態直接定義了網格點i+1、i+2和i+3上的虛擬流體狀態及網 格 點i上 的真實流體 狀 態,上述即給定 了 氣體的界面邊界條件。對于水可進行類似操作,同樣的,虛擬點上的密度都為水的密度。

上述2種方法給出了水和氣體的界面邊界條件,在這個基礎上,分別采用不同的離散方法求解水和氣的控制方程,并通過求解level set方程更新level set函數隱式追蹤界面。

圖2 利用RGFM方法對氣體定義界面邊界條件Fig.2 The definition of interface boundary condition for gas by RGFM method

3 數值模擬

對一維水/氣多介質問題進行數值模擬,時間方向上采用三階 TVD Runge-Kutta方法,空間上用三階ENO格式離散歐拉方程組,對于不可壓方程組采用標準的二階投影算法離散求解。

3.1 水在空氣中高速運動

在長度為1 m,網格點為200的區域內,中心處有長度為0.2 m 的水被周圍的靜止的氣體包圍,其中氣 體 的 狀 態 為γ=1.4,ρ=1.226 kg/m3,u=0,p=100 kPa,水 的 狀 態 為ρ=1 000 kg/m3,u= 100 m/s,p=100 k Pa。水在靜止的空氣中突然以一定速度向右運動,導致水右側的氣體中產生激波,在水的左側氣體中產生稀疏波。圖3~6所示為在時間t=0.75 ms時,利用 MGFM 和 RGFM 和 newGFM 這3種方法[1]計算得到的流場狀態結果對比圖。

由圖3~6中可以看出,3種方法都能較準確地捕捉到界面和激波的位置,在激波處,由于密度變化過小而很難在圖中體現。3種方法計算結果對于界面位置的捕捉結果基本相同,這是由于水/氣 Riemann問題的解給出的界面處的速度為界面附近水的速度,這一結論與new GFM 方法中界面的速度直接取界面附近水的速度相同。而 MGFM 方法和 RGFM 方法計算結果中,激波位置比new GFM 方法領先1個網格步長,這是由于界面處的壓力是由水/氣 Riemann問題得到,沒有簡單的取界面附近的氣體的壓力,充分考慮了界面兩側水/氣物質屬性差異對界面處流動狀態的非線性影響。此外,通過網格加密方法對 RGFM 方法進行了收斂速度測試,在200、400 和 800 個網格點上計算得到在時間t= 0.75 ms時水的速度依次為99.693、99.689 和 99.687 m/s,計算得到收斂速度為0.992 7。

圖3 密度1 000 kg/m3的水在空氣中向右運動時流場密度Fig.3 Density profile of the water movement in air while water density is 1 000 kg/m3

圖4 密度1 000 kg/m3的水在空氣中向右運動時流場速度Fig.4 Velocity profile of the water movement in air while water density is 1 000 kg/m3

圖5 密度1 000 kg/m3的水在空氣中向右運動時流場壓力Fig.5 Pressure profile of the water movement in air while water density is 1 000 kg/m3

圖6 流場壓力細節對比圖Fig.6 Detail comparison of pressure profile

為了進一步考察水的密度對于計算結果的影響,將上述算例中的水密度變為10 kg/m3,圖7~10給出了在時間t=0.75 ms時,3種方法得到的計算結果對比。

對于界面位置的捕捉,3種方法計算的結果基本相同,在激波處,密度變化明顯可以分辨。由于水的密度變小,更容易被空氣減慢速度,使得在水與右側的激波間,和水與左側的稀疏波之間都產生1個二次稀疏波。在這2個二次稀疏波中,new GFM 方法計算結果的速度大于 RGFM 和 MGFM 方法的計算結果,激波位置落后1個網格步長,MGFM 和 RGFM 方法只在稀疏波區域有細微差異,在區域其他點上基本相同。利用網格加密方法對 RGFM 方法進行收斂速度測試,在200、400和800個網格點上計算得在t=0.75 ms時,水的速度分別為73.589、73.315和73.171 m/s,計算出收斂速度為0.953 1。

3.2 激波與水/氣界面相互作用

圖7 密度10 kg/m3的水在空氣中向右運動時流場密度Fig.7 Density profile of the water movement in air while water density is 10 kg/m3

圖8 密度10 kg/m3的水在空氣中向右運動時流場速度Fig.8 Velocity profile of the water movement in air while water density is 10 kg/m3

圖9 密度10 kg/m3的水在空氣中向右運動時流場壓力Fig.9 Pressure profile of the water movement in airwhile water density is 10 kg/m3

圖10 流場壓力細節對比圖Fig.10 Detail comparison of pressure profile

在長度為1 m,網格數 為200 的區域 內,有 狀 態 為ρ=1 000 kg/m3,u=0,p=98.067 k Pa,長 度 為0.2 m的 水 在 區 域 中 心 ,周 圍 的 空 氣 狀 態 為ρ=1.58 kg/m3,u=0,p=98.067 kPa,在0.1 m 處 的 空 氣 中有一向右運動的激波,其 中波 后 狀 態 為ρ=2.124 kg/m3,u=89.98 m/s,p=148.407 k Pa。 圖 11~14所示為在時間t=1.75 ms時,分別采用 MGFM、RGFM 和 new GFM 這3種方法的激波與水/氣界面相互作用的計算結果對比圖。激波與水/氣界面相互作用時,由于水的剛性很強會產生1個反射激波,和1個很弱的入射波,入射波以極快的速度穿過水的區域,進入右側的空氣,入射波由于強度太小而很難在圖中分辨。MGFM 和 RGFM 方法的計算結果中,激波的位置領先new GFM 方法中激波位置1個網格步長。由于入射波的作用,水滴產生1個向右的速度,同時由于水密度較大,導致產生的速度過小。通過網格加密方法對 RGFM 方法進行收斂速度測試,在200、400和800個網格點上計算得到在時間t=1.75 ms時水的速度分別為0.531、0.538和0.542 m/s,計算得到收斂速度為0.9835。

將上述算例中的水密度變為ρ=10 kg/m3,圖15~18所示為采用3種方法計算得到結果對比圖。由結果可以看出,激波與水/氣界面相互作用以后,產生1個向左的反射激波,在這個反射激波與水之間形成了1個二次稀疏波,于此同時,向右運動的入射波以極快的速度穿過水進入右側的氣體中形成一個稀疏波,這一結果也符合波在水中的傳播速度無窮大的性質。

圖12 水密度1 000 kg/m3時激波與水/氣界面相互作用后的流場速度Fig.12 Velocity profile of shock impact with water-gas interface while water density is 1 000 kg/m3

圖13 水密度1 000 kg/m3時激波與水/氣界面相互作用后的流場壓力Fig.13 Pressure profile of shock impact with water-gas interface while water density is 1 000 kg/m3

圖14 流場壓力細節對比圖Fig.14 Detail comparison of pressure profile

圖15 水密度10 kg/m3時激波與水/氣界面相互作用后的流場密度Fig.15 Density profile of shock impact with water-gas interface while water density is 10 kg/m3

圖16 水密度10 kg/m3時激波與水/氣界面相互作用后的流場速度Fig.16 Velocity profile of shock impact with water-gas interface while water density is 10 kg/m3

由圖15~18中的對比結果可知,在激波位置方面,MGFM 和 RGFM 方法的計算結果仍然比new GFM 方法的計算結果要領先,在二次稀疏波范圍內,new GFM 方法的計算結果中,速度比另兩種方法的結果稍小,而壓力比另外2種方法的結果大。

由于 RGFM 方法中,不僅用 Riemann問題的解定義了界面處虛擬點的狀態值,同時也修改了界面處真實流體的速度和壓力,而 MGFM 方法中,只是修改了虛擬流體點的速度和壓力,對于界面處真實流體點只是通過等熵修正了密度,并沒有修改壓力和速度,因此在得到的結果中,RGFM 方法中水的速度相對另外兩種方法要略大。通過網格加密方法對 RGFM 方法進行了收斂速度測試,在200、400和800個網格點上計算得到在時間t=1.75 ms時水的速度分別為40.117、40.431和40.589 m/s,計算得到收斂速度為0.9922。

圖17 水密度10 kg/m3時激波與水/氣界面相互作用后的流場壓力Fig.17 Pressure profile of shock impact with water-gas interface while water density is 10 kg/m3

圖18 流場壓力細節對比圖Fig.18 Detail comparison of pressure profile

4 結 論

針對水/氣多介質問題,考慮將水作為不可壓流體,氣體作為可壓縮流體,在基于不可壓無黏流體是可壓縮流體在馬赫數趨向于零的收斂極限的結論上,定義并求解了水/氣 Riemann問題,Riemann問題的解反映出界面處真實的流動狀態,并采用 MGFM 以及RGFM 這2種方法給出了更合理的反應水/氣物質屬性差異的界面邊界條件。通過一維數值算例結果以及與相應的new GFM 方法計算結果的對比分析,證明了算法的有效性及合理性,同時對 RGFM 方法進行了收斂速度測試,充分考慮了水的不可壓縮性質以及氣體的可壓縮性質,適合針對水/氣具有不同的運動特性的水/氣多介質問題,例如空氣處于高速運動或者有強激波存在而水處于低速運動狀態,此時不能被忽略水的不可壓性質,同樣的適合于需要考慮水中波的傳播速度等問題,如水下爆炸。

[1]Caiden R,Fedkiw R P,Anderson C.A numerical method for two-phase flow consisting of separate compressible and incompressible regionsg[J].Journal of Computational Physics,2001,166(1):1-27.

[2]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.

[3]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.

[4]Liu T G,Khoo B C,Wang C W.The ghost fluid method for gas-water simulation[J].Journal of Computational Physics,2005,204(1):193-221.

[5]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.

[6]王春 武,趙 寧.基于求 解 Riemann問題的界面 處理方法[J].計算物理 ,2006,22(4):306-310. Wang Chun-wu,Zhao Ning.An interface treating method based on Riemann problems[J].Chinese Journal of Computational Physics,2006,22(4):306-310.

[7]Agemi R.The incompressible limit of compressible fluid motion in a bounded domain[C]∥Proceedings of the Japan Academy Series A:Mathematical Sciences.1981.

[8]Steve S.The compressible Euler equations in a bounded domain:Existence of solutions and the incompressible limit [J].Communications in Mathematical Physics,1986,104(1):49-75.

[9]Asano K.On the incompressible limit of the compressible Euler equation[J].Japan Journal of Applied Mathematics,1987,4(3):455-488.

[10]Jiang G S,Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics, 1996,126(1):202-228.

Interface treating methods for the gas-water multi-phase flows

Xu Shuang1,Zhao Ning1,Wang Chun-wu2,Wang Dong-hong2
(1.College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics, Nanjing 210016,Jiangsu,China;
2 College of Science,Nanjing University of Aeronautics and Astronautics, Nanjing 210016,Jiangsu,China)

A new interface treating method is presented for the compressible-incompressible gas-water multi-phase flow.The Riemann problem is constructed at the compressible gas-water interface,and then solved according to the hypothesis that the sound speed tends to infinity in the water.The solution of Riemann problem provides the fluid states for compressible gas and incompressible water at the interface.Those states can then be used to define the interface boundary condition by coupling the ghost fluid method.The level set method is employed to track the interface.The numerical examples of one-dimension case are given in this paper,furthermore,several comparisons are made with other results to verify the algorithm.Numerical results show that the provided algorithm can capture the discontinuities accurately,which demonstrates the robustness and efficiency.

mechanics of explosion;Riemann problem;ghost fluid method;level set;gas-water interface;multi-phase flows

O382.1國標學科代碼:13035

:A

10.11883/1001-1455-(2015)03-0326-09

(責任編輯 王易難)

2013-11-14;

2014-02-28

國家 自然科學基 金項目(11271188,91130030);北京理工大學爆炸科學與技術國家重點實驗室開放基金項目(KFJJ11-4 M)

徐 爽(1984— ),男,博士研 究生,shuangxu@nuaa.edu.cn。