胡立軍， 袁海專

(1.衡陽(yáng)師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，衡陽(yáng) 421002；2.湘潭大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，湘潭 411105)

1 引 言

近幾十年，基于Godunov方法的近似黎曼解法器在涉及激波和其他間斷的可壓縮流問(wèn)題中取得了巨大的成功。主要包括Roe[1]，HLL[2]，HLLC[3]，HLLEM[4]和AUSM-類格式[5-7]等。其中，HLL-類格式由于簡(jiǎn)單、精確、滿足熵條件和保正性并且易于推廣到不同的雙曲系統(tǒng)而受到青睞。根據(jù)能否分辨線性波，可將其分成兩類。一類是能夠分辨線性波的全波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]等。在計(jì)算剪切占優(yōu)的流動(dòng)現(xiàn)象、混合流體、燃燒面和物質(zhì)界面時(shí)，需要使用全波HLL-類格式來(lái)求流通量。然而，在模擬某些超聲速?gòu)?qiáng)激波問(wèn)題時(shí)，全波HLL-類格式會(huì)出現(xiàn)各種形式的不穩(wěn)定現(xiàn)象，尤其是carbuncle現(xiàn)象[16]。

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

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

雖然研究人員對(duì)激波不穩(wěn)定性開(kāi)展了大量研究，但是仍然沒(méi)有普遍接受的機(jī)理分析，這導(dǎo)致了很多改進(jìn)方法都有一些局限性。本文結(jié)合小擾動(dòng)分析方法和數(shù)值實(shí)驗(yàn)來(lái)研究HLLC格式的激波不穩(wěn)定性，據(jù)此構(gòu)造一種激波穩(wěn)定的HLLC通量。一系列數(shù)值實(shí)驗(yàn)證明了新格式的魯棒性和精度。

2 控制方程組

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

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

(1)

式中

(2)

(3)

(4)

(5)

3 HLLC格式

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

(6)

左右波速分別為

(7)

(8)

式中

(9)

(10)

式中k=L,R。

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

4.1 小擾動(dòng)分析

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

(11)

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

(12)

亞聲速(0

(13)

接下來(lái)，考慮在y-方向(橫向)添加如下奇偶對(duì)稱擾動(dòng)

(14)

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

(15)

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

圖1 超聲速條件下縱向擾動(dòng)的演化趨勢(shì)

圖2 亞聲速條件下縱向擾動(dòng)的演化趨勢(shì)

5 一種魯棒的HLLC格式

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

(16)

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

在熵波粘性項(xiàng)的單獨(dú)作用下，橫向擾動(dòng)的演化方程為

圖3 橫向擾動(dòng)的演化趨勢(shì)

圖4 二維Sedov爆轟波問(wèn)題的密度等值線

(17)

此時(shí)，密度擾動(dòng)可以有效衰減。在剪切波粘性項(xiàng)的單獨(dú)作用下，橫向擾動(dòng)的演化方程為

(18)

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

5.1 開(kāi)關(guān)函數(shù)

在計(jì)算中，為了保留原HLLC格式精確分辨接觸面和剪切層的優(yōu)點(diǎn)，采用不需要引入經(jīng)驗(yàn)參數(shù)的界面壓力比來(lái)探測(cè)激波的位置[10]，

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

(19)

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

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)

利用余弦函數(shù)來(lái)定義開(kāi)關(guān)函數(shù)

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

(21)

在激波層內(nèi)，由于存在壓力差，所以0

5.2 R-HLLC格式

Xu等[27]認(rèn)為激波層超聲速區(qū)的數(shù)值結(jié)構(gòu)在擾動(dòng)下是穩(wěn)定的，但是亞聲速區(qū)的擾動(dòng)會(huì)放大，進(jìn)而導(dǎo)致激波失穩(wěn)，因此定義亞聲速區(qū)探測(cè)函數(shù)[25]

(22)

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

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

(23)

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

(24)

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

6 數(shù)值實(shí)驗(yàn)

通過(guò)典型算例來(lái)檢驗(yàn)R-HLLC格式的性能，尤其是計(jì)算強(qiáng)激波問(wèn)題時(shí)的魯棒性。

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

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

(25)

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

6.2 隨機(jī)數(shù)值噪聲問(wèn)題

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

(26)

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

6.3 高超聲速繞柱流問(wèn)題

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

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

圖6 隨機(jī)數(shù)值噪聲問(wèn)題的密度等值線

6.4 激波后臺(tái)階繞射問(wèn)題

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

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

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

(27)

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

圖8 高超聲速繞柱流問(wèn)題的密度等值線

圖9 激波后臺(tái)階繞射問(wèn)題的密度等值線

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

6.6 激波/邊界層相互作用問(wèn)題

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

(28)

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

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

圖12 底部壁面的密度分布

7 結(jié) 論

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