劉旭亮，范召林，張樹(shù)海1，，李 虎1，，羅 勇1，，孫曉峰

(1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽(yáng) 621000；3.北京航空航天大學(xué) 能源與動(dòng)力工程學(xué)院，北京 100191)

0 引 言

雙曲守恒律方程的空間離散需要構(gòu)造數(shù)值通量，而數(shù)值通量的構(gòu)造方法主要分為兩大類(lèi)：一類(lèi)方法是通量分裂，包括Lax-Friedrichs分裂[1]、Steger-Warming分裂[2]和Van Leer分裂[3]等；另一類(lèi)方法是Riemann求解器，或稱(chēng)為Riemann算子。

Riemann求解器是數(shù)值求解雙曲系統(tǒng)的重要組成部分。對(duì)于守恒形式的雙曲系統(tǒng)，已有許多學(xué)者提出了多種著名的Riemann求解器。精確的Riemann求解器由Godunov[4]在1959年提出，具有耗散小、精度高等優(yōu)點(diǎn)，但計(jì)算量過(guò)大。后續(xù)學(xué)者發(fā)展出了多種近似Riemann求解器并對(duì)其進(jìn)行改良。1981年，Roe[5]提出了著名的近似Riemann求解器，該求解器是對(duì)非線(xiàn)性Euler方程組的特殊線(xiàn)化。但是原始Roe算子的主要缺陷是在特定問(wèn)題中違反熵增條件，因此Roe算子必須進(jìn)行熵修正才能保證計(jì)算的準(zhǔn)確性。1983年，Harten等[6]提出了HLL型Riemann求解器，HLL算子非常高效且有很好的魯棒性，滿(mǎn)足熵增條件，并保持正定性，但HLL算子耗散較大且不能完全解決接觸間斷的問(wèn)題。在HLL算子的基礎(chǔ)上，Einfeldt等[7-8]提出了HLLEM型近似Riemann求解器，這類(lèi)求解器具有正定保持性質(zhì)并且是熵增的。Einfeldt等[8]通過(guò)構(gòu)造反擴(kuò)散系數(shù)來(lái)降低HLL算子的耗散，并證明了HLLEM求解器和Roe求解器的區(qū)別僅在于數(shù)值信號(hào)速度的不同，因此，HLLEM求解器依然可能在計(jì)算激波時(shí)出現(xiàn)不穩(wěn)定。為了克服HLL算子無(wú)法模擬接觸間斷的缺陷，Toro等[9]提出了HLLC型Riemann求解器，恢復(fù)了HLL算子中丟失的接觸波和剪切波。HLLC算子保留了HLL算子的熵增特性和正定保持性質(zhì)，但在計(jì)算激波中仍然會(huì)出現(xiàn)非物理的數(shù)值振蕩。基于由相空間中的路徑積分來(lái)得到數(shù)值耗散項(xiàng)的構(gòu)造思路， Osher等[10]提出了Osher型Riemann求解器，其優(yōu)點(diǎn)是數(shù)值通量光滑并且滿(mǎn)足熵增特性。在具體計(jì)算時(shí)，HLLEM求解器與Roe求解器僅需求解特征值和相關(guān)的左右特征向量，而Osher型Riemann求解器需要計(jì)算整套的特征系統(tǒng)，因此Osher求解器計(jì)算量較大。為了將通量中的線(xiàn)性項(xiàng)和非線(xiàn)性項(xiàng)加以區(qū)別，Liou[11]等提出了一種結(jié)合通量分裂和Riemann求解器的AUSM系列格式，這種方法把壓力項(xiàng)和對(duì)流項(xiàng)進(jìn)行分裂，對(duì)部分問(wèn)題能夠消除計(jì)算激波的不穩(wěn)定性現(xiàn)象。

通過(guò)與激波捕捉格式結(jié)合，近似Riemann求解器廣泛應(yīng)用于超聲速流動(dòng)的數(shù)值模擬。但是，在高維計(jì)算中，使用近似Riemann求解器可能會(huì)遇到激波計(jì)算的不穩(wěn)定性問(wèn)題。1988年，Peery和Imlay[12]首先報(bào)道了計(jì)算激波不穩(wěn)定的紅玉現(xiàn)象（carbuncle phenomenon），他們使用Roe的近似Riemann求解器計(jì)算了鈍頭體周?chē)某曀倭鲌?chǎng)，發(fā)現(xiàn)Riemann求解器的數(shù)值不穩(wěn)定性可能會(huì)嚴(yán)重影響激波計(jì)算的準(zhǔn)確性。從那時(shí)起，許多學(xué)者開(kāi)始研究和解決紅玉現(xiàn)象問(wèn)題。防止Riemann算子計(jì)算激波時(shí)出現(xiàn)不穩(wěn)定性，主要有三類(lèi)方法：第一類(lèi)是Harten熵修正方法[13]，該方法是對(duì)Riemann算子的特征值在接近零時(shí)強(qiáng)制增大。熵修正方法引入了自由參數(shù)，在實(shí)際應(yīng)用中自由參數(shù)設(shè)置過(guò)大或者過(guò)小時(shí)都可能導(dǎo)致計(jì)算崩潰，并且不同的問(wèn)題要選取不同的自由參數(shù)，因此熵修正方法的經(jīng)驗(yàn)性非常強(qiáng)。第二類(lèi)是Quirk提出的混合方法[14]。Quirk注意到，某些耗散小的Riemann求解器總是出現(xiàn)紅玉現(xiàn)象，而其他耗散較大的Riemann求解器則沒(méi)有這種不穩(wěn)定性。Quirk建議在激波區(qū)域使用耗散大的Riemann求解器，而在其他區(qū)域使用耗散較小的Riemann求解器，基于此類(lèi)混合方法來(lái)計(jì)算激波問(wèn)題可以克服數(shù)值不穩(wěn)定現(xiàn)象。盡管具體的組合方式可能有著顯著差異[15-17]，但類(lèi)似的思路已經(jīng)被廣泛采用。第三類(lèi)是構(gòu)造高維Riemann求解器方法。由于在高維問(wèn)題中計(jì)算激波更容易出現(xiàn)不穩(wěn)定性問(wèn)題，一些研究者發(fā)展了高維近似Riemann求解器[18-19]或旋轉(zhuǎn)Riemann求解器[20-21]。此類(lèi)方法本質(zhì)上是高維的，可以部分抑制激波不穩(wěn)定性。根據(jù)Huang等[16]的研究，在實(shí)際計(jì)算中使用混合方法的計(jì)算效率高于旋轉(zhuǎn)Riemann求解器，他們認(rèn)為混合方法更高效，同時(shí)也具有較好的魯棒性。目前在求解高維雙曲守恒律方程時(shí)，更常用的策略是采用局部一維Riemann求解器來(lái)計(jì)算數(shù)值通量，高維Riemann求解器方法應(yīng)用還不夠廣泛。

混合方法能夠解決激波計(jì)算不穩(wěn)定性問(wèn)題，同時(shí)也能避免熵修正方法中自由參數(shù)的經(jīng)驗(yàn)性，因此混合方法在實(shí)際計(jì)算中比較常用，但必須選擇適當(dāng)?shù)幕A(chǔ)Riemann求解器和混合因子。近年來(lái)，Dumbser等[22]和Xie等[23]基于HLLEM算子的反擴(kuò)散矩陣來(lái)構(gòu)造混合Riemann求解器，他們認(rèn)為以HLLEM算子作為基礎(chǔ)，對(duì)接觸波和剪切波分量進(jìn)行修正的混合方法比較合理。

為了消除近似Riemann求解器在數(shù)值模擬激波時(shí)出現(xiàn)的計(jì)算不穩(wěn)定性現(xiàn)象，本文采用混合方法對(duì)反擴(kuò)散矩陣進(jìn)行修正，發(fā)展了一類(lèi)新型具有自適應(yīng)反擴(kuò)散的近似Riemann求解器。該求解器應(yīng)用到高階格式時(shí)，能夠保持差分格式的高階精度，并且計(jì)算穩(wěn)定性較好。

1 數(shù)值格式的構(gòu)造

1.1 控制方程和差分格式

本文的數(shù)值方法主要應(yīng)用于雙曲守恒率方程。以二維可壓縮Euler方程為例：

其中，守恒變量為：

x、y方向的通量分別為：

本文的數(shù)值方法主要是離散Euler方程的對(duì)流項(xiàng)，下面以 ?F/?x的離散為例進(jìn)行說(shuō)明。

選取圖1的網(wǎng)格模板，網(wǎng)格結(jié)點(diǎn)處的一階導(dǎo)數(shù)由線(xiàn)性中心緊致格式[24]來(lái)得到。

圖1 緊致格式的模板Fig.1 Stencil for compact scheme

為了能夠計(jì)算包含激波等復(fù)雜流動(dòng)的非線(xiàn)性問(wèn)題，緊致格式（4）中可以采用Riemann求解器來(lái)近似。通過(guò)半結(jié)點(diǎn)處的通量把迎風(fēng)和耗散性引入到差分格式，下標(biāo)L表示半結(jié)點(diǎn)左側(cè)方向，R表示半結(jié)點(diǎn)右側(cè)方向。其中半結(jié)點(diǎn)處的值和由混合加權(quán)非線(xiàn)性插值[25]得到。因?yàn)榈南禂?shù)關(guān)于與是對(duì)稱(chēng)的，所以只需給出的構(gòu)造方法。

其中非線(xiàn)性權(quán)為：

其中線(xiàn)性權(quán)為：

光滑因子的詳細(xì)公式參考文獻(xiàn)[25]。優(yōu)化參數(shù)的選擇為：α=0.25， σ =0.67。

本文記這類(lèi)格式為混合優(yōu)化非線(xiàn)性緊致格式（hybrid optimized nonlinear compact scheme, HONCS），本文主要采用五階混合優(yōu)化非線(xiàn)性緊致格式，簡(jiǎn)單記為HONCS5。

1.2 數(shù)值信號(hào)速度

由于數(shù)值信號(hào)速度在Riemann求解器的計(jì)算中起著影響耗散性的重要作用，同時(shí)構(gòu)造所有的Riemann求解器都要先指定數(shù)值信號(hào)速度。因此在說(shuō)明Riemann求解器之前，先給出幾種常用的數(shù)值信號(hào)速度。

1.2.1 Roe類(lèi)型[5]

1.2.2 Einfeldt 1988類(lèi)型[7]

1.2.3 Einfeldt 1991類(lèi)型[8]

1.2.4 Batten類(lèi)型[26]

1.2.5 Davis類(lèi)型[27]

1.2.6 Toro類(lèi)型[28]

其中K=L或者K=R。

數(shù)值信號(hào)速度為：SL=qL-ηLcL，SR=qR+ηRcR。

1.3 傳統(tǒng)Riemann求解器

為了行文方便，選取如下符號(hào)：S-=min(SL,0)和S+=max(SR,0)。

定義左右狀態(tài)通量之差為：

則有如下關(guān)系式：

1.3.1 Roe 型Riemann求解器[5]

特征值為：

原始Roe算子可以表示成如下幾種形式：

直接采用原始Roe算子求解激波問(wèn)題，通常會(huì)出現(xiàn)熵違反解[28]，必須對(duì)特征值進(jìn)行熵修正才能得到正確解。但是一般來(lái)說(shuō)，熵修正的經(jīng)驗(yàn)性很強(qiáng)，并且特定問(wèn)題必須采用特定的熵修正值，因此目前還沒(méi)有標(biāo)準(zhǔn)的做法。

1.3.2 HLL型Riemann求解器[6]

HLL通量用S-、S+可以表達(dá)為：

把關(guān)系式（9）代入式（15），則可簡(jiǎn)潔地表達(dá)為：

1.3.3 HLLC型Riemann求解器[9,28]

中間星區(qū)守恒變量為：

其中：

1.3.4 HLLEM型Riemann求解器[8]

以二維問(wèn)題為例，HLLEM通量[8]為：

把關(guān)系式（9）代入式（21），則可表示為：

其中：

反擴(kuò)散系數(shù)的定義為：

1.3.5 通量分裂型Riemann求解器

一些學(xué)者[29-31]采用通量分裂的方式來(lái)構(gòu)造Riemann算子，即通過(guò)來(lái)計(jì)算數(shù)值通量。但是Steger-Warming分裂[2]、Van Leer分裂[3]等是逐點(diǎn)分裂的，無(wú)法表示成Riemann算子的通用形式。只有局部Lax-Friedrichs分裂可以變換成通用形式其中局部分裂系數(shù)為表示局部分裂點(diǎn)坐標(biāo)。通量分裂型Riemann算子不是本文關(guān)注的重點(diǎn)，因此不做詳細(xì)比較和討論。

1.4 新型自適應(yīng)耗散RAD型Riemann求解器

由于HLLEM算子和原始Roe算子在數(shù)值信號(hào)速度相同時(shí)完全等價(jià)[8]，所以HLLEM算子中的反擴(kuò)散系數(shù)在激波附近取值較大，導(dǎo)致在計(jì)算激波時(shí)出現(xiàn)不穩(wěn)定現(xiàn)象。

為了能夠保持HLLEM算子在光滑區(qū)的低耗散特性，同時(shí)能夠穩(wěn)定地計(jì)算激波，本文構(gòu)造了新型混合反擴(kuò)散矩陣，以二維問(wèn)題為例：

其中， φi∈ [0,1]是自適應(yīng)混合因子：

流場(chǎng)光滑量度為：

本文中取無(wú)量綱參數(shù) ε =1×10-40。本文中若無(wú)特殊說(shuō)明，則無(wú)量綱參數(shù) δ =1×10-4。

顯然，矩陣D的特征值可以寫(xiě)成：

當(dāng)S+=SR＞0且S-=SL＜0時(shí)，即流動(dòng)是亞聲速時(shí)，令則可以得到反擴(kuò)散矩陣

通過(guò)對(duì)比RAD求解器與Roe、HLL和HLLEM等Riemann求解器的公式，可以得到如下的等價(jià)關(guān)系：

因此，傳統(tǒng)的Roe、HLL、HLLEM等Riemann求解器都是本文新型RAD求解器的特殊形式。

根據(jù)自適應(yīng)混合因子的定義，當(dāng)流場(chǎng)處于光滑區(qū)時(shí)， φi≈ 1，即當(dāng)流場(chǎng)處于間斷區(qū)時(shí)，所以，新型RAD求解器在光滑區(qū)耗散較小，而在間斷區(qū)能夠抑制激波計(jì)算的不穩(wěn)定性。

2 格式精度驗(yàn)證和頻譜分析

2.1 數(shù)值信號(hào)速度的選擇

為了確定在Riemann算子中的數(shù)值信號(hào)速度，本文采用修正Sod激波管問(wèn)題[26,28]來(lái)驗(yàn)證和比較計(jì)算結(jié)果。該問(wèn)題包含激波、膨脹波、接觸間斷等，能夠很好的評(píng)估數(shù)值方法的熵滿(mǎn)足特性。計(jì)算網(wǎng)格采用100個(gè)點(diǎn)，計(jì)算到無(wú)量綱時(shí)間t= 0.2，間斷左右兩邊參數(shù)為：

圖2給出了基于一階迎風(fēng)格式和HLLEM算子的不同數(shù)值信號(hào)速度的計(jì)算結(jié)果，可以看出，Roe類(lèi)型、Einfeldt 1988類(lèi)型和Einfeldt 1991類(lèi)型的信號(hào)速度的計(jì)算結(jié)果明顯有振蕩。Toro類(lèi)型的信號(hào)速度比較復(fù)雜，不易從數(shù)學(xué)上證明正定保持特性[26]。Davis類(lèi)型信號(hào)速度耗散比Batten類(lèi)型的大。因此，本文全部采用Batten類(lèi)型數(shù)值信號(hào)速度。

圖2 數(shù)值信號(hào)速度的比較Fig.2 Comparison of numerical signal velocities

2.2 精度驗(yàn)證

為了驗(yàn)證數(shù)值方法的精度，選取如下的格式：時(shí)間離散采用三階TVD Runge-Kutta格式[32]，空間離散采用五階HONCS差分格式，數(shù)值通量計(jì)算采用RAD型Riemann求解器。

以一維Euler方程的對(duì)流密度波問(wèn)題為驗(yàn)證標(biāo)準(zhǔn)[33]，其精確解為： ( ρ,u,p)=(1+0.1sin[π(x-t)],1,1)。計(jì)算域取 [0 ,2]，邊界取周期性邊界條件，計(jì)算的最終時(shí)刻為t=20，初始計(jì)算網(wǎng)格為15個(gè)點(diǎn)，CFL= 0.5。隨著網(wǎng)格點(diǎn)數(shù)增多，CFL數(shù)為CFL2=CFL1(N2/N1)3-n/3。其中：下標(biāo)“1”表示上一個(gè)時(shí)間層的值；下標(biāo)“2”表示下一個(gè)時(shí)間層的值；N表示網(wǎng)格點(diǎn)數(shù)；n表示空間離散精度。計(jì)算得到的數(shù)值結(jié)果如表1所示。結(jié)果表明五階HONCS格式達(dá)到了設(shè)計(jì)精度，驗(yàn)證了RAD型Riemann求解器能夠應(yīng)用于高階精度格式。

表1 五階HONCS格式的數(shù)值精度驗(yàn)證Table 1 Numerical orders of accuracy for HONCS5 scheme

2.3 頻譜分析

本文采用Pirozzoli[34]的ADR（approximate dispersion relation）方法來(lái)分析五階HONCS格式的分辨率和耗散。該方法現(xiàn)在被廣泛地應(yīng)用于非線(xiàn)性激波捕捉格式的頻譜分析[35-36]。計(jì)算得到的修正波數(shù)實(shí)部對(duì)應(yīng)于格式的色散（分辨率），而虛部對(duì)應(yīng)于格式的耗散。

為了更清楚地說(shuō)明五階HONCS格式的頻譜特性，對(duì)比了經(jīng)典五階WENO格式[37]。從圖3、圖4中可以看出，五階HONCS格式的分辨率高于五階WENO格式的分辨率，并且五階HONCS格式的耗散遠(yuǎn)小于五階WENO格式。

圖3 數(shù)值格式的色散特性Fig.3 Dispersion of numerical schemes

圖4 數(shù)值格式的耗散特性Fig.4 Dissipation of numerical schemes

3 數(shù)值實(shí)驗(yàn)與分析

空間離散采用五階HONCS格式結(jié)合各類(lèi)Riemann求解器，時(shí)間離散采用三階TVD Runge-Kutta格式，計(jì)算了Titarev-Toro問(wèn)題、激波衍射問(wèn)題、激波雙馬赫反射問(wèn)題、鈍頭體繞流問(wèn)題算例。對(duì)于二維問(wèn)題，也給出了一階迎風(fēng)格式的計(jì)算結(jié)果，其目的是為了更清晰地對(duì)比Riemann求解器計(jì)算激波的穩(wěn)定性。需要說(shuō)明的是，Einfeldt等[8]證明了HLLEM求解器和Roe求解器的區(qū)別僅僅在于數(shù)值信號(hào)速度的不同。在本文的數(shù)值算例中，所有Riemann求解器采用的數(shù)值信號(hào)速度是相同的，因此HLLEM求解器和Roe求解器是等價(jià)的，本文只給出HLLEM求解器的計(jì)算結(jié)果。

3.1 Titarev-Toro問(wèn)題[38]

該問(wèn)題流場(chǎng)中含有豐富的密度波結(jié)構(gòu)，能夠考察計(jì)算格式對(duì)流場(chǎng)細(xì)節(jié)的分辨能力，是驗(yàn)證數(shù)值格式的標(biāo)準(zhǔn)算例。該問(wèn)題的初始條件為：

計(jì)算區(qū)域取為[-5, 5]，計(jì)算網(wǎng)格采用1 000個(gè)點(diǎn)，初始間斷位于x=-4.5處，最終計(jì)算時(shí)刻取t= 5。

采用五階HONCS格式結(jié)合四種不同的Riemann求解器進(jìn)行計(jì)算，密度波的分布如圖5所示。從計(jì)算結(jié)果可以看出，新型RAD算子與HLLC算子和HLLEM算子的表現(xiàn)十分接近，這是因?yàn)楸締?wèn)題的密度波不是間斷的，流場(chǎng)比較光滑，使得RAD算子的效果接近于HLLEM算子，從而說(shuō)明RAD算子的混合因子取值是合理的。HLL算子因?yàn)楹纳⒈容^大，在本問(wèn)題中對(duì)高頻密度波的分辨率結(jié)果相對(duì)較差。

圖5 HONCS5格式求解Titarev-Toro問(wèn)題Fig.5 Solutions of the Titarev-Toro problem obtained by HONCS5 scheme

3.2 激波衍射問(wèn)題[14]

激波衍射是由Quirk[14]提出的膨脹波問(wèn)題，被廣泛用來(lái)驗(yàn)證Riemann求解器計(jì)算激波的穩(wěn)定性[21，23]。該問(wèn)題描述的是馬赫數(shù)5.09的激波從90°拐角的臺(tái)階角點(diǎn)處運(yùn)動(dòng)，激波沿x方向從左向右傳播，計(jì)算區(qū)域取[0, 1]× [0, 1]。激波波前靜止空氣的密度為1.4，壓力為1，波后按運(yùn)動(dòng)激波關(guān)系給定初始條件。初始物理參數(shù)為：

左邊界在[0, 0.5]的區(qū)域給定反射壁面邊界條件，在(0.5, 1]的區(qū)域給定波后值。上邊界區(qū)域根據(jù)激波運(yùn)動(dòng)所在的位置給定波前值或波后值。下邊界和右邊界給定波前值。本文采用800 × 800的均勻網(wǎng)格，計(jì)算終止時(shí)刻為t= 0.18。

圖6和圖7分別給出了一階迎風(fēng)格式和五階HONCS格式的計(jì)算結(jié)果，取密度為0.5～7.3共30條等值線(xiàn)作圖。

圖6 一階迎風(fēng)格式求解激波衍射問(wèn)題Fig.6 Numerical results of the shock diffraction problem obtained by first-order upwind scheme

圖7 HONCS5格式求解激波衍射問(wèn)題Fig.7 Numerical results of the shock diffraction problem obtained by HONCS5 scheme

因?yàn)閺?qiáng)激波附近是間斷區(qū)域，根據(jù)混合因子的定義可知，新型RAD算子在激波附近的計(jì)算效果接近于HLL算子。從計(jì)算結(jié)果可以看出，新型RAD算子和HLL算子都能夠準(zhǔn)確地計(jì)算運(yùn)動(dòng)激波、膨脹波和再生二次激波。

對(duì)于一階迎風(fēng)格式，從圖6中可以看出，采用HLLC算子和HLLEM算子會(huì)導(dǎo)致計(jì)算正激波出現(xiàn)紅玉現(xiàn)象。對(duì)于五階HONCS格式，從圖7中可以看出，采用HLLC算子會(huì)出現(xiàn)嚴(yán)重的紅玉現(xiàn)象。而采用HLLEM算子會(huì)出現(xiàn)負(fù)密度，導(dǎo)致計(jì)算發(fā)散，無(wú)法得到最終解。從圖6和圖7的對(duì)比中還可以發(fā)現(xiàn)，高階五階HONCS格式對(duì)于接觸面的分辨率明顯優(yōu)于迎風(fēng)格式。

3.3 激波雙馬赫反射問(wèn)題[39]

雙馬赫反射問(wèn)題包含強(qiáng)激波和滑移線(xiàn)，非常適合于考察格式的激波捕捉能力和流場(chǎng)精細(xì)結(jié)構(gòu)的分辨率。本問(wèn)題描述的是馬赫數(shù)為10的強(qiáng)運(yùn)動(dòng)斜激波以與x軸方向呈60°角的方向入射，入射點(diǎn)在（1/6, 0），計(jì)算區(qū)域取[0, 4]× [0, 1]。激波波前靜止空氣的密度為1.4，壓力為1，波后按激波關(guān)系給定初始條件。初始物理參數(shù)為：

下邊界在[1/6, 4]的區(qū)域給定壁面反射邊界條件，其他邊界按照激波運(yùn)動(dòng)所在的位置分別給定波前或波后的值。采用1 920×480的均勻網(wǎng)格計(jì)算到無(wú)量綱時(shí)間t= 0.2。

圖8和圖9分別給出了一階迎風(fēng)格式和五階HONCS格式的計(jì)算結(jié)果，取密度為1.731～20.92共30條等值線(xiàn)作圖。從計(jì)算結(jié)果可以看出，對(duì)于一階迎風(fēng)格式和五階HONCS格式，新型RAD算子和HLL算子都能夠準(zhǔn)確地計(jì)算激波的馬赫桿。對(duì)于一階迎風(fēng)格式，采用HLLC算子和HLLEM算子都會(huì)導(dǎo)致彎曲馬赫桿現(xiàn)象[14]，HLLEM算子的計(jì)算結(jié)果比HLLC算子更不穩(wěn)定。對(duì)于五階HONCS格式，從圖9中可以看出，采用HLLC算子在馬赫桿處會(huì)出現(xiàn)紅玉現(xiàn)象，而采用HLLEM算子會(huì)導(dǎo)致計(jì)算發(fā)散。從圖8和圖9的對(duì)比中還可以發(fā)現(xiàn)，高階五階HONCS格式對(duì)于滑移線(xiàn)的分辨率明顯優(yōu)于迎風(fēng)格式。

圖8 一階迎風(fēng)格式求解雙馬赫反射問(wèn)題Fig.8 Numerical results of double Mach reflection problem for first-order upwind scheme

圖9 HONCS5格式求解雙馬赫反射問(wèn)題Fig.9 Numerical results of the double Mach reflection problem obtained by HONCS5 scheme

滑移線(xiàn)附近是相對(duì)光滑的區(qū)域，根據(jù)混合因子的機(jī)制可知，新型RAD算子在滑移線(xiàn)附近耗散應(yīng)當(dāng)?shù)陀贖LL算子。圖9的計(jì)算結(jié)果驗(yàn)證了新型RAD算子比HLL算子在滑移線(xiàn)附近的分辨率高。

3.4 鈍頭體繞流問(wèn)題

超聲速鈍頭體繞流問(wèn)題是檢驗(yàn)數(shù)值格式是否會(huì)遭遇紅玉現(xiàn)象的典型算例。該問(wèn)題的紅玉現(xiàn)象通常是指鈍頭體繞流的弓形激波在駐點(diǎn)線(xiàn)附近發(fā)生異常的凸起。平頭鈍頭體高度為0.4，長(zhǎng)度為3.4，前緣位置x= 0.6，中心線(xiàn)位置y= 0。計(jì)算域?yàn)?[0, 4]× [-1, 1]。初始條件為馬赫數(shù)3的強(qiáng)運(yùn)動(dòng)激波沿著x方向向右傳播。流場(chǎng)初始物理參數(shù)為：

左邊界取來(lái)流條件，上下邊界和右邊界為出流邊界條件，鈍頭體壁面采用滑移壁面邊界條件。采用800 × 200的均勻網(wǎng)格計(jì)算截止到無(wú)量綱時(shí)間t= 5的時(shí)刻。

圖10和圖11分別給出了一階迎風(fēng)格式和五階HONCS格式的計(jì)算結(jié)果，取馬赫數(shù)0.2～3.7共30條等值線(xiàn)作圖。

圖11 HONCS5格式求解鈍頭體繞流問(wèn)題Fig.11 Flow fields around a blunt body obtained by HONCS5 scheme

鈍頭體繞流的強(qiáng)激波計(jì)算需要魯棒的數(shù)值方法來(lái)克服數(shù)值振蕩，由于五階HONCS格式的耗散較小，需要耗散較大的Riemann算子來(lái)保證計(jì)算方法的穩(wěn)定性。對(duì)于本問(wèn)題，新型RAD算子參數(shù)調(diào)整為δ=1×10-10。

從計(jì)算結(jié)果可以看出，對(duì)于一階迎風(fēng)格式和五階HONCS格式，新型RAD算子和HLL算子都能夠準(zhǔn)確地捕捉弓形激波，沒(méi)有任何的非物理振蕩現(xiàn)象。從圖10中可以發(fā)現(xiàn)，對(duì)于一階迎風(fēng)格式，采用HLLC算子和HLLEM算子都會(huì)導(dǎo)致紅玉現(xiàn)象，即沿著駐點(diǎn)線(xiàn)附近出現(xiàn)了明顯的異常凸起結(jié)構(gòu)，并且HLLEM算子比HLLC算子的紅玉現(xiàn)象更嚴(yán)重。對(duì)于五階HONCS格式，采用HLLC算子在弓形激波處出現(xiàn)紅玉現(xiàn)象，而采用HLLEM算子會(huì)導(dǎo)致計(jì)算發(fā)散。

4 結(jié) 論

針對(duì)近似Riemann求解器在計(jì)算激波時(shí)的數(shù)值不穩(wěn)定性問(wèn)題，通過(guò)合理設(shè)計(jì)反擴(kuò)散矩陣，發(fā)展出一類(lèi)新型RAD近似Riemann求解器。

通過(guò)適當(dāng)?shù)淖儞Q，可以發(fā)現(xiàn)傳統(tǒng)的Roe、HLL、HLLEM等Riemann求解器都是RAD求解器的某種特殊形式。本文構(gòu)造了自適應(yīng)混合因子，得到的新型RAD算子能夠抑制HLLEM算子和HLLC算子出現(xiàn)的計(jì)算激波不穩(wěn)定現(xiàn)象。根據(jù)混合因子的機(jī)制，新型RAD算子比HLL算子在滑移線(xiàn)等光滑區(qū)域的耗散小。

新型RAD求解器不但能夠應(yīng)用于低階迎風(fēng)格式，并且能夠廣泛應(yīng)用到高階格式，同時(shí)可以保持原有差分格式的高階精度。

數(shù)值實(shí)驗(yàn)結(jié)果表明，新的RAD求解器克服了傳統(tǒng)近似Riemann求解器的缺陷，既能精確捕捉接觸間斷和激波，又能大幅提高對(duì)剪切層等精細(xì)結(jié)構(gòu)的分辨率。

后續(xù)的研究中，將針對(duì)黏性流動(dòng)問(wèn)題驗(yàn)證新型RAD求解器的適用性，同時(shí)探索更好的自適應(yīng)因子，避免對(duì)特定問(wèn)題需要調(diào)整經(jīng)驗(yàn)參數(shù)的問(wèn)題。