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

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

0 引 言

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

Riemann求解器是數值求解雙曲系統(tǒng)的重要組成部分。對于守恒形式的雙曲系統(tǒng)，已有許多學者提出了多種著名的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求解器依然可能在計算激波時出現(xiàn)不穩(wěn)定。為了克服HLL算子無法模擬接觸間斷的缺陷，Toro等[9]提出了HLLC型Riemann求解器，恢復了HLL算子中丟失的接觸波和剪切波。HLLC算子保留了HLL算子的熵增特性和正定保持性質，但在計算激波中仍然會出現(xiàn)非物理的數值振蕩?；谟上嗫臻g中的路徑積分來得到數值耗散項的構造思路， Osher等[10]提出了Osher型Riemann求解器，其優(yōu)點是數值通量光滑并且滿足熵增特性。在具體計算時，HLLEM求解器與Roe求解器僅需求解特征值和相關的左右特征向量，而Osher型Riemann求解器需要計算整套的特征系統(tǒng)，因此Osher求解器計算量較大。為了將通量中的線性項和非線性項加以區(qū)別，Liou[11]等提出了一種結合通量分裂和Riemann求解器的AUSM系列格式，這種方法把壓力項和對流項進行分裂，對部分問題能夠消除計算激波的不穩(wěn)定性現(xiàn)象。

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

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

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

則有如下關系式：

1.3.1 Roe 型Riemann求解器[5]

特征值為：

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

直接采用原始Roe算子求解激波問題，通常會出現(xiàn)熵違反解[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算子中的反擴散系數在激波附近取值較大，導致在計算激波時出現(xiàn)不穩(wěn)定現(xià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求解器的公式，可以得到如下的等價關系：

因此，傳統(tǒng)的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格式的分辨率和耗散。該方法現(xiàn)在被廣泛地應用于非線性激波捕捉格式的頻譜分析[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算子的表現(xiàn)十分接近，這是因為本問題的密度波不是間斷的，流場比較光滑，使得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算子會導致計算正激波出現(xiàn)紅玉現(xiàn)象。對于五階HONCS格式，從圖7中可以看出，采用HLLC算子會出現(xiàn)嚴重的紅玉現(xiàn)象。而采用HLLEM算子會出現(xiàn)負密度，導致計算發(fā)散，無法得到最終解。從圖6和圖7的對比中還可以發(fā)現(xiàn)，高階五階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算子都會導致彎曲馬赫桿現(xiàn)象[14]，HLLEM算子的計算結果比HLLC算子更不穩(wěn)定。對于五階HONCS格式，從圖9中可以看出，采用HLLC算子在馬赫桿處會出現(xiàn)紅玉現(xiàn)象，而采用HLLEM算子會導致計算發(fā)散。從圖8和圖9的對比中還可以發(fā)現(xiàn)，高階五階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 鈍頭體繞流問題

超聲速鈍頭體繞流問題是檢驗數值格式是否會遭遇紅玉現(xiàn)象的典型算例。該問題的紅玉現(xiàn)象通常是指鈍頭體繞流的弓形激波在駐點線附近發(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算子都能夠準確地捕捉弓形激波，沒有任何的非物理振蕩現(xiàn)象。從圖10中可以發(fā)現(xiàn)，對于一階迎風格式，采用HLLC算子和HLLEM算子都會導致紅玉現(xiàn)象，即沿著駐點線附近出現(xiàn)了明顯的異常凸起結構，并且HLLEM算子比HLLC算子的紅玉現(xiàn)象更嚴重。對于五階HONCS格式，采用HLLC算子在弓形激波處出現(xiàn)紅玉現(xiàn)象，而采用HLLEM算子會導致計算發(fā)散。

4 結 論

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

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

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

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

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