一種簡單的精確捕捉接觸間斷的黎曼求解器

胡立軍， 杜玉龍

(1.衡陽師范學院 數學與統計學院，衡陽 421002；2.北京航空航天大學 數學科學學院，北京 100191)

1 引 言

基于Godunov型數值格式[1]的有限體積法由于其良好的守恒性和網格適應性，已經成為求解雙曲守恒律系統最具代表性的數值方法。通過求解局部黎曼問題來得到網格界面的數值通量是實現有限體積方法的關鍵步驟。從物理角度來看，采用迭代方法求解黎曼問題的精確解似乎是最合理的，但是由于其效率低且實現難度大，在實際問題的計算中精確黎曼求解器很少使用。因此構造性能良好的近似黎曼求解器一直都是計算流體力學研究的熱點問題。

過去幾十年，研究人員構造了許多不同性能的近似黎曼求解器。根據捕捉接觸間斷的能力，可以將其分為兩類。非全波求解器，如Rusanov格式[2]、HLLE格式[3]和HLL-CPS格式[4]，過高的耗散行為在計算中不能精確分辨接觸波或者剪切波；全波求解器，如Roe格式[5]、Osher格式[6]和HLLC格式[7]，在計算中能夠精確捕捉接觸間斷和剪切波。但是，在計算強激波問題時，這些低耗散的求解器會遭遇嚴重的不穩定現象，這大大限制了它們在高超聲速流動問題中的應用。

研究人員嘗試在保留全波求解器精確分辨接觸間斷優點的同時來消除它們的激波不穩定性，其中最流行的方法是根據當地流場在耗散格式和低耗散格式之間進行切換的混合方法[8-12]。盡管混合格式可以成功地抑制激波失穩現象，但是對于復雜的流動問題，不恰當的開關函數可能會影響接觸間斷的分辨率，并且兩種格式間的突然切換也可能會影響到格式的收斂速度[13]。此外，增加多維耗散也是抑制全波黎曼求解器激波失穩現象的一種常用方法，主要通過增加格式的剪切粘性[14,15]、構造旋轉黎曼求解器[16,17]或者法向速度重建技術[18]來實現。旋轉格式在每個網格界面處需要計算兩次數值通量，因此計算效率低。而單純地增加剪切粘性在某些情形下并不能完全消除激波異?，F象[19]。最近，Chen等[20]通過在動量通量中引入剪切粘性構造了一種具有良好激波穩定性的HLLC+格式。

此外提高非全波黎曼求解器的接觸捕捉精度也是構造性能良好的數值格式的一種策略。本文利用雙曲正切函數[21]來重構界面兩側的密度值，并且結合邊界變差遞減算法[22]來減小單波Rusanov格式耗散項中的密度差，從而提高格式分辨接觸間斷的能力。數值結果表明,本文構造的黎曼求解器比全波的HLLC格式具有更高的接觸分辨率和更好的激波穩定性。

2 背景知識

守恒形式的二維歐拉方程組為

(1)

式中

(2)

(3)

本文的比熱比γ取1.4。

采用有限體積法對式(1)進行空間離散，可以得到半離散方程為

(4)

式中Ui,j為守恒向量U的單元平均值，Fi +1/2,j和Gi,j +1/2為數值通量。式(4)可以改寫成ODE系統為

dU/dt=L(U)

(5)

2.1 時間離散

采用優化的三階TVD龍格-庫塔格式[23]來進行時間離散,

U(1)=Un+ΔtL(Un)

(6)

式中Δt為時間步長。

2.2 五階WENO格式

采用五階WENO格式[24]重構得到界面左右兩側的狀態值為

(7)

(8)

非線性權重系數定義為

(9)

(10)

式(7)的線性加權系數為

(11)

而式(8)的線性加權系數為

(12)

為防止分母為0，取ε=10-6。

利用WENO重構得到界面左右兩側的狀態值，然后利用近似黎曼求解器求解局部黎曼問題，從而得到每個網格界面處的數值通量。

2.3 Rusanov求解器

最簡單的單波Rusanov黎曼求解器[2]的數值通量函數為

(13)

波速S定義為

S=max(|uL-aL|,|uR-aR|,|uL+aL|,|uR+aR|)

(14)

2.4 HLLC求解器

全波的HLLC黎曼求解器的數值通量函數為

(15)

式中

(16)

波速估計為

(17)

3 一種精確分辨接觸間斷的Rusanov黎曼求解器

Rusanov格式的通量函數式(13)可以改寫為

(18)

式中右端第一項為中心差分項，D1/2為耗散項，其表達式為

(19)

(20)

(21)

為了計算方便，本文直接給出網格界面左右兩側密度重構值的計算公式為

(22)

式中

(23)

式中β為控制密度跳躍大小的參數，本文取β=1.6。

(24)

將由式(24)確定的界面兩側的密度值代入式(19)計算耗散項D1/2，從而得到一種低耗散的Rusanov格式(命名為LD -Rusanov，Low Diffusion Rusanov)。此外,本文為了提高非線性波的分辨率，采用Roe平均值來計算波速，

(25)

在計算多維問題時，采用算子分裂方法來逐維計算數值通量。因此本文所述的一維算法可以直接應用于每個坐標方向。雙曲正切重構作為一種代數方法，在進行多維計算時不會涉及到任何的幾何重構，因此實施起來比較簡單。此外，與HLLC格式通過修改HLL格式的波系結構來提高接觸波的分辨率不同的是，LD-Rusanov格式采用代數方法而不是修改原Rusanov格式的波系結構來提高接觸間斷的捕捉能力，因此其很好地保留了原Rusanov格式的激波穩定性。接下來，通過一系列的數值實驗來證明LD-Rusanov格式的高分辨率和強魯棒性。

4 數值結果

通過一些典型的數值算例來驗證LD -Rusanov格式的表現。

4.1 孤立的接觸間斷

首先考慮一個孤立的慢行接觸波，計算區域為[0,1]，其初始條件為

(26)

計算中使用網格數為100的均勻網格。從 圖1 可以看出，本文構造的LD -Rusanov格式不僅可以更加精確地捕捉該接觸間斷，而且還消除了原Rusanov格式在間斷附近出現的過沖。圖2展示了取不同β值的LD -Rusanov格式計算該算例時的計算結果?？梢钥闯?，β=1.3,1.6和1,9都可以得到精確的結果。為了定量分析不同β值對于結果的影響，采用式(27)來計算間斷的厚度，

(27)

圖1 孤立的接觸間斷在t =3時的密度分布

圖2 取不同β值的LD -Rusanov格式的密度分布

由表1可知，β值越大，間斷的厚度越小，其分辨率也會越高，但是過大的β會卷起平行于速度方向的界面[21]。本文發現1.2～2.0的β值可以得到令人滿意的結果。

4.2 一維爆轟波問題

計算區域為[0,1]，初始條件為

(28)

計算中使用網格數為200的均勻網格。該算例的解由兩個相互作用的爆轟波形成的復雜間斷結構構成。從圖3可以看出，相比原始的Rusanov格式和全波的HLLC格式，本文構造的LD -Rusanov格式能夠更加精確地捕捉各個波系。

圖3 一維爆轟波問題在t =0.038時的密度分布

4.3 Titarev-Toro問題

考慮Titarev等[24]提出的激波-熵波相互作用問題。計算區域[-5,5]均勻劃分成500個網格，初始條件為

(29)

該算例涉及到高頻振蕩的正弦波和激波的相互作用。從圖4可以看出，Rusanov格式的計算結果具有較大的數值耗散，而本文構造的LD -Rusanov格式得到了更加精確的解，在捕捉高頻波時比全波的HLLC格式具有更高的分辨率。

圖4 Titarev-Toro問題在t =5時的密度分布

4.4 二維爆炸問題

二維爆炸問題是一維Sod問題的推廣，計算該算例來檢驗格式捕捉不同波系的能力。計算區域[-1,1]×[-1,1]均勻劃分成101×101的正方形網格，初始條件為

(30)

計算時間為t=0.25。此時，該算例的解由一個向外運動的圓形激波和接觸面以及向中心運動的稀疏波構成。圖5展示了時間t=0.25時沿徑向y=0的密度分布?？梢钥闯?，本文構造的 LD -Rusanov 格式能夠精確地捕捉各個波系，特別是對接觸間斷的分辨率明顯優于原始的Rusanov格式和HLLC格式。

圖5 二維爆炸問題在t =0.25時的密度分布

4.5 移動的圓形接觸面問題

計算區域[0,1]×[0,1]劃分成100×100的均勻網格，其初始分布為

(31)

該算例的解由一個以勻速(u,v)=(1,1)運動的圓形接觸面構成。在t=0.3時刻該問題的精確解為

(32)

從圖6可以看出，相比于原始的Rusanov格式和全波的HLLC格式，本文構造的LD -Rusanov格式對于該圓形接觸面具有更高的分辨率。

圖6 移動的圓形接觸面問題的計算結果

4.6 奇偶失聯問題

本文數值實驗表明，在捕捉接觸間斷時LD -Rusanov格式比HLLC格式具有更高的分辨率。文獻[8-20]曾報道精確分辨接觸間斷的數值格式(如HLLC格式)在計算多維強激波問題時會出現嚴重的不穩定現象。接下來，計算幾個典型的強激波問題來檢驗LD -Rusanov格式的魯棒性。

首先,計算Quirk[8]提出的奇偶失聯算例來評估格式的激波穩定性。馬赫數為20的平面運動激波從左向右傳播，其初始位置為x=5。右側的波前狀態為(ρ,u,v,p)R=(1.4,0,0,1)，左側的波后狀態由Rankine-Hugoniot關系式計算得到。計算區域[0,600]×[0,20]分割成600×20的矩形網格，且y-方向的網格中心線上存在±10-3的奇偶小擾動。圖7展示了t=20時的密度等值圖，為了全面比較格式的激波穩定性，本文還展示了兩種穩定版本的HLLC格式(HLLCM[13]和HLLC+[20])的計算結果?？梢钥闯?，HLLC格式表現出了明顯的激波失穩現象，而LD -Rusanov格式和另外兩種穩定版本的HLLC格式均得到了穩定的激波面。

圖7 奇偶失聯問題的密度等值圖

4.7 穩態斜激波問題

接下來考慮穩態斜激波問題[26]。計算區域[0,50]×[0,30]劃分成50×30的正方形網格，沿著直線y=2(x-12)有一個穩態激波，其波前狀態為(ρ,u,v,p)L=(1,447.26,-223.5,3644.31)，波后狀態(ρ,u,v,p)R=(5.444,82.15,-41.05,207725.94)。圖8展示了時間迭代1000步的計算結果?？梢钥闯?，全波的HLLC格式出現了明顯的激波不穩定現象，而LD-Rusanov格式保留了原始Rusanov格式的穩定性，得到了清晰的激波面。

圖8 穩態斜激波問題的密度等值線

4.8 圓柱繞流問題

計算馬赫數為20的無粘圓柱繞流問題來評估不同數值格式在高超聲速流動問題中的激波穩定性。該算例具體的計算區域和初邊值條件的描述參考文獻[12]。本文采用40×80的非規則結構四邊形網格來劃分計算區域。圖9展示了t=4時的密度等值圖?？梢钥闯?，HLLC格式表現出明顯的激波失穩現象，出現了嚴重的紅玉現象，而LD -Rusanov格式和另外兩種穩定版本的HLLC格式均得到了穩定的數值解。

圖9 無粘圓柱繞流問題的密度等值圖

4.9 正激波繞射問題

計算馬赫數為5.09的正激波繞射問題來驗證格式的魯棒性。采用200×200的均勻網格來劃分計算區域[0,1]×[0,1]，其中角點位于(x,y)=(0.05,0.6)處。初始時正激波右側靜止流體的密度為1.4，壓力為1。角點上方左側的入口邊界條件可以由Rankine -Hugoniot關系式計算得到，區域底部和角點下方左側設置為反射邊界條件，在右邊界所有變量的梯度設置為零，區域頂部的邊界條件會隨著激波的運動而調整。圖10展示了t=0.15 時的密度等值圖?？梢钥闯?，HLLC格式在正激波的上方出現了明顯的激波失穩現象，而LD -Rusanov格式和另外兩種穩定版本的HLLC格式均得到了穩定的正激波。值得注意的是，在該算例中本文構造的LD -Rusanov格式表現出了最好的穩定性。

圖10 正激波繞射問題的密度等值圖

5 結 論

基于單波的Rusanov求解器構造了一種低耗散的LD -Rusanov黎曼求解器。使用雙曲正切函數和五階WENO格式來重構界面兩側的密度值。利用邊界變差下降算法從密度重構值的不同組合中挑選出使界面兩側密度差最小的密度值，從而最大限度地減小Rusanov求解器耗散項中的密度差，進而提高格式對于接觸間斷的分辨率。全波的HLLC格式通過修改HLL格式的波系結構來提高接觸間斷的分辨率，因此無法保留HLL格式的激波穩定性。本文構造的LD -Rusanov格式使用一種代數方法來提高接觸間斷的捕捉能力，并沒有修改原Rusanov求解器的波系結構，因此很好地保留了Rusanov的激波穩定性。一系列數值實驗表明LD -Rusanov格式比全波的HLLC格式具有更好的接觸間斷捕捉能力和魯棒性。本文構造的低耗散數值格式展現出了良好的應用前景，因此將其應用于復雜流動問題(如可壓縮湍流和化學反應流)的數值模擬值得未來進一步的研究。