一種求解多介質流問題的正保護數值方法

陳凌蕙，徐 偉

（南昌航空大學 數學與信息科學學院，南昌 330063）

引 言

在對多介質流進行數值模擬時，有時會遇到流體的密度、壓力或內能非常小的情況，如爆炸、高速射流及慣性約束核聚變等問題。此時用傳統的高精度格式，如WENO，DG等進行求解，常會由于數值誤差而得到負的密度或壓力值，而這顯然是不符合實際情況的，同時這也會使得算法不穩定，從而導致計算失敗。針對這類問題，目前已有一些研究工作，發展出了一些具有正保護功能的數值算法，但大多是針對單介質流體的，而對于多介質流的相關工作則相對較少。這是因為，多介質流體在介質界面兩邊的密度與壓力等狀態常常會相差巨大，其所使用的狀態方程也不盡相同，甚至不同介質使用的控制方程有時也不一樣，這些都使得構造求解多介質流動問題的正保護數值算法要困難得多。解決這類問題的一個簡單的方法是強行將負的密度或壓力值等賦值為零，但這會使得算法的精度和守恒性都無法得到保證。Einfeldt在早期給出了一類求解可壓縮Euler方程的，具有正保護功能的Godunov格式[1]，但是只有一階精度。Perthame則分別給出了具有正保護功能的一階Lax-Friedrichs格式[2]和二階Boltzmann格式[3]。然而要想將這些格式中給出的正保護方法在不破壞精度的條件下，推廣運用于任意高精度格式是非常困難的。Zhang和Shu在上述Lax-Friedrichs格式的基礎上給出了一種更易于推廣到高精度的正保護格式構造方法[4]，并給出了判斷格式是否具有正保護性質的充分條件，但這種方法只是適用于理想氣體。Cheng和Shu則結合HLLC方法和Lagrange方法，給出了一種針對多介質流動問題，且具有正保護功能的數值算法[5]。上述方法進行正保護的大致思路都是通過構造正保護限制器，進而利用限制器對計算出現的負密度負壓力等值進行修正，同時保證不破壞格式原有的精度及守恒性。

另外，在求解多介質流問題時，若直接使用傳統的高精度格式，由于介質界面兩邊流體狀態方程的不同，數值結果在界面附近會出現非物理振蕩或計算失敗。針對該問題，Fedkiw等給出了一種GFM方法[6]，這種方法通過在界面附近定義虛擬流體，將多介質問題轉化為多個單介質問題，然后再使用已有的高精度方法進行計算。最初的GFM方法在計算時，將虛擬流體的壓力與速度值取為當地數值，再由熵修正得到密度。但是在計算水-氣多介質流時，該算法的穩定性較差，從而產生了一些變形，如nGFM法[7]、mGFM法[8]以及rGFM法[9]等。其中mGFM法與rGFM法都是通過在介質界面附近構造黎曼問題，并利用黎曼問題的解來對界面附近的單元進行特殊處理。不同的是，mGFM法只對虛擬流體中的單元進行定義，而在真實流體中只對密度進行等熵修正。這種做法實際上使得真實流體與虛擬流體在界面附近的單元沒有得到同步更新，從而在計算某些問題時還是會出現虛假振蕩。而rGFM法則利用黎曼問題的解同時對真實流體與虛擬流體進行更新，能更真實的反映流體的特性。

由上述可知，在計算多介質流問題時，有兩處需要進行正保護，一是在界面處計算黎曼問題的時候，一是對各單介質進行計算的時候。本文基于正保護限制器[4]與正保護HLLC方法[5]，給出了一種適用于氣-氣及水-氣多介質問題的正保護數值算法，并通過多個算例驗證了算法的有效性。

1 方 程

1.1 控制方程

一維無粘可壓縮Euler方程：

其中U=(ρ,ρu,E)T，F(U)=(ρu,ρu2+p,(E+p)u)T，E=ρe+ρu2/2，式中的ρ是密度，u是速度，p是壓強，e是單位質量的內能，E是單位體積的總能。

二維無粘可壓縮Euler方程：

其中

式中的ρ，e與E的含義與一維情況一樣，而u是v則分別 表示x和y方向上的速度。

1.2 狀態方程

為了求解Euler方程，需要給出流體的狀態方程。當流體為理想氣體時，其方程為：

而水的狀態方程有多種，本文取

其中γ與B為熱力學常數，將在具體算例中給出。不難看出氣體與水的狀態方程可以表示成統一的形式。此外，式（4）也可寫為p+B=(γ?1)(ρe?B)。

2 正保護數值方法

2.1 正保護有限體積法

對于一維情況，首先取點集G：

利用Jensen不等式易證G為凸集。而所謂正保護格式，即要使得當Un∈G時，能夠有Un+1∈G。

考慮有限體積法的一般形式：

其中，m=ρu。接著取單元Ij上的Gauss-Lobatto點集

且有2N?3≥k。

由文獻[4]可知，對于有限體積格式（5），若在tn時刻，?α,j，有以及。則在CFL條件：

可得tn+1時刻的∈G。

根據上述結論，下面將給出正保護限制器對單元Ij上的重構多項式qj(x)進行修正，使得修正后的滿足：

修正將分兩步進行，具體如下：

首先，密度修正：在各單元Ij上，將密度函數ρj(x)修正為，即取

其中

于是

本文采用WENO方法進行空間離散，但由于WENO格式與DG格式不同，它不給出重構多項式的具體形式，而是直接得到點值與。因此為了使用上述限制器對重構多項式進行修正，可以按文獻[10]中的方式得到重構多項式，再對其進行修正，繼而利用修正后的多項式再計算出所需的點值與。而時間離散則使用三階Runge-Kutta時間離散，式（5）實際上就是其中的第一步。由于其后兩步都是凸組合的形式，因此格式仍然具有正保護性質。

對于二維情況，不妨設區域被劃分為若干矩形網格單元Iij。在每個小矩形網格單元中，取如下點集：

與一維情況類似，若在tn時刻，?α,β,i,j，有。則在CFL條件：

可得tn+1時刻的。其中qij(x,y)為單元Iij上 各守恒量的k次近似多項式向量(ρij(x,y),mij(x,y),ni j(x,y),Eij(x,y))T。因此，二維情形下正保護限制器的構造及修正方式也與一維情形均基本相同。

2.2 正保護界面處理方法

本文采用rGFM方法來重定義界面兩側的流體狀態。使用該方法需要在界面處求解黎曼問題，因此也需要一種具有正保護功能的算法。

首先考慮一維情況。假設界面位于單元j與j+1之間，如圖1所示。

圖1 一維界面邊界條件定義

于是構造如下黎曼問題：

對于求解黎曼問題，已有許多成熟的求解器。本文采用基于雙激波近似的HLLC求解器，如圖2所示。

圖2 雙激波近似黎曼問題

圖中的u?是介質界面的速度，sL與sR分別是左右激波的波速，是該黎曼問題的解。以左激波為例，由于流體在激波兩側的狀態滿足Rankine-Hugoniot條件，則有

其中

而對于右側激波則可以得到類似的式子。又由于在介質界面處的速度和壓力連續，即有。于是由上述等式可得黎曼問題的近似解：

可以證明，當取

sL=min(uL?cL,?),sR=max(uR+cR,+),

而對于二維多介質問題，則需要沿介質界面的法向構造并求解黎曼問題[9]。由于此時得到的黎曼問題仍是一維形式，因此可以使用上面給出的正保護方法進行求解。與一維問題相比，二維問題中介質界面的追蹤要復雜的多，本文采用Level Set方法[11-12]對界面進行追蹤。

3 數值實驗

下面將通過幾個數值算例來驗證算法的有效性。這些算例由于初始條件比較極端，因此會出現低密度低內能區域，傳統高精度方法容易算得負值，從而導致程序崩潰。

算例1 氣?氣雙稀疏波問題

該問題的計算區域為[?1,1]，介質界面位于區域中心。其左側流體的初始狀態為：

右側的初始狀態為：

在該問題中，介質界面兩側的流體同時向兩側反向運動，從而產生兩個強稀疏波，使得中心區域的壓力大幅下降，幾乎形成真空。圖3給出了網格數為300，t=0.6時的數值結果。

圖3 t=0.6時的數值解

算例2 氣?液激波管問題

這是一個氣?液多介質流問題，計算區域為[0,1]，介質界面位于x=0.3處。界面左邊氣體的初始狀態為：

右邊液體的狀態為：

該問題是一個極端條件問題，兩側流體的初始壓力p與狀態方程中的參數B相差非常大。隨著時間推進，將形成一個向左的激波和一個向右的稀疏波，圖4給出了t=0.000 24時刻，網格數取為300的數值結果?？梢钥吹?，即使兩種流體的物理特性及初始條件都相差較大，本文給出的正保護算法依然能準確的捕捉到各種間斷，并且介質界面非常清晰。

圖4 t=0.000 24時的數值解

算例3 高壓氣泡爆炸問題

這是一個二維氣?氣多介質流問題。計算區域為[?1,1]×[?1,1]，其中充滿靜止的理想氣體，中心處有一個半徑為0.3的靜止高壓氣泡。氣泡內部的流體初始狀態為：

外部流體的初始狀態為：

該問題中，氣泡爆炸會產生一個向內的稀疏波，使得內部的密度與壓力值急劇減小。圖5、圖6給出了網格數為200×200，t=0.000 7時刻的數值結果?？梢钥吹?，本文算法保證了結果的正確性，且能準確的追蹤到介質界面。

圖5 t=0.000 7時的等值線圖

圖6 沿中軸線y=0的數值解

4 結 論

針對一些極端條件下的多介質流問題，會因算法等誤差導致出現不合理負值的情況，本文給出了一種具有正保護功能的求解方法。該方法首先使用具有正保護功能的雙激波近似黎曼問題求解器求解黎曼問題，并用得到的近似解對界面附近的真實和虛擬流體狀態進行重定義。接著在分別對各單介質流體進行求解時，將針對理想氣體的正保護WENO算法推廣運用于液體，從而得到了一種高精度正保護數值算法。最后通過數值算例驗證了該算法能有效保正結果的正確性，且分辨率高，能準確捕捉到各類間斷。