球坐標系下多介質混合物模型的數值模擬*

吳宗鐸，趙 勇，嚴 謹，宗 智，高 云

（1.廣東海洋大學海洋工程學院，廣東 湛江 524088；2.大連海事大學交通運輸工程博士后流動站，遼寧 大連 116026；3.大連理工大學船舶工程學院，遼寧 大連 116024；4.西南石油大學油氣藏地質及開發工程國際重點實驗室，四川 成都 610500）

多介質的激波間斷問題一直以來都是計算流體力學方面的一個備受關注的問題，目前廣泛應用的一個處理方法是利用Level Set 函數追蹤氣水界面，同時結合近似Riemann 解下的虛擬流體質點完成對沖擊波參數的計算[1]。雖然得到的結果較理想，但是計算過程復雜，且計算依賴近似Riemann 解，無法在JWL (Jones-Wilkins-Lee)及Mie-Grüneisen 等復雜形式狀態方程下展開更深入的研究[2-3]。

相比之下，利用多介質混合模型，可以將流場視為一個整體并用體積分數或質量分數來區分流體，從而在不直接求近似Riemann 解的前提下完成計算。多介質混合物模型早期由Abgrall 提出[4-5]，通過非守恒方程來計算由氣體或低壓縮度的固液體組成的流場[5-6]。隨后，Abgrall 將該混合物模型擴展成一個七方程模型，能適應一般形式狀態方程的流場[7]。柏勁松等[8]也利用此類多介質混合物方程來模擬兩種以上介質的沖擊過程，且配合等效方程可以適應一些復雜形式的狀態方程[9]。梁珊等[10]和Liang 等[11]將Abgrall 的七方程模型與多GPU 的并行計算相結合來提高計算精確度。為了讓算法更好的適應并行計算，劉娜等[12]將譜體積方法運用到多介質混合物模型中，提升了計算效率。

雖然該多介質混合物模型在采用直角坐標系處理平面波的問題時，取得不錯的進展，但是對于常用于處理球面波的球坐標系下的多介質問題，仍然具備一定的局限性。第一，多介質混合物模型建立在速度和壓力在界面處維持平衡的基礎，僅對平面波適用，但是球面波的壓力和速度都會隨著傳播半徑的增加而衰減；第二，球坐標系下，原點位置處為沖擊波運動的奇點位置，必須作出合適的處理；第三，為了增加工程實用性，計算模型需要能適應復雜形式的狀態方程，而不僅僅是理想氣體方程等形式簡單的狀態方程。

考慮到以上因素，本文將直角坐標系下的基于Mie-Grüneisen 狀態方程下的多介質混合模型延伸到球坐標系下。為了壓力和速度不平衡問題以及坐標原點的奇點問題，在原計算模型下作了參數修正并用質量分數替換體積分數，并對奇點處數值差分作了特殊處理—利用相鄰網格點參數賦值給奇點參數但奇點處速度設定為0。最后，將通過數值算例驗證其數值穩定性和計算準確性。

1 直角坐標系下的計算模型

考慮多介質的可壓縮流體組成的混合流場，其運動形式可以用歐拉方程來描述：

式中：ρ、p、u、E 分別為密度、壓力、速度和單位質量能量[13]。單位質量的總能量為內能和動能的總和，且E=e+u2/2。

流場中介質的物理特性，可以用Mie-Grüneisen 狀態方程來表述：式中：Γ 為Mie-Grüneisen 系數，pref和eref為參考點的壓力和單位質量內能，與密度ρ 有關。

1.1 界面處熱力學關系

界面處，流體的壓力p 和速度u 均保持平衡[4]，此時p 與u 在界面處不發生變化，而密度ρ 與內能e 則出現數值間斷。因此?p/?x 和?u/?x 約為0，式(1)中的質量守恒與能量守恒可以分別表示為：

將式(2)代入式(4)，有：

由于Γ、pref、eref都只和密度ρ 有關的函數，且式(5)對于任意的表達式pref(ρ)和eref(ρ)均成立。這里首先假設pref、eref都為0，那么有[2,9]

同理，可依次類推得到關于pref和eref的非守恒形式的方程：

1.2 波面處熱力學關系

由于式(6)～(8)的推導過程需要基于p 和u 在界面處連續這一前提條件，但是在左右波面處，p 和u 均出現間斷面，式(3)和(4)不再成立。此時，對(6)作如下考慮：

因此，在左右波面處雖然p 和u 均間斷，但由于在同類介質中1/Γ 等參數的導數保持連續，因此，式(6)依然成立。同理，式(7)和(8)在左右波面處也成立。

1.3 輸送方程

對于多介質流場，必須利用輸送方程來捕捉界面的運動。流場中，某種介質的界面隨時間變化的運動方程，稱為輸送方程。一般利用體積分數zi來建立介質的輸送方程：

式中：i 為介質的種類序號。非守恒變量1/Γ、pref/Γ、ρeref均為獨立變量，但它們對的ρ 偏導數可表達成關于體積分數zi的加權疊加。這樣，綜合式(1)、(6)～(8)和(10)，可完成直角坐標系下的多介質流場的求解。

2 球坐標系下的計算模型

球坐標系下，則需要考慮更多的問題。

第一，在采用直角坐標系下計算平面波問題時，p 與u 在界面處維持不變且偏導數?p/?x 與?u/?x 基本為0。該數值特性一直從界面處向左右延伸，直到波面處才形成間斷，如圖1 所示。但是采用球坐標系處理球面波問題時，當沖擊波沿著半徑r 向外擴散時，傳播的形式則不再是平面沖擊波，而是球面沖擊波。此時，雖然p 與u 依然在界面處連續，但是球面波的擴散過程中p 與u 會隨著半徑增加而發生衰減，并非保持不變。因此，球面沖擊波在界面處的偏導數?p/?r 與?u/?r 不再是0，而是一個穩定的負數。這種情況下，式(3)～(4)非守恒性方程的建立的理論基礎并不十分準確。

第二，球坐標系下的內外介質產生相互作用時，內部的稀疏波(或激波)匯聚到一起后又反向朝著外部擴散，原點處便形成了一個奇點。奇點處的稀疏波形成了獨特的反射機制。因此奇點的物理參數變化需仔細考慮。

圖 1 平面沖擊波界面處p 的數值特征Fig.1 The numerical property of p at the interfacefor plane shock wave

2.1 球坐標系下控制方程

這里將球坐標系下的控制方程表示成守恒變量U 與通量F 的關系為[13-14]:

其中：

式中：系數α=2。式(11)可以改寫成：

其中：

此時，將非守恒形式的方程組添加到式(11)中，即可得到球坐標系下的控制方程。

2.2 熱力學參數修正

需要注意的是，非守恒形式的方程式(6)～(8)的推導，是以界面處于平衡狀態為前提的。在球坐標系下，由于沖擊波是沿著球面向四周擴散。隨著擴散范圍的增加，幾乎所有的沖擊波參數都會隨著半徑增加而降低。此時，p 與u 在界面處將呈現出持續的減小趨勢。這樣，式(6)～(8)無法直接使用到式(12)中。在直角坐標系和球坐標系下，界面間斷的關系分別為

考慮到三個非守恒變量1/Γ、pref/Γ、ρeref均為密度ρ 的函數，這里每隔一個時間步對它們做修正:

式中：ρi為介質i 的密度。而的其他參數如ρ、p、u、E 的計算，仍然借助當前時間步下的Γ、pref、eref。而式(6)～(8)則保留在系統方程中，并且每隔一個時間步，非守恒變量1/Γ、pref/Γ、ρeref都會得到修正。當時間步足夠密時，非守恒形式在球坐標系下帶來的誤差可以控制在一定的范圍內。

2.3 輸送方程變換

由于每隔一個時間步，都會需要利用ρi進行修正。這里用質量分數代替體積分數來構造輸送方程。利用式(12)中的質量守恒關系可得：

根據式(14)得到的質量分數yi可以很容易得到各介質的密度ρi(i=1, 2, ···)。而非守恒方程(6)～(8)中的偏導數，則轉換成與質量分數相關的加權疊加：

式中：φ、φ、ψ 為非守恒變量1/Γ、pref/Γ、ρeref偏導數，下角標1、2 代表流場中介質的序號。

相比式(10)給出的基于體積分數的輸送方程，基于質量分數的輸送方程能考慮到每一種介質的密度ρi，從而得到更為合理的偏導數φ、φ、ψ。

2.4 奇點處理

對于球坐標系來說，奇點是一個無法避免的問題。參考文獻[15]中對奇點的處理方式，奇點處質點的運動速度為0，即u=0，但是奇點處的質量和能量守恒關系依然成立：

展開后，將u=0 代入，有：

由于在奇點處主要涉及沖擊波的匯聚和發散，不涉及界面的運動，因此，可不考慮輸送方程。參數Γ、pref、eref由密度ρ 直接計算得到。

文獻[15]中，對偏導數?u/?r 采用了如下的計算方式：

式中：下角標1 和2 代表網格點序號。但考慮到本文的方程體系非守恒性較強，將u1和u2都替換為奇點處的速度0，這樣偏導數?u/?r=0。再根據式(17)計算得到奇點處的其他參數。

2.5 差分格式

在有限體積差分下，第j 個網格點的變量Uj，可按如下格式計算：

式中：Δt 和Δr 分別為時間和空間步長，上標(n)代表第n 步時間序號。對于物理量第j 和j+1 個網格點之間的通量和 源項分別為：

式中：角標j+1/2 表示通量中的參數，角標j 表示網格點的參數；在球面沖擊波中，通量表示沿半徑r 方向通過j 點的物理變量。

2.6 時間空間步長

時間步長Δt，滿足收斂的條件為：

式中：CCFL是為了保證迎風型計算格式收斂而設置的一個0～1 之間的系數，滿足收斂的條件為CCFL＜1,并且當CCFL取得越小時間步長越小，計算越精細, 這里取CCFL=0.4； uj和cj表示第j 個網格點的質點速度和聲速。聲速cj的表達式為：

式中：H=E+p/ρ。

3 計算結果

3.1 無因次參數的氣水作用算例

這里考慮一個無因次的氣水作用問題[16]。初始時刻，中心氣團的物理狀態為：p=8.3×103，ρ=1.27，u=0，外部水的物理狀態為：p=1.0，ρ=1.27，u=0。氣體和水的狀態方程分別為

式中：γ 和B 為描述水和氣體在低壓縮度時的熱力學參數。氣團初始半徑為R0。由于中心氣團的壓力和密度高于周邊水的。因此，高壓的氣團會在水中形成沖擊波并推動界面向外部擴散。對沖擊波波面和介質界面的運動位置的變化情況進行記錄（波面為水中壓力的最大值位置，界面位置為從外向內氣體的質量分數剛超過0.5 的位置），并繪制了沖擊波和界面隨時間變化的曲線，如圖2 所示。從圖2 可以看出，本文所得到的沖擊波和界面的運動位置與Liu 等[16]和Flores 等[17]的結果都比較接近。三個數據的計算結果，在氣水作用的初期差別不大。但是隨著沖擊波和界面向外擴散，本文的計算結果與Flores 等[17]的計算結果相比，誤差略偏大，但仍然在一個合理范圍內。

圖 2 沖擊波與界面位置隨時間的變化Fig.2 Time evolutions of position of shock wave and interface

圖 3 沖擊波到達3 倍R0 時，壓力與速度的分布Fig.3 The distributions of pressure and velocity when shock wave reaches 3R0

圖3 為沖擊波到達3 倍氣團初始半徑時候的壓力和速度曲線。此時，向內收縮的稀疏波已達到中心原點處。由于氣團已擴散開，中心處的壓力開始逐漸下降，但中心處的速度仍然為0。速度曲線中，由中心向外的第一個拐點為稀疏波的端點，第二個拐點則為界面所在位置。另外，從圖3 中可以看出，處于數值奇點的原點處，數值穩定性較好。證明本文所用的奇點處理方式能效果良好。

3.2 圓形藥球的水下爆炸算例

假設一TNT 的實心藥球質量為50 kg，將空間步長取為炸藥球半徑的1/100。利用改進后的多介質混合模型計算球坐標系下的沖擊波運動情況。

炸藥的狀態方程為JWL 狀態方程[3]：

式中：A、B、R1、R2、ω 均為常數，由圓筒試驗標定得到[3]；θ=ρ/ρ0，ρ0為物質的初始密度。水的多項式狀態方程為[18]：

式中：μ=θ-1；a1～a3, b0～b3為展開后的多項式系數。狀態方程(23)～(24)的參數取值見表1[19]。

表 1 狀態方程參數Table 1 Coefficients of EOS

本算例中，將計算范圍擴大至20 倍以上的藥球半徑。這樣，可以得到沖擊波運動到不同位置處的壓力，如圖4 所示。從圖4 可以看出，由于在偏導數φ、φ、ψ 的處理方式差異，兩種計算模型得到的壓力分布有差別。但都具備較好的數值穩定性。當沖擊波運動到較遠位置時，會出現多個壓力峰值。

在實際工程中，最受重視的是第一次的峰值。圖5 為沖擊波峰值壓力隨運動距離變化的曲線。作為對比，將給初始的體積分數模型[2]和改進后的質量分數模型同Zamyshlyayev 的經驗公式[20]曲線作了對比。對比發現，改進后的質量分數模型可以得到更為準確的計算結果。在壓力隨運動距離發生急劇下降的過程中，利用體積分數得到的計算結果，得到的結果比經驗公式的結果略微偏大。

圖 4 沖擊波到達不同位置處的壓力曲線Fig.4 Pressure curves at the time instants when the shock wave reaches different positions

圖 5 沖擊波在不同位置處的壓力峰值變化Fig.5 Peak value variation of pressure of the shock wave at different positions

圖6 為利用體積分數與質量分數計算模型分別得到的壓力時程曲線。從不同的計算結果中可以看出，利用改進后的質量分數模型，得到的計算結果好于體積分數模型。各個半徑處的壓力時程，質量分數模型均吻合得比較好。由于Zamyshlyayev 的經驗公式僅考慮沖擊波壓力最大峰值未考慮壓力的二次峰值，本文也僅模擬了沖擊波壓力的最大峰值。

圖 6 不同計算模型下壓力隨時間變化曲線Fig.6 Time evolutions of pressure for different numerical models

4 結 語

針對球坐標下多介質界面捕捉困難的問題，本文對原始的Mie-Grüneisen 多介質混合物模型作了熱力學參數修正、輸送方程變換、特殊奇點處理等多項數值改進。改進后的基于質量分數的模型，可以在常規狀態方程下準確捕捉界面和沖擊波的運動，且可以在適應JWL 等形式復雜的狀態方程。而適當的奇點處理也保證了數值計算在奇點處的穩定性。另外，數值算例證明，相比原始的體積分數模型，改進后的質量分數模型可以得到較為準確的計算結果。