胡立軍，翟 健，袁 禮

(1.衡陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，衡陽 421002； 2.曙光信息產(chǎn)業(yè)(北京)有限公司，北京 100193； 3.中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院 計算數(shù)學(xué)與科學(xué)工程計算研究所，北京 100190)

1 引 言

近幾十年，求解可壓縮流的Godunov類型黎曼解法器的研究吸引了大量計算流體力學(xué)研究人員的關(guān)注[1]。中心格式[2]、通量差分裂格式[3]、通量向量分裂格式[4,5]、HLL類格式[6-12]和AUSM-類格式[13]已經(jīng)廣泛地應(yīng)用于流體力學(xué)的數(shù)值計算中。這些數(shù)值方法都是基于方向分裂的思路來設(shè)計的，在計算網(wǎng)格界面數(shù)值通量時，只考慮了法向的波系，忽略了橫向波系的影響[14]。在求解多維問題時，不僅會降低格式的分辨率，而且穩(wěn)定性CFL數(shù)也會減小，從而影響格式的計算效率[15]。許多研究人員做了大量的研究來克服這些問題，其基本思路是在數(shù)值通量的計算中增加橫向波的影響，進(jìn)而設(shè)計一種真正多維的數(shù)值格式。

一些研究人員嘗試對一維黎曼解法器進(jìn)行復(fù)雜的組合來實現(xiàn)多維迎風(fēng)性。如Rumsey等[16]將多維性引入到一維Roe黎曼解法器中，Colella[17]提出的多維角點輸運格式和LeVeque[18]設(shè)計的多維波傳播算法。盡管這樣處理可以增加穩(wěn)定性CFL數(shù)，但是常會以求解大量的一維黎曼問題作為代價。另外一些研究人員嘗試去建立更加復(fù)雜的多維波傳播模型，從而獲得一種線性化的多維Roe格式[15,19]。但是該解法器只能用來計算歐拉方程，很難推廣到其他的守恒律系統(tǒng)。此外，Wendroff[20]利用HLLE格式適用于任何守恒律系統(tǒng)這一特性來構(gòu)造一種多維的HLLE格式，但是該格式引入了九個常數(shù)狀態(tài)，形式非常復(fù)雜。

Balsara[21]基于一維HLLE格式構(gòu)造了一種簡單高效的真正二維黎曼解法器。與其他多維格式不同，該格式具有簡單的封閉形式，其算法很容易實現(xiàn)。為了克服二維HLLE格式無法捕捉接觸間斷和剪切波的缺陷，引入一個子波系，在HLLC格式的基礎(chǔ)上對二維HLLE格式進(jìn)行了改進(jìn)。文獻(xiàn)[22]的結(jié)果表明，真正二維的HLLC格式能夠精確分辨接觸間斷。Balsara等[23]對具有自相似內(nèi)部結(jié)構(gòu)的多維黎曼問題提出了一種求解策略，并且構(gòu)造了一種可以用來求解守恒律和非守恒律問題的多維HLLI格式，該格式在數(shù)值計算中不需要詳細(xì)區(qū)分接觸間斷的方向。

采用Balsara的思路，Mandal等[24]利用HLL-CPS格式的優(yōu)點并且結(jié)合Zha-Bilgen分裂方法，構(gòu)造了一種真正二維的HLL格式。胡立軍等[14]基于 Toro -Vázquez 分裂方法，在二維結(jié)構(gòu)網(wǎng)格上構(gòu)造了一種真正二維的HLL黎曼解法器。Qu等[25]將該黎曼解法器推廣到了二維曲線坐標(biāo)系。結(jié)合通量分裂方法設(shè)計的多維黎曼解法器，形式簡單并且計算代價小，為構(gòu)造真正多維黎曼解法器提供了一條新的思路。

本文以傳統(tǒng)的一維TV格式[26]為基礎(chǔ)，構(gòu)造一種簡單的真正二維TV通量分裂格式。

2 Toro -Vázquez 分裂方法

2.1 控制方程組

直角坐標(biāo)系下的二維歐拉方程組為

(1)

式中

(2)

ρ為密度，u和v為x方向和y方向的流體速度，p為壓力，E為總能。狀態(tài)方程為

(3)

本文比熱比γ=1.4。

用守恒型數(shù)值方法求解式(1)，

(4)

式中Fi + 1/2,j和Gi,j + 1/2為x和y方向的數(shù)值通量。

2.2 對流/壓力通量分裂

Toro等[26]將歐拉方程的能量方程分裂成動能項和壓力項，從而將通量分裂成對流通量和壓力通量。以x方向為例：

F(U)=Ax(U)+Px(U)=

(5)

式中 內(nèi)能e=p/[(γ-1)ρ]。

對流通量Ax(U)的Jacobian矩陣為

(6)

矩陣B的特征值為

(7)

壓力通量Px(U)的Jacobian矩陣為

(8)

矩陣M的特征值為

(9)

3 二維TV通量分裂格式

3.1 界面中點一維通量的計算

根據(jù)對流/壓力通量分裂式(5)，中點處的一維通量可以分裂成兩部分之和。

(10)

3.1.1 中點對流數(shù)值通量

使用TV格式[26]來計算網(wǎng)格界面中點處的對流數(shù)值通量。以x方向為例

(11)

式中 界面速度

(12)

3.1.2 中點壓力數(shù)值通量

同樣使用TV格式[26]來計算中點處的壓力數(shù)值通量。將狀態(tài)方程p=(γ-1)ρe代入壓力通量表達(dá)式(5)，可以得到壓力通量的計算公式為

(13)

界面速度u*的定義見式(12)，界面壓力p*的計算公式為

(15)

3.2 界面角點二維通量的計算

根據(jù)分裂式(5)，角點處的二維通量可以分裂成兩部分之和，即

(16)

3.2.1 角點對流數(shù)值通量

使用迎風(fēng)格式來求解角點處的對流數(shù)值通量。通過考慮角點處橫向波對數(shù)值通量的影響，來實現(xiàn)格式真正二維的特性。其具體的計算公式為

(17)

圖1 結(jié)構(gòu)化網(wǎng)格和角點處黎曼問題的四個初態(tài)

Fig.1 Structured grid and four initial states at the corner

(18)

迎風(fēng)方向k1和k2按如下方式來選取。

(19)

3.2.2 角點壓力數(shù)值通量

為了更好地保留一維TV格式的優(yōu)點，與文獻(xiàn)[14,25]中基于 Toro -Vázquez 分裂方法構(gòu)造二維HLL格式不同的是，本文壓力數(shù)值通量使用TV格式而不是HLL格式來求解：

(20)

式中 右端第一項為考慮了影響域面積加權(quán)的TV通量，其計算公式為

(21)

(22)

式(20)右端第二項為橫向壓力通量對角點數(shù)值通量的貢獻(xiàn)。其中，Py為y方向的壓力通量，定義為

(23)

如圖1(b)所示，波速SU,SD,SL和SR分別代表沿著上下左右四個方向傳播的最快波速，本文按照文獻(xiàn)[14]的定義來計算。

3.3 界面數(shù)值通量

(24)

4 數(shù)值結(jié)果

使用本文構(gòu)造的二維TV通量分裂格式求解幾個典型的數(shù)值算例來驗證格式的精度、魯棒性和計算效率。

4.1 二維黎曼問題

計算文獻(xiàn)[27]的二維黎曼問題來比較三種數(shù)值格式，即本文構(gòu)造的GT-TV格式、一維TV格式[26]和真正二維的GM-TV-HLL格式[14]。計算區(qū)域[0,1]×[0,1]劃分成400×400的矩形網(wǎng)格，初始條件為

(25)

該算例沒有精確解，本文采用HLLC格式在800×800的密網(wǎng)格下所得到的計算結(jié)果作為參考解，三種格式的數(shù)值解均采用400×400的粗網(wǎng)格來計算。圖2展示了計算時間T=0.25時，參考解和三種數(shù)值解的密度輪廓。為了更加精細(xì)地比較三種數(shù)值格式之間的差異，計算數(shù)值解與參考解之間的誤差，其表達(dá)式為

(26)

在計算數(shù)值通量時，真正二維的數(shù)值格式不僅考慮界面法向波系的影響，還考慮了界面橫向波系的影響，因此比傳統(tǒng)的一維數(shù)值格式具有更高的分辨率[15,17]。而GM-TV-HLL格式的壓力通量采用耗散較大的HLL格式來計算，從而會帶來更大的數(shù)值誤差。本節(jié)數(shù)值試驗的結(jié)果也證明了這 一點。

圖2 二維黎曼問題的計算結(jié)果

Fig.2 Solutions of two -dimensional Riemann problem

4.2 隨機擾動問題

進(jìn)行數(shù)值模擬時，由于各種原因，在初始分布和計算過程中會引入很小的隨機擾動。在計算中，隨機擾動的發(fā)展過程，常用來檢驗數(shù)值格式的魯棒性。馬赫數(shù)為10的平面激波從左向右移動，初始物理量上增加10-10倍隨機擾動。激波右側(cè)具體的初始分布為

(27)

式中αk(k=1,2,3,4)是(0,1)之間的隨機數(shù)，計算中由相應(yīng)的庫函數(shù)生成。計算區(qū)域[0,35]×[0,1]劃分成700×20的矩形網(wǎng)格。

圖3展示了在計算時間T=3時，使用TV格式和GT-TV格式求解隨機擾動問題所得到的計算結(jié)果?？梢钥闯?，真正二維的GT-TV格式不僅能夠精確地捕捉激波的位置(x=30)，而且消除了一維TV格式的波后不穩(wěn)定現(xiàn)象。在計算中，一維TV格式的穩(wěn)定性CFL數(shù)不能超過0.4，而GT-TV格式的穩(wěn)定性CFL數(shù)最大可以取到0.95。

表1 不同數(shù)值格式的誤差比較

Tab.1 Comparison of errors between different numerical schemes

SchemeError‖L2(ρ)‖TV4.8325e-5GM-TV-HLL6.3139e-5GT-TV4.0652e-5

圖3 隨機擾動問題的計算結(jié)果

Fig.3 Solutions of random perturbation problem

4.3 后臺階激波衍射問題

后臺階激波衍射問題對于數(shù)值格式的魯棒性要求很高，許多數(shù)值格式在計算該算例時都會失效。因此，文獻(xiàn)多使用該算例來測試格式的魯棒性。馬赫數(shù)為5.09的右行超聲速氣流繞過90°的角，計算區(qū)域[0,1]×[0,1]劃分成400×400的矩形網(wǎng)格，角點位置(x,y)=(0.05,0.625)。激波右側(cè)的初始條件為

(28)

圖4展示了計算時間T=0.15時，一維TV格式和二維GT-TV格式的計算結(jié)果。二維GT-TV格式消除了一維TV格式在區(qū)域上邊界正激波處的不穩(wěn)定現(xiàn)象。在計算中，一維 TV格式的穩(wěn)定性CFL數(shù)不能超過0.4，二維GT-TV格式的穩(wěn)定性CFL數(shù)最大可以取到0.9。

文獻(xiàn)[11,12]指出，能夠精確分辨接觸間斷的數(shù)值格式都會遭受不同程度的激波不穩(wěn)定現(xiàn)象，并且橫向擾動的增長未能得到有效的抑制是造成激波不穩(wěn)定的根本原因。真正二維的GT-TV格式的數(shù)值通量考慮了橫向波系的影響，因此橫向的數(shù)值粘性可以有效地抑制橫向擾動的增長，從而消除一維TV格式的激波不穩(wěn)定性。4.2節(jié)和4.3節(jié)的數(shù)值試驗也證明了這一結(jié)論。

圖4 激波衍射問題的計算結(jié)果

Fig.4 Solutions of shock diffraction problem

4.4 計算效率的比較

比較一維TV格式和二維GT-TV格式的計算效率。表2展示了用兩種格式計算本文三個數(shù)值算例的CPU時間和CFL數(shù)取格式允許的最大值。如計算二維黎曼問題時，TV格式的CFL數(shù)取0.5，GT-TV格式取0.95。在計算中，GT-TV格式不僅要計算界面中點的數(shù)值通量，還需要計算界面角點的數(shù)值通量，因此單個時間步的計算量會大于TV格式。但是GT-TV格式可以提高穩(wěn)定性CFL數(shù)[17]，從而增加時間步長，減少CPU計算時間。從表2可以看到，相比于一維TV格式，GT-TV格式的計算效率分別提高了 25.20%，38.62%和36.88%。因此，GT-TV格式具有更高的計算效率。

表2 數(shù)值格式計算效率的比較

Tab.2 Comparisons of computational efficiency for numerical schemes

Example1Example2Example3SchemeTime/CFLTime/CFLTime/CFLTV86.58/0.552.17/0.492.64/0.4GT-TV64.76/0.9532.02/0.9558.47/0.9

5 結(jié) 論

一維TV通量分裂格式不僅形式簡單,而且可以精確捕捉一維激波和膨脹激波。但是,在求解界面數(shù)值通量時，只考慮了界面法向的波系，忽視了橫向波系的影響。因此,使用一維TV格式求解二維問題時，不僅會減小穩(wěn)定性CFL數(shù)，而且無法保證格式的魯棒性。本文基于一維TV格式構(gòu)造的二維TV通量分裂格式,通過求解考慮了橫向波系影響的角點數(shù)值通量,來克服一維TV格式的缺點。數(shù)值結(jié)果表明，與一維TV格式相比，本文構(gòu)造的二維TV格式不僅有更高的分辨率和計算效率，而且有更好的魯棒性。構(gòu)造非結(jié)構(gòu)網(wǎng)格上的二維TV通量分裂格式并且將其應(yīng)用到其他的守恒律系統(tǒng),可以作為未來的研究工作。