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

(1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000；2.中國空氣動力研究與發(fā)展中心，綿陽 621000；3.北京航空航天大學 能源與動力工程學院，北京 100191)

0 引 言

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

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

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

混合方法能夠解決激波計算不穩(wěn)定性問題，同時也能避免熵修正方法中自由參數的經驗性，因此混合方法在實際計算中比較常用，但必須選擇適當的基礎Riemann求解器和混合因子。近年來，Dumbser等[22]和Xie等[23]基于HLLEM算子的反擴散矩陣來構造混合Riemann求解器，他們認為以HLLEM算子作為基礎，對接觸波和剪切波分量進行修正的混合方法比較合理。

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

1 數值格式的構造

1.1 控制方程和差分格式

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

其中，守恒變量為：

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

本文的數值方法主要是離散Euler方程的對流項，下面以 ?F/?x的離散為例進行說明。

選取圖1的網格模板，網格結點處的一階導數由線性中心緊致格式[24]來得到。

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

為了能夠計算包含激波等復雜流動的非線性問題，緊致格式（4）中可以采用Riemann求解器來近似。通過半結點處的通量把迎風和耗散性引入到差分格式，下標L表示半結點左側方向，R表示半結點右側方向。其中半結點處的值和由混合加權非線性插值[25]得到。因為的系數關于與是對稱的，所以只需給出的構造方法。

其中非線性權為：

其中線性權為：

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

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

1.2 數值信號速度

由于數值信號速度在Riemann求解器的計算中起著影響耗散性的重要作用，同時構造所有的Riemann求解器都要先指定數值信號速度。因此在說明Riemann求解器之前，先給出幾種常用的數值信號速度。

1.2.1 Roe類型[5]

1.2.2 Einfeldt 1988類型[7]

1.2.3 Einfeldt 1991類型[8]

1.2.4 Batten類型[26]

1.2.5 Davis類型[27]

1.2.6 Toro類型[28]

其中K=L或者K=R。

數值信號速度為：SL=qL-ηLcL，SR=qR+ηRcR。

1.3 傳統Riemann求解器

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

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

則有如下關系式：

1.3.1 Roe 型Riemann求解器[5]

特征值為：

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

直接采用原始Roe算子求解激波問題，通常會出現熵違反解[28]，必須對特征值進行熵修正才能得到正確解。但是一般來說，熵修正的經驗性很強，并且特定問題必須采用特定的熵修正值，因此目前還沒有標準的做法。

1.3.2 HLL型Riemann求解器[6]

HLL通量用S-、S+可以表達為：

把關系式（9）代入式（15），則可簡潔地表達為：

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

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

其中：

1.3.4 HLLEM型Riemann求解器[8]

以二維問題為例，HLLEM通量[8]為：

把關系式（9）代入式（21），則可表示為：

其中：

反擴散系數的定義為：

1.3.5 通量分裂型Riemann求解器

一些學者[29-31]采用通量分裂的方式來構造Riemann算子，即通過來計算數值通量。但是Steger-Warming分裂[2]、Van Leer分裂[3]等是逐點分裂的，無法表示成Riemann算子的通用形式。只有局部Lax-Friedrichs分裂可以變換成通用形式其中局部分裂系數為表示局部分裂點坐標。通量分裂型Riemann算子不是本文關注的重點，因此不做詳細比較和討論。

1.4 新型自適應耗散RAD型Riemann求解器

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

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

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

流場光滑量度為：

本文中取無量綱參數 ε =1×10-40。本文中若無特殊說明，則無量綱參數 δ =1×10-4。

顯然，矩陣D的特征值可以寫成：

當S+=SR＞0且S-=SL＜0時，即流動是亞聲速時，令則可以得到反擴散矩陣

通過對比RAD求解器與Roe、HLL和HLLEM等Riemann求解器的公式，可以得到如下的等價關系：

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

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

2 格式精度驗證和頻譜分析

2.1 數值信號速度的選擇

為了確定在Riemann算子中的數值信號速度，本文采用修正Sod激波管問題[26,28]來驗證和比較計算結果。該問題包含激波、膨脹波、接觸間斷等，能夠很好的評估數值方法的熵滿足特性。計算網格采用100個點，計算到無量綱時間t= 0.2，間斷左右兩邊參數為：

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

圖2 數值信號速度的比較Fig.2 Comparison of numerical signal velocities

2.2 精度驗證

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

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

表1 五階HONCS格式的數值精度驗證Table 1 Numerical orders of accuracy for HONCS5 scheme

2.3 頻譜分析

本文采用Pirozzoli[34]的ADR（approximate dispersion relation）方法來分析五階HONCS格式的分辨率和耗散。該方法現在被廣泛地應用于非線性激波捕捉格式的頻譜分析[35-36]。計算得到的修正波數實部對應于格式的色散（分辨率），而虛部對應于格式的耗散。

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

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

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

3 數值實驗與分析

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

3.1 Titarev-Toro問題[38]

該問題流場中含有豐富的密度波結構，能夠考察計算格式對流場細節(jié)的分辨能力，是驗證數值格式的標準算例。該問題的初始條件為：

計算區(qū)域取為[-5, 5]，計算網格采用1 000個點，初始間斷位于x=-4.5處，最終計算時刻取t= 5。

采用五階HONCS格式結合四種不同的Riemann求解器進行計算，密度波的分布如圖5所示。從計算結果可以看出，新型RAD算子與HLLC算子和HLLEM算子的表現十分接近，這是因為本問題的密度波不是間斷的，流場比較光滑，使得RAD算子的效果接近于HLLEM算子，從而說明RAD算子的混合因子取值是合理的。HLL算子因為耗散比較大，在本問題中對高頻密度波的分辨率結果相對較差。

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

3.2 激波衍射問題[14]

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

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

圖6和圖7分別給出了一階迎風格式和五階HONCS格式的計算結果，取密度為0.5～7.3共30條等值線作圖。

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

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

因為強激波附近是間斷區(qū)域，根據混合因子的定義可知，新型RAD算子在激波附近的計算效果接近于HLL算子。從計算結果可以看出，新型RAD算子和HLL算子都能夠準確地計算運動激波、膨脹波和再生二次激波。

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

3.3 激波雙馬赫反射問題[39]

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

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

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

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

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

滑移線附近是相對光滑的區(qū)域，根據混合因子的機制可知，新型RAD算子在滑移線附近耗散應當低于HLL算子。圖9的計算結果驗證了新型RAD算子比HLL算子在滑移線附近的分辨率高。

3.4 鈍頭體繞流問題

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

左邊界取來流條件，上下邊界和右邊界為出流邊界條件，鈍頭體壁面采用滑移壁面邊界條件。采用800 × 200的均勻網格計算截止到無量綱時間t= 5的時刻。

圖10和圖11分別給出了一階迎風格式和五階HONCS格式的計算結果，取馬赫數0.2～3.7共30條等值線作圖。

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

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

從計算結果可以看出，對于一階迎風格式和五階HONCS格式，新型RAD算子和HLL算子都能夠準確地捕捉弓形激波，沒有任何的非物理振蕩現象。從圖10中可以發(fā)現，對于一階迎風格式，采用HLLC算子和HLLEM算子都會導致紅玉現象，即沿著駐點線附近出現了明顯的異常凸起結構，并且HLLEM算子比HLLC算子的紅玉現象更嚴重。對于五階HONCS格式，采用HLLC算子在弓形激波處出現紅玉現象，而采用HLLEM算子會導致計算發(fā)散。

4 結 論

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

通過適當的變換，可以發(fā)現傳統的Roe、HLL、HLLEM等Riemann求解器都是RAD求解器的某種特殊形式。本文構造了自適應混合因子，得到的新型RAD算子能夠抑制HLLEM算子和HLLC算子出現的計算激波不穩(wěn)定現象。根據混合因子的機制，新型RAD算子比HLL算子在滑移線等光滑區(qū)域的耗散小。

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

數值實驗結果表明，新的RAD求解器克服了傳統近似Riemann求解器的缺陷，既能精確捕捉接觸間斷和激波，又能大幅提高對剪切層等精細結構的分辨率。

后續(xù)的研究中，將針對黏性流動問題驗證新型RAD求解器的適用性，同時探索更好的自適應因子，避免對特定問題需要調整經驗參數的問題。