一種魯棒的HLLC格式及其穩(wěn)定性分析

胡立軍， 袁海專

(1.衡陽師范學院 數學與統(tǒng)計學院，衡陽 421002；2.湘潭大學 數學與計算科學學院，湘潭 411105)

1 引 言

近幾十年，基于Godunov方法的近似黎曼解法器在涉及激波和其他間斷的可壓縮流問題中取得了巨大的成功。主要包括Roe[1]，HLL[2]，HLLC[3]，HLLEM[4]和AUSM-類格式[5-7]等。其中，HLL-類格式由于簡單、精確、滿足熵條件和保正性并且易于推廣到不同的雙曲系統(tǒng)而受到青睞。根據能否分辨線性波，可將其分成兩類。一類是能夠分辨線性波的全波HLL-類格式，包括HLLC，HLLEM，HLLE+[8]和HLLI[9]，HLLEMS[10]，HLLC-ADC[11]和HLLC-SWM[12]等,另一類是非全波HLL-類格式，包括HLL，HLLE[13]，HLLS,HLLCM[14]和HLL-CPS[15]等。在計算剪切占優(yōu)的流動現象、混合流體、燃燒面和物質界面時，需要使用全波HLL-類格式來求流通量。然而，在模擬某些超聲速強激波問題時，全波HLL-類格式會出現各種形式的不穩(wěn)定現象，尤其是carbuncle現象[16]。

Quirk[16]首先對數值計算中的激波異?，F象進行系統(tǒng)分類，認為基于任何單一黎曼解法器的Godunov方法都會存在缺陷，并且提出了一種自適應的組合黎曼解法器。Gressier等[17]基于特殊的奇偶失聯擾動分析和數值觀察發(fā)現，格式的嚴格穩(wěn)定性和精確捕捉接觸波這兩種特性是無法兼容的。Liou[18]認為質量通量獨立于壓力項的數值格式不會出現carbuncle現象，并且通過修正線性波的波速來設計一些激波穩(wěn)定的接觸捕捉格式。Dumbser等[19]利用矩陣方法研究穩(wěn)態(tài)激波的穩(wěn)定性，發(fā)現凸激波比凹激波具有更好的穩(wěn)定性，并且上游馬赫數越大穩(wěn)定性越差。Shen等[14]認為格式的剪切粘性不足是誘發(fā)不穩(wěn)定現象的唯一原因。

增加耗散已經成為改善激波異?，F象的主流方法。Kim等[20]提出了一種HLLC-HLL混合方法，借助于開關函數，使得激波層內的橫向通量由不穩(wěn)定的HLLC轉換成穩(wěn)定的HLL，并且使用加權平均通量方法(WAF)來獲得高階精度。Huang等[21]在部分分量上對HLLC和HLL進行加權來消除不穩(wěn)定現象。一些研究人員致力于開發(fā)真正多維的HLLC通量來抑制不穩(wěn)定現象[22,23]。Rodionov[24]將N-S方程的粘性機制引入到無粘歐拉方程來消除Goudnov型數值格式的carbuncle現象。Chen等[25]在激波層亞聲速區(qū)增加剪切波粘性來消除低耗散迎風格式的carbuncle現象。Xie等[26]通過壓力控制方法來限制壓力擾動的傳播，進而治愈HLLC格式的不穩(wěn)定性。Simon等[11,12]使用反擴散控制和選波修正兩種策略來控制HLLC格式的耗散機制，進而抑制其不穩(wěn)定性。

雖然研究人員對激波不穩(wěn)定性開展了大量研究，但是仍然沒有普遍接受的機理分析，這導致了很多改進方法都有一些局限性。本文結合小擾動分析方法和數值實驗來研究HLLC格式的激波不穩(wěn)定性，據此構造一種激波穩(wěn)定的HLLC通量。一系列數值實驗證明了新格式的魯棒性和精度。

2 控制方程組

守恒形式的二維無粘可壓縮歐拉方程組：

?U/?t+?F(U)/?x+?G(U)/?y=0

(1)

式中

(2)

(3)

(4)

(5)

3 HLLC格式

Harten等[2]利用積分關系和 Rankine -Hugoniot 條件提出了一種簡單的兩波近似的HLL格式，其通量表達式為

(6)

左右波速分別為

(7)

(8)

式中

(9)

(10)

式中k=L,R。

4 HLLC格式的激波穩(wěn)定性分析

4.1 小擾動分析

首先，在x-方向(縱向)添加奇偶對稱擾動，擾動后的狀態(tài)為

(11)

超聲速(M0> 1)條件下，擾動的演化方程為

(12)

亞聲速(0

(13)

接下來，考慮在y-方向(橫向)添加如下奇偶對稱擾動

(14)

采用HLLC格式計算y-方向的數值通量Gi,j + 1/2和Gi,j - 1/2，從而得到橫向擾動的演化方程為

(15)

4.2 橫向數值通量與不穩(wěn)定性之間的關系

圖1 超聲速條件下縱向擾動的演化趨勢

圖2 亞聲速條件下縱向擾動的演化趨勢

5 一種魯棒的HLLC格式

基于4.1節(jié)的結論，在計算橫向界面的數值通量時，增加熵波粘性和剪切波粘性來抑制HLLC格式的不穩(wěn)定現象，具體表達式為

(16)

式中EV0和SV0分別為熵波粘性項和剪切波粘性項，Δq=qR-qL。

在熵波粘性項的單獨作用下，橫向擾動的演化方程為

圖3 橫向擾動的演化趨勢

圖4 二維Sedov爆轟波問題的密度等值線

(17)

此時，密度擾動可以有效衰減。在剪切波粘性項的單獨作用下，橫向擾動的演化方程為

(18)

此時，剪切速度的擾動可以有效衰減。因此只有在兩者共同的作用下，橫向所有物理量的擾動才能有效衰減，從而不會誘發(fā)不穩(wěn)定現象。

5.1 開關函數

在計算中，為了保留原HLLC格式精確分辨接觸面和剪切層的優(yōu)點，采用不需要引入經驗參數的界面壓力比來探測激波的位置[10]，

h= min. (pL/pR,pR/pL)

(19)

上述分析表明,不穩(wěn)定現象僅與橫向界面的通量有關。在計算x-方向網格界面(i+1/2,j)的數值通量時，本文探測與其相鄰的y-方向的四個界面是否處在激波層，反之亦然。即

hx= min. (hi,j - 1/2,hi,j + 1/2,hi + 1,j - 1/2,hi + 1,j + 1/2)

hy= min. (hi - 1/2,j,hi + 1/2,j,hi - 1/2,j + 1,hi + 1/2,j + 1)

(20)

利用余弦函數來定義開關函數

f= [1+cos(πh)]/2

(21)

在激波層內，由于存在壓力差，所以0

5.2 R-HLLC格式

Xu等[27]認為激波層超聲速區(qū)的數值結構在擾動下是穩(wěn)定的，但是亞聲速區(qū)的擾動會放大，進而導致激波失穩(wěn)，因此定義亞聲速區(qū)探測函數[25]

(22)

式中M為馬赫數。顯然，在超聲速區(qū)g=0；在亞聲速區(qū)0

考慮式(21,22)，激波面橫向通量上添加的熵波與剪切波粘性項可表示為

(23)

從而，一種魯棒的HLLC(R-HLLC，Robust HLLC)格式的通量可表示為

(24)

由TEV和TSV的定義可知,計算強激波亞聲速區(qū)的橫向數值通量時，R-HLLC格式接近于兩波近似的HLL格式；在其余的大部分地方，R-HLLC格式等價于HLLC格式。由于HLL格式和HLLC格式對于強激波的捕捉能力基本一致，因此局部增加橫向通量的粘性不會導致計算結果的失真，后面數值實驗的結果也說明了這一點。此外，Chen等[25]引入剪切粘性來抑制不穩(wěn)定性。與其相比，本文在三個方面進行了改進,(1) 利用小擾動分析探究了不穩(wěn)定現象的根源；(2) 僅在橫向通量上添加粘性，而文獻[25]在橫向和縱向都添加了粘性；(3) 引入熵波和剪切波粘性來消除不穩(wěn)定現象，而文獻[25]只引入了剪切波粘性。數值實驗的結果表明，單純添加剪切粘性無法完全抑制不穩(wěn)定現象。

6 數值實驗

通過典型算例來檢驗R-HLLC格式的性能，尤其是計算強激波問題時的魯棒性。

6.1 二維穩(wěn)態(tài)激波問題

考慮類似于4.1節(jié)小擾動分析的二維穩(wěn)態(tài)激波問題，馬赫數為7的穩(wěn)態(tài)激波位于x=0.5處，區(qū)域[0,1]×[0,0.5]劃分成100×50的正方形網格。波前和波后的初始分布為

(25)

左右邊界固定為波前和波后狀態(tài)，上下使用周期性邊界條件。從圖5可以看出，原始的HLLC格式和文獻[25]構造的HLLC+SV格式都出現了明顯的不穩(wěn)定現象，而R-HLLC格式有效地抑制了不穩(wěn)定現象的發(fā)生,表明僅增加剪切粘性不能完全消除HLLC格式的不穩(wěn)定性。

6.2 隨機數值噪聲問題

馬赫數為20的正激波沿x-方向傳播，波前和波后的初始狀態(tài)為

(26)

區(qū)域[0,1000]×[0,20]劃分成步長Δx= Δy= 1的正方形網格。在波前流場的初始分布上添加取值從-0.5×10-6到0.5×10-6的隨機擾動。圖6 展示了t=40時的密度等值線?？梢钥闯?，HLLC格式形成了明顯的carbuncle，激波結構完全破壞，而R-HLLC格式展示了清晰的激波面。圖7給出了整個區(qū)域內y-方向速度最大值隨時間的變化情況。HLLC格式計算得到的 max(|v|) 由10-6量級迅速增長到100量級，而R-HLLC格式計算得到的 max(|v|) 一直維持在初始時的10-6量級。

6.3 高超聲速繞柱流問題

計算高超聲速繞柱流問題的穩(wěn)態(tài)解是研究數值格式是否會遭遇致命的carbuncle現象的一個常規(guī)測試問題。馬赫數為20的高超聲速流流經一圓柱體，問題的詳細描述參見文獻[21]。本文采用40×80的結構化貼體四邊形網格。圖8展示了t=4時的計算結果?？梢钥闯?，HLLC格式的激波面在流場中心線附近出現了輕微的不穩(wěn)定現象，而本文構造的R-HLLC格式得到了清晰的穩(wěn)態(tài)弓形激波。

圖5 二維穩(wěn)態(tài)激波問題的密度等值線

圖6 隨機數值噪聲問題的密度等值線

6.4 激波后臺階繞射問題

馬赫數為5.09的右行超聲速激波繞過一個90°的角，區(qū)域[0,1]×[0,1]劃分成步長Δx= Δy=1/400的正方形網格，角點位于(0.05,0.625)處。初始條件和邊界條件參見文獻[11]。圖9展示了t= 0.1561時的計算結果。可以看出，HLLC格式出現了明顯的激波失穩(wěn)現象，而R-HLLC格式則完全消除了偽振蕩，得到了清晰的激波面。對于次激波、接觸面、膨脹波和入射激波兩者具有相同的分辨率。

6.5 Rayleigh-Taylor不穩(wěn)定性問題

初始時，區(qū)域[0,0.25]×[0,1]上下部分兩種流體的初始狀態(tài)為

(27)

圖7 y-方向速度的最大取值

圖8 高超聲速繞柱流問題的密度等值線

圖9 激波后臺階繞射問題的密度等值線

圖10 Rayleigh-Taylor不穩(wěn)定性問題的密度等值線

6.6 激波/邊界層相互作用問題

計算復雜的激波/邊界層相互作用問題來檢驗格式的精度。計算區(qū)域[0,1]×[0,0.5]劃分成步長Δx= Δy= 1/512的均勻網格。初始條件為

(28)

無粘通量的計算采用二階MUSCL重構，粘性通量的計算采用二階中心差分，時間離散采用二階 Runge -Kutta 格式。計算中，雷諾數取220。從 圖11 可以看出，R-HLLC和HLLC格式得到了幾乎相同的解。圖12表明，兩種格式的計算結果沿著底部壁面的密度分布完全吻合。該算例證明了R-HLLC格式在計算粘性流問題時的有效性。

圖11 激波/邊界層相互作用問題的密度等值線

圖12 底部壁面的密度分布

7 結 論

構造了一種激波穩(wěn)定的HLLC格式。穩(wěn)定性分析表明,HLLC格式橫向通量的密度與剪切速度的擾動不會發(fā)生衰減。在橫向數值通量上增加熵波與剪切波粘性來抑制激波不穩(wěn)定現象的發(fā)生。為了保留HLLC格式低耗散性的優(yōu)點，定義開關函數來自動控制添加粘性的位置和大小。新的格式實施簡單，僅需在現有HLLC程序的基礎上增加若干行代碼，并且不會引入依賴于問題的可調參數。數值算例的計算結果展示了新的HLLC格式具有高分辨率和強魯棒性的優(yōu)點。此外，構造其他魯棒的Godunov型激波捕捉格式、開發(fā)三維和非結構網格上的R-HLLC格式，并將其應用到其他的高超聲速復雜流動問題的數值模擬值得進一步研究。