李妍慧, 陳琮巍, 任玉新,*, 孫振生, 王秋菊

(1. 清華大學(xué) 航天航空學(xué)院, 北京 100084; 2. 火箭軍工程大學(xué), 西安 710025; 3. 北京應(yīng)用物理與計算數(shù)學(xué)研究所, 北京 100094)

0 引 言

超聲速流中，多尺度結(jié)構(gòu)(如湍流)、間斷(如激波)以及這些結(jié)構(gòu)之間的相互作用同時存在，給數(shù)值模擬帶來很大的挑戰(zhàn)。一方面，多尺度結(jié)構(gòu)要求數(shù)值格式具有較小的色散耗散誤差以精確地捕捉幅值和相位，另一方面，間斷附近需要格式具有足夠的耗散以抑制數(shù)值振蕩。這兩種矛盾的要求使得構(gòu)造高分辨率的激波捕捉格式非常困難。目前最具代表性的高精度激波捕捉格式是WENO格式，它最初由Liu等[1]提出，之后被Jiang和Shu[2]改進。WENO格式可以根據(jù)流場局部的光滑性，動態(tài)調(diào)整候選模板的權(quán)重來無振蕩地捕捉激波，但是這一非線性機制也使得格式在光滑區(qū)的耗散過大，特別是在計算可壓縮流動的多尺度結(jié)構(gòu)時，數(shù)值耗散可能會高于物理耗散[3]。

解決這一問題的一種方法是改進WENO或ENO格式的非線性機制，更好地實現(xiàn)光滑區(qū)和間斷的尺度分離。Henrick等[4]提出了WENO-M格式，通過映射非線性權(quán)解決了WENO-JS在極點處不滿足五階收斂條件的問題。Borgers等[5]提出的WENO-Z格式采用了一種歸一化的光滑因子，可以使非線性權(quán)在光滑區(qū)更快地恢復(fù)到線性權(quán)。該格式在極點處也滿足五階收斂條件，且計算量遠小于WENO-M格式。劉朋欣等[6]參考Henrick等[4]提出的映射函數(shù)的思想，將一種五次分段多項式映射函數(shù)應(yīng)用于一種由NND格式改型得到的中心型三階格式，構(gòu)造出SWENO3-PPM5格式，算例測試結(jié)果表明該格式具有比WENO3更高的分辨率。Taylor等[7]基于一種相對光滑限制器構(gòu)造得到WENO-RL格式，當光滑因子間的比值較小時直接采用線性權(quán)，這種方法可以使格式在中低波數(shù)段的色散耗散特性和對應(yīng)的線性格式一致，從而改善格式的譜特性。WENO-SYMBO格式[8]將對應(yīng)的線性格式由迎風(fēng)改為中心格式，并通過求解一個積分形式目標函數(shù)的最小值，對格式的譜特性進行優(yōu)化。Hu和Adams等[9]發(fā)展了一種自適應(yīng)的中心迎風(fēng)六階WENO格式(WENO-CU6)，采用新的光滑因子控制格式在中心和迎風(fēng)之間切換。最近發(fā)展的TENO格式[10]采用新的模板選擇和加權(quán)策略，來實現(xiàn)光滑區(qū)和間斷的有效區(qū)分，從而改善分辨率。

另外一種方法是利用間斷探測器，將激波捕捉格式和譜特性良好的格式結(jié)合起來，構(gòu)造混合格式?；旌细袷皆诹鲌龉饣瑓^(qū)采用低耗散的格式以提高分辨率，而在間斷附近采用激波捕捉格式以避免數(shù)值振蕩?；旌细袷绞紫扔葾dams和Shariff[11]提出，他們將非守恒形式的緊致格式與ENO格式相結(jié)合，用于激波湍流相互作用流場的模擬。Pirozzoli[12]將守恒的緊致格式與WENO格式相結(jié)合，構(gòu)造了一種守恒的緊致-WENO混合格式。該工作由Ren等[13]進一步發(fā)展，他設(shè)計了一種新的間斷探測器，并采用加權(quán)的方式將守恒的緊致格式和WENO格式混合，避免了兩種格式的突然切換。武從海等[14]將Ren等[13]提出的間斷探測器進行改進，構(gòu)造了基于WENO-Z格式和五階守恒緊致格式的混合格式，并且通過僅對混合格式中WENO部分進行特征投影處理，提高了計算效率。對于混合格式，間斷探測器的精度對格式的譜特性有很大影響。傳統(tǒng)的間斷探測器通常采用流場解的一階導(dǎo)數(shù)或二階導(dǎo)數(shù)來判斷流場的光滑性[13,15-16]。此外，近年來還發(fā)展了基于非正交小波基函數(shù)[15,17]以及WENO非線性權(quán)[18-21]的間斷探測器。然而，如何將間斷和高頻成分以及極點進行區(qū)分仍然是一個困難的問題。

上述兩種方法改善計算效果的主要思路都是在流場光滑區(qū)盡可能地優(yōu)先采用線性格式，因此，線性格式的性能對最終計算結(jié)果有直接影響。例如，WENO-M和WENO-Z等的線性格式是五階迎風(fēng)格式，而WENO-SYMBO格式對應(yīng)的線性格式是六階中心格式?；旌细袷降木€性部分則多采用緊致格式和中心格式。當線性格式是迎風(fēng)格式時，對湍流直接數(shù)值模擬等應(yīng)用而言，耗散還是較大；當線性格式為中心格式時，耗散則不足。Pirozzoli[12]指出，為了減少數(shù)值振蕩，格式需要具有一定的耗散以抑制色散誤差較大的高波數(shù)成分。但是，如何合理地確定線性格式的耗散大小仍然是一個開放問題。

我們最近的工作，就是解決混合格式中線性部分的色散優(yōu)化和耗散控制問題。關(guān)于色散優(yōu)化的方法，已經(jīng)開展了較多研究，但是線性格式的耗散控制相關(guān)的研究還比較少。本文綜述了我們近期的相關(guān)研究工作。首先，Sun等[22]提出了一類色散耗散相互獨立的有限差分格式，在優(yōu)化色散同時，耗散可以通過一個參數(shù)調(diào)節(jié)，從而得到了色散最小、耗散可控(Minimized Dispersion and Controllable Dissipation, MDCD)格式。但是MDCD格式耗散參數(shù)的設(shè)置依賴于具體問題，需要基于經(jīng)驗人工進行選取，而在大型的數(shù)值計算如直接數(shù)值模擬中，參數(shù)的調(diào)試會花費大量的時間。Hu等[23]認為，格式的耗散大小應(yīng)該與色散誤差相關(guān)。他們據(jù)此得到一個色散耗散條件，來確定足夠的耗散以抑制非物理的高頻振蕩。Sun等[24]根據(jù)這一關(guān)系，確定了六階MDCD格式的耗散參數(shù)，在多個測試算例中，格式均表現(xiàn)出較好的魯棒性。這種方法雖然可以得到具有合適耗散的格式，但它只能調(diào)節(jié)格式耗散的整體水平，而沒有考慮解的局部尺度。

理想的格式耗散應(yīng)該是與數(shù)值解的局部尺度相關(guān)的。在解的局部尺度遠大于網(wǎng)格尺度時，線性格式的耗散應(yīng)該很小甚至為0。當解的局部尺度與網(wǎng)格尺度接近時，此時格式的色散誤差已非常大，應(yīng)當將耗散增加到足以抑制高波數(shù)振蕩的水平。根據(jù)這一思想，Li等[25]發(fā)展了一類具有最小色散、自適應(yīng)耗散特性(Minimized Dispersion and Adaptive Dissipation, MDAD)的格式。

自適應(yīng)耗散的差分格式包含以下要素。首先，為了保證耗散的自適應(yīng)不影響格式的色散特性，選擇色散和耗散相互獨立的四階MDCD格式作為自適應(yīng)耗散格式的基礎(chǔ)。其次，為了識別數(shù)值解的局部尺度，Li等[25]提出了一種基于流場解的導(dǎo)數(shù)的尺度識別器，可用來得到數(shù)值解局部的等效無量綱波數(shù)。最后，采用Hu等[23]提出的色散耗散條件，設(shè)計出耗散參數(shù)與局部等效無量綱波數(shù)之間的關(guān)系，從而構(gòu)造得到具有自適應(yīng)耗散特性的MDAD格式。

為了使格式具有激波捕捉能力，通過將MDAD格式的數(shù)值通量表示為一系列候選模板上的低階多項式的線性組合，并引入WENO格式的非線性權(quán)，發(fā)展了MDAD-WENO格式。之后將MDAD格式與MDAD-WENO格式相結(jié)合，構(gòu)造了一種混合格式，記為MDAD-HY。采用近似色散關(guān)系[26]分析MDAD-HY格式與其他激波捕捉格式的近似譜特性，發(fā)現(xiàn)MDAD-HY格式相比基于MDCD的混合格式MDCD-HY，在中低波數(shù)段的耗散更小，而在高波數(shù)段有相近的耗散，說明格式在數(shù)值解尺度較大的區(qū)域分辨率較高，而在間斷區(qū)也有較好的魯棒性。構(gòu)造混合格式的一個關(guān)鍵問題是合理的間斷探測器。為此，Li等[25]將傳統(tǒng)間斷探測器[13]與尺度識別器相結(jié)合，提出了一種新的間斷探測器，能夠有效減少對光滑區(qū)的極點的誤判。

本文介紹了MDCD格式、MDAD格式、混合格式MDCD-HY和MDAD-HY，以及MDCD格式的應(yīng)用情況和MDAD格式的測試算例，最后作出了總結(jié)。

1 色散最小、耗散可控有限差分格式

本文以一維線性波動方程為例，介紹高精度有限差分格式的色散優(yōu)化和耗散控制方法。方程形式如下:

其中,u為守恒量，f=au為通量且a為常數(shù)。考慮在均勻網(wǎng)格上，采用半離散有限差分格式對式(1)在(xj,t)處進行離散，得到:

其中fj+m=f(xj+mΔx)。式(3)可改寫為守恒形式，即:

常規(guī)的數(shù)值格式色散和耗散特性是相互關(guān)聯(lián)的，但是MDCD格式具有色散、耗散相互獨立的特性，這就使得它可以對色散進行單獨優(yōu)化，而耗散也可以根據(jù)不同工況進行調(diào)節(jié)。為了簡潔，本文將只介紹四階MDCD格式[22]，六階MDCD格式的具體細節(jié)參見文獻[24]。

1.1 MDCD格式的半離散特性

當空間有限差分格式包含對稱的(2r+1)個模板點時，式(3)可寫為:

當r=3時，四階的MDCD格式[22]形式如下：

不論γdisp和γdiss如何選取，格式都有四階精度。

下面分析式(6)所示四階MDCD的半離散色散耗散特性?？紤]單一的簡諧波，

則空間一階偏導(dǎo)數(shù)有精確解

其中k=ωΔx為無量綱波數(shù)。將式(8)代入式(6)并寫為：

其中k′=R(k′)+iI(k′)為修正的無量綱波數(shù)，R(k′)和I(k′)為k′的實部和虛部，分別控制半離散格式的色散和耗散特性。則四階MDCD格式的半離散色散和耗散特性為:

I(k′)=γdiss(cos3k-6cos2k+15cosk-10)

(11)

顯然格式的色散特性只與γdisp有關(guān)，耗散特性只與γdiss有關(guān)，從而實現(xiàn)了對格式色散和耗散的單獨控制。

MDCD格式的色散和耗散特性是相互獨立的，因此可以對色散進行單獨優(yōu)化而不影響耗散。通常認為格式的色散誤差應(yīng)當在選定的標準下達到最小。文獻[8]設(shè)計了一種加權(quán)積分形式的優(yōu)化目標函數(shù):

(12)

其中e-νk用于控制不同波數(shù)下色散誤差的相對權(quán)重，通過求解使目標函數(shù)式(12)最小的優(yōu)化問題，可以得到最優(yōu)的色散參數(shù)γdisp。為了使格式在盡可能大的波數(shù)范圍內(nèi)都有較小的色散誤差，文獻[22]最終采用了ν=8下的最優(yōu)色散參數(shù),

γdisp=0.046 378 3

(13)

圖1 不同格式的色散特性

圖2 不同格式的色散誤差

對于耗散特性，為了使格式保持穩(wěn)定，所有波數(shù)下的耗散都應(yīng)非負。由式(11)所示，MDCD格式的耗散大小只與γdiss有關(guān)，虛部I(k′)可改寫為:

I(k′)=γdissg(cosk)

(14)

其中,

g(x)=4x3-12x2+12x-4

(15)

由于g′(x)=12(x-1)2≥0，為單調(diào)非減函數(shù)，易證g(cosk)≤0。因此當a>0時，如果γdiss≥0，則MDCD格式有非負耗散，格式是穩(wěn)定的，且隨著γdiss的增大，格式的耗散增加。當耗散參數(shù)取0時，MDCD為中心格式，此時格式分辨率最高，但是無法抑制高波數(shù)的非物理振蕩。

在實際計算中，正如文獻[12]提到的，格式需要具有一定的耗散。對于四階MDCD格式，Sun等[22]對不同算例采用試錯的方法確定合適的耗散參數(shù)γdiss進行計算。而在文獻[24]中，Sun等采用Hu等[23]提出的色散耗散條件確定了六階MDCD格式的耗散參數(shù)。多個數(shù)值實驗結(jié)果顯示該六階MDCD格式表現(xiàn)出較好的魯棒性，沒有出現(xiàn)明顯的數(shù)值振蕩。Hu等[23]提出的色散耗散條件具體形式如下:

Hu等認為格式的耗散應(yīng)當與色散誤差相關(guān)，在式(16)中，給定波數(shù)下的色散耗散比值r越小表示格式對該波數(shù)的波耗散越大。當r<10時，格式的耗散基本可以抑制非物理的高頻波動。這里我們采用與文獻[24]相同的做法確定四階MDCD格式的通用耗散參數(shù)γdiss?？紤]到k=π時格式的耗散最大，我們采用r(π)來確定γdiss。本文選擇r≤9作為限制條件，此時滿足條件的最小耗散參數(shù)為：

γdiss=0.012

(17)

圖3展示了不同耗散參數(shù)下的MDCD格式與UW5以及C6格式的半離散耗散特性對比，其中括號內(nèi)的數(shù)字為耗散參數(shù)的大小。可以看到C6和MDCD(0)均為零耗散，MDCD(0.035)的耗散特性與UW5接近，而MDCD(0.012)的耗散介于C6和UW5之間。MDCD(0.012)在中低波數(shù)段耗散較小，格式有較高的分辨率，在高波數(shù)段，格式的色散誤差變大時，也能提供足夠的耗散來抑制振蕩，有較好的魯棒性。因此，對于大多數(shù)工況，我們推薦采用MDCD(0.012)。在本文中，如無特殊說明，MDCD格式的耗散參數(shù)均采用γdiss=0.012。

圖3 不同格式的耗散特性

1.2 MDCD格式的全離散特性分析

上一節(jié)中我們提到，MDCD格式的半離散色散耗散特性相互獨立，分別受色散參數(shù)γdisp和耗散參數(shù)γdiss控制。但實際計算中，結(jié)果同時受到空間離散和時間離散的影響，因此分析格式的全離散特性是有必要的。

考慮線性對流方程：

其中c=aΔt/Δx為庫朗數(shù)，k=ωΔx為無量綱波數(shù)。而半離散方程式(4)的解為：

其中k′為1.1節(jié)中所述的無量綱修正波數(shù)。對比式(19)和式(20)，可以發(fā)現(xiàn)空間離散格式引入的誤差為：

rs=e-i(k′-k)c=eI(k′)ce-i(R(k′)-k)c

(21)

其中|rs|=eI(k′)c為半離散耗散誤差，而δs=(R(k′)-k)c為半離散色散誤差。若采用顯式龍格庫塔格式進行時間積分，則在t=Δt時刻有：

其中G(k,c)為放大因子，具體形式為一個s階多項式，其中s為龍格庫塔格式的級數(shù)。對比式(22)和式(19)，可以發(fā)現(xiàn)G(k,c)為精確的放大因子e-ikc的近似。因此，格式的全離散誤差為：

其中|rt|=|G(k,c)|為全離散耗散誤差，δt=-δG-kc為全離散色散誤差?？紤]到1.1節(jié)中分析半離散色散耗散特性時采用的是R(k′)和I(k′)，為了使全離散特性分析與半離散量級一致，對比|rs|和|rt|，以及δs和δt，我們得到：

其中kt為全離散情況下的等效無量綱修正波數(shù)。由式(24)可知，全離散的色散耗散特性均與庫朗數(shù)c有關(guān)。為了不失一般性，我們采用c=0.3、0.6、1.0三種條件，對全離散色散耗散特性關(guān)于γdisp和γdiss的靈敏度進行測試，結(jié)果如圖4～圖6所示。

(a) 不同γdisp下的全離散色散特性

(a) 不同γdisp下的全離散色散特性

(a) 不同γdisp下的全離散色散特性

可以看到，當庫朗數(shù)較小(c=0.3、0.6)時，格式的全離散色散耗散特性受時間離散影響尚不明顯，主要由空間離散決定[24]，此時全離散色散特性主要受色散參數(shù)γdisp控制，對耗散參數(shù)γdiss不敏感，而全離散耗散特性主要受耗散參數(shù)γdiss影響，對色散參數(shù)γdisp不敏感。但當庫朗數(shù)c=1.0時，時間離散帶來的誤差與空間離散相當，此時全離散色散耗散特性表現(xiàn)出同時受色散參數(shù)γdisp和耗散參數(shù)γdiss影響。但是色散參數(shù)γdisp對耗散的影響仍比耗散參數(shù)γdiss小得多，而耗散參數(shù)γdiss對色散的影響也比色散參數(shù)γdisp小得多。因此可以認為，在全離散意義下，MDCD格式仍是可用的。

2 色散最小、耗散自適應(yīng)有限差分格式

在1.1節(jié)中，我們采用Hu等[23]提出的色散耗散條件確定了一個通用的四階MDCD格式的耗散參數(shù)。但是我們也注意到，當流場尺度較大時，可以采用更小的數(shù)值耗散，因為此時物理耗散可被準確地預(yù)測以使計算穩(wěn)定。在某些情況下，對于可精確分辨的解，采用中心格式計算已經(jīng)足夠[9]。因此，當全場采用相同的耗散參數(shù)時，在流場尺度較大的光滑區(qū)，可能出現(xiàn)耗散過大的情況。

理想的耗散參數(shù)應(yīng)該隨著解的局部尺度而自動調(diào)節(jié)。為了實現(xiàn)這一目標，必須發(fā)展數(shù)值解的尺度識別器，并在此基礎(chǔ)上構(gòu)造具有最小色散、自適應(yīng)耗散特性(Minimized Dispersion and Adaptive Dissipation, MDAD)的格式[25]。本節(jié)綜述這方面工作的最新進展。

2.1 數(shù)值解尺度識別器

原則上，在解的局部尺度遠大于網(wǎng)格尺度時，線性格式的耗散應(yīng)該很小甚至為0，以提高格式的分辨率。而當解的局部尺度與網(wǎng)格尺度接近時，應(yīng)該將耗散增加到足以抑制高波數(shù)振蕩的水平。如果能夠?qū)崿F(xiàn)這一點，則MDCD格式在魯棒性不下降時，分辨率將得到進一步改善。為此，本節(jié)將提出一種數(shù)值解局部尺度識別器，用于識別數(shù)值解局部的等效長度尺度，作為構(gòu)造自適應(yīng)耗散格式的第一步。

數(shù)值解尺度的確定是一個開放性的問題，目前相關(guān)的研究很少。Li等[28]提出了一種尺度識別器，用于優(yōu)化一種五點非線性格式，以計算擴散方程中的二階導(dǎo)數(shù)項。而文獻[25]提出了一種新的數(shù)值解局部波數(shù)識別器，用于自適應(yīng)耗散格式的構(gòu)造。

設(shè)計尺度識別器的難點在于定量描述尺度時看似矛盾的兩個要求。一方面，考慮到我們關(guān)注的是解的局部特性，尺度識別器應(yīng)該作用于局部的物理空間。另一方面，為了將尺度識別器用于耗散控制，需要將其表示為波數(shù)的形式，而波數(shù)是傅里葉空間中的概念。為了滿足上述兩個條件，文獻[25]提出一種兩步的構(gòu)造校準方法，即首先基于物理空間構(gòu)造一種尺度識別器，之后再采用正弦波來校準尺度識別器中的參數(shù)。

第一步：構(gòu)造?？紤]光滑函數(shù)f(x)，為了估計其在x0附近的長度尺度，將其在x0處作Taylor展開，有：

當數(shù)值解的尺度較大時，低階導(dǎo)數(shù)對函數(shù)變化的影響比高階導(dǎo)數(shù)大。例如，當解為直線分布時，只有Δf1對解的變化有貢獻。而x0附近如果出現(xiàn)小尺度結(jié)構(gòu)，則說明高階導(dǎo)數(shù)項產(chǎn)生的影響足夠大。因此，原則上數(shù)值解的無量綱尺度可以用不同Δfp的比值來表示：

其中δ為局部無量綱尺度，且p>q。將式(26)代入式(27)得到：

其中β是待定參數(shù)?？紤]到四階MDCD格式數(shù)值通量的模板點為6個，可以數(shù)值計算的導(dǎo)數(shù)f(p)的階數(shù)為p=1～5，這里選取p=1～4構(gòu)造尺度識別器，但此時仍然有多種可能的p-q組合方式。為了解決這個問題，可以采用如下的校準步驟。

第二步：校準。在這一步中，采用正弦波來確定式(28)中的參數(shù)。若將式(28)定義的長度尺度校準為精確的正弦波無量綱波長，最簡單的方式是令p=q+2，此時p-q只剩如下兩種組合方式：

證明過程如下?？紤]單一正弦波，

f=Asin(ωx+φ)

(30)

可以得到它的各階導(dǎo)數(shù)：

f(1)=Aωcos(ωx+φ)

f(2)=-Aω2sin(ωx+φ)

f(3)=-Aω3cos(ωx+φ)

f(4)=Aω4sin(ωx+φ)

(31)

因此，當β=2π時，無量綱尺度δ即為精確的無量綱波長，

基于式(33)，定義解的等效無量綱波數(shù)k,

將式(29)代入式(34)，并令β=2π，則有：

需要注意的是，式(35)定義的等效無量綱波數(shù)并不是傅里葉空間中的波數(shù)，而是局部長度尺度的倒數(shù)，是一個物理空間中局部的定義。對于單一正弦波，它可以得到精確的正弦波無量綱波數(shù)。為了避免式(35)在極值點和拐點附近的奇異性，可以把該式的兩種情況相結(jié)合，從而有:

其中ε為防止分母為0的小量，取ε=1×10-3。

為了驗證尺度識別器在流場中的效果，Li等[25]采用如下四組指定的函數(shù)進行靜態(tài)測試，其中包括:

(a)f=sin(16πx) , -1≤x≤1

(d)f=sin(2πex+1x) , -1≤x≤1

(a)

由測試結(jié)果可知，尺度識別器能夠準確地識別標準正弦波的無量綱波數(shù)，對于幅值增加和頻率增加的正弦波識別效果也很好。式(38)通過各階導(dǎo)數(shù)的恰當組合，有效抑制了等效無量綱波數(shù)在極值點和拐點處的振蕩。而對于間斷，尺度識別器往往將其識別為一個高頻成分，這也與物理事實相吻合。

2.2 MDAD格式

2.1節(jié)已經(jīng)介紹了可以用于判斷流場局部等效無量綱波數(shù)的尺度識別器，為了構(gòu)造耗散自適應(yīng)的MDAD格式，需要根據(jù)等效無量綱波數(shù)對式(7)中的耗散參數(shù)γdiss進行動態(tài)的調(diào)整。此時，γdiss在全場不再是一個常數(shù)。為了更清楚地介紹MDAD格式的構(gòu)造過程，可以將式(7)改寫為：

Li等[25]采用Hu等[23]提出的色散耗散條件指導(dǎo)式(40)的設(shè)計。需要注意的是色散耗散條件是在傅里葉空間中的關(guān)系，然而自適應(yīng)耗散應(yīng)該在物理空間中實施。因此，使用該條件時需要格外小心。不過對于單一正弦波，等效無量綱波數(shù)即為精確的正弦波波數(shù)，此時色散耗散條件可以很容易的使用。從這個角度考慮，使用色散耗散條件設(shè)計式(40)也可以看作是一個校準的過程，即以單一正弦波為例，采用色散耗散條件構(gòu)造耗散參數(shù)與等效無量綱波數(shù)之間的關(guān)系，并應(yīng)用于更一般的情況。需要注意的是，式(40)也可以基于經(jīng)驗或者物理意義進行設(shè)計，這一問題將在接下來的工作中繼續(xù)研究。

1.1節(jié)中提到，當限制色散耗散比值r(π)≤9時，滿足條件的最小耗散參數(shù)為γdiss=0.012?？紤]到MDCD格式的色散和耗散特性均為已知函數(shù)，則當r固定時，可以反推不同波數(shù)下所需的耗散參數(shù)。將式(14)代入式(16)可以得到耗散參數(shù)與波數(shù)的關(guān)系式：

圖8 耗散參數(shù)與波數(shù)的關(guān)系(實線為自適應(yīng)耗散參數(shù)，虛線為固定的耗散參數(shù)0.012)

圖9 色散耗散比值r隨波數(shù)的變化曲線

3 混合格式

雖然MDCD格式和MDAD格式有很好的色散耗散特性，但它們不能直接用于含間斷流場的數(shù)值模擬。為了使這些格式具有激波捕捉的能力，Sun等[22,24]和Li等[25]的方案是首先將這些格式與WENO格式相結(jié)合，發(fā)展相應(yīng)的非線性格式，之后將線性格式與非線性格式相結(jié)合，構(gòu)造混合格式。由于MDCD格式與MDAD格式在混合格式構(gòu)造方面基本一致，本文我們只介紹MDAD格式對應(yīng)的混合格式MDAD-HY。本節(jié)將首先介紹MDAD-WENO格式，之后介紹一種新的間斷探測器用于MDAD-HY格式的構(gòu)造，并用近似色散關(guān)系(ADR)[26]對混合格式的譜特性進行分析。需要注意的是，引入自適應(yīng)耗散機制后，MDAD格式實際上是非線性的，但為了簡便，在后文我們?nèi)匀环Q其為線性格式。

3.1 MDAD-WENO格式

WENO格式[2]是一類代表性的非線性激波捕捉格式，在含間斷流場的數(shù)值模擬中有廣泛的應(yīng)用。WENO格式的基本思想是通過對一系列候選模板上低階多項式的非線性組合，實現(xiàn)無振蕩的激波捕捉能力。為了構(gòu)造MDAD-WENO格式，可以采用包含對稱的6個模板點的MDAD格式作為它的線性部分，這些模板點被劃分如圖10所示的四條候選模板{S0,S1,S2,S3}，則MDAD格式的數(shù)值通量可以改寫為:

圖10 MDAD格式數(shù)值通量的候選模板Sk

表1 線性權(quán)Ck關(guān)于γdisp和的表達式

其中qk與式(43)意義相同，ωk為WENO格式的非線性權(quán)，與對應(yīng)候選模板Sk上解的光滑性有關(guān)。MDAD-WENO的非線性權(quán)采用與Jiang和Shu[2]相同的形式:

其中Ck采用表1中的關(guān)系計算，因此MDAD-WENO格式也是耗散自適應(yīng)的。ε為防止分母為0的小量，p與候選模板自適應(yīng)的敏感性有關(guān)，通常取p=2，而βk為對應(yīng)候選模板上解的光滑因子:

代入每個候選模板得到:

(48)

然而需要注意的是，與WENO-JS不同，MDAD-WENO是采用對稱的模板，正如文獻[8]中所說，它的候選模板S3包含的全部為順風(fēng)模板點，當該模板的非線性權(quán)ω3大于其他非線性權(quán)時，有可能會引起格式的不穩(wěn)定，因此，在式(48)計算完所有非線性權(quán)之后，需要對β3增加額外的限制,

以此來限制非線性權(quán)ω3不大于其他線性權(quán)，保證格式的穩(wěn)定。

3.2 MDAD-HY格式

WENO格式能夠很好地捕捉流場中的間斷，但是相比對應(yīng)的線性格式，由于引入了非線性機制，色散耗散特性往往會變差[26]。為了在提高格式分辨率的同時仍然保持良好的間斷捕捉能力，一種思路是采用合適的方法將WENO格式與譜特性良好的線性格式進行混合，在流場光滑區(qū)采用線性格式，在間斷附近采用WENO格式。文獻[25]對MDAD格式和MDAD-WENO格式進行組合，構(gòu)造了具有耗散自適應(yīng)、低色散特性的四階有限差分混合格式(MDAD-HY)。

MDAD-HY格式的數(shù)值通量可表示為：

其中，

ε為防止分母為0的小量，

為了驗證新的間斷探測器的效果，我們將其與Harten[15]、Jameson等[16]以及Ren等[13]提出的基于一階導(dǎo)或者二階導(dǎo)的間斷探測器進行比較。Harten的間斷探測器同樣基于流場梯度的變化，具體形式如下:

(56)

其中,閾值ψc=0.3，ε為防止分母為0的小量。而Jameson提出的間斷探測器基于歸一化的二階導(dǎo)數(shù):

采用2.1節(jié)中的(b)算例對間斷探測器進行靜態(tài)測試，結(jié)果如圖11所示?？梢钥吹?，Harten和Jameson的間斷探測器能夠識別間斷，但是在正弦波光滑區(qū)都出現(xiàn)大量誤判，這將嚴重影響混合格式的分辨率。Ren提出的間斷探測器識別效果相對較好，但在光滑區(qū)部分極點附近也存在誤判。Li提出的新的間斷探測器不僅能夠準確識別間斷，而且由于引入了尺度識別器的結(jié)果進行截斷，光滑區(qū)的識別結(jié)果得到很大改善。

(a) Harten

圖12 間斷探測器在不同波數(shù)下測試得到的的平均值

3.3 MDAD-HY格式的譜特性

線性格式的色散耗散特性可由傅里葉分析得到，但是對于非線性格式，目前還沒有解析的分析譜特性的方法。Pirozzoli[26]提出一種近似色散關(guān)系(Approximate Dispersion Relation, ADR)，可以數(shù)值地計算非線性的激波捕捉格式的修正波數(shù)，從而得到非線性格式近似的色散耗散特性。本節(jié)將利用ADR分析MDAD-HY格式的色散耗散特性，并與原始的MDCD-HY格式[22]進行對比。

圖13展示了不同格式的近似耗散特性，Jiang和Shu[2]提出的WENO格式(即圖中的WENO-JS)，Borgers[5]提出的WENO-Z格式，以及Sun等[22]提出的線性MDCD格式也被引入作為參考?？梢钥吹?，WENO-Z的耗散小于WENO-JS，但是仍比兩種混合格式要大。混合格式MDCD-HY和MDAD-HY的耗散特性在低波數(shù)和高波數(shù)段與線性MDCD格式接近，但在中波數(shù)段耗散變大。通過局部放大，可以觀察到MDAD-HY在k<1.04時耗散為0，小于MDCD-HY的耗散，之后隨著波數(shù)的增加，耗散逐漸增大，直到k>2.5后MDAD-HY的耗散特性曲線與MDCD-HY基本接近。這就說明MDAD-HY格式在低波數(shù)段相比MDCD-HY有更好的分辨率，而在高波數(shù)段能提供和MDCD-HY相似的耗散來抑制數(shù)值振蕩，具有較好的魯棒性。

(a)

圖14為不同格式的近似色散特性對比，可以看到MDAD-HY和MDCD-HY的色散特性非常接近，在高波數(shù)段稍差于MDCD，但遠好于WENO-JS和WENO-Z。

圖14 不同格式的近似色散特性

4 MDCD格式的發(fā)展應(yīng)用與MDAD格式的算例測試

前面我們以線性波動方程為例介紹了MDCD和MDAD格式，這些算法很容易推廣到Euler方程和N-S方程。具體方案是：采用特征分解[29]的方法計算無黏數(shù)值通量，尺度識別器和間斷探測器都作用于特征分解后的局部特征通量。更多細節(jié)可以參考文獻[22,30]。而黏性通量采用文獻[31]中的方法進行計算。時間積分采用低存儲的四階龍格-庫塔格式[32]，也可以采用常見的隱式格式。

MDCD格式自提出以來，得到了持續(xù)的發(fā)展和廣泛的應(yīng)用，而MDAD格式及其對應(yīng)的混合格式MDCD-HY最近才被提出。因此在本節(jié)中，對于MDCD格式我們主要介紹其發(fā)展和應(yīng)用，而對于MDAD格式，我們主要介紹其在求解一些標準算例中的表現(xiàn)。

4.1 MDCD格式的發(fā)展及應(yīng)用

MDCD格式[22]自提出以來，得到了一定的發(fā)展和應(yīng)用。原始MDCD格式是四階精度的，Sun等[24]研究了六階精度的MDCD格式，并初步研究了耗散控制問題。Wang等[33]將MDCD格式推廣到有限體積格式，通過采用多塊結(jié)構(gòu)網(wǎng)格，使算法具備了模擬復(fù)雜幾何外形的能力。圖15是用有限體積MDCD格式得到的AIAA CFD HiLiftPW-1多段翼模型標準算例的壓力分布云圖。圖16是41%翼展的壓力系數(shù)分布，對比同樣網(wǎng)格上MUSCL格式的計算結(jié)果，可以看到MDCD格式預(yù)測的壓力分布更接近實驗結(jié)果。

(a) 機身壓力分布云圖

圖16 41%翼展壓力系數(shù)分布

文獻[34-35]將MDCD格式發(fā)展到求解多介質(zhì)復(fù)雜流動，并研究了激波與多物質(zhì)界面作用產(chǎn)生的界面不穩(wěn)定性和湍流混合問題。圖17是激波與氦氣泡作用后，界面失穩(wěn)的形態(tài)。

圖17 激波與氦氣泡作用后的界面不穩(wěn)定性

此外，Li等[36]、Chen等[37]利用MDCD格式耗散可調(diào)的特點，研究了MDCD格式的耗散對隱式大渦模擬的影響。Jiang等[38]把MDCD格式與SIMPLE算法相結(jié)合，將MDCD格式推廣到求解不可壓縮流動，顯著提高了算法的分辨率。Duan等[39-40]把MDCD格式應(yīng)用于高超聲速邊界層粗糙元誘發(fā)轉(zhuǎn)捩問題的直接數(shù)值模擬，Jiang等[41]把MDCD格式應(yīng)用于低壓渦輪分離誘導(dǎo)轉(zhuǎn)捩的預(yù)測，Yang等[42]把MDCD格式應(yīng)用于基于渦流發(fā)生器的流動控制研究。這些應(yīng)用充分說明了MDCD格式解決實際科學(xué)與工程問題的能力。

4.2 MDAD格式算例測試

MDAD格式是MDCD格式的進一步發(fā)展，本節(jié)通過典型算例說明MDAD格式具有比采用傳統(tǒng)耗散控制方法的MDCD格式具有更好的性能。由于篇幅的限制，我們只介紹二維無黏和黏性流動兩個算例，更多的測試算例請參見文獻[25]。

4.2.1 雙馬赫反射問題

考慮文獻[43]中提到的雙馬赫反射問題作為二維驗證算例。該算例的控制方程為二維Euler方程。一個馬赫數(shù)為10、與x軸夾角為60°的激波向右移動，并與下壁面發(fā)生反射。初始條件如下:

其中x0=1/6，計算域為(x,y)∈[0,4]×[0,1]。邊界條件方面，在下邊界，為了使激波附著在壁面前緣，x∈[0,x0]段采用激波后的參數(shù)，而x∈[x0,4]段為滑移壁面。入口和出口分別采用來流和出流邊界條件。上邊界采用精確的波前波后參數(shù)，并隨著激波的移動而變化。

圖18展現(xiàn)了計算時間t=0.18時不同格式的密度等值線，網(wǎng)絡(luò)量采用800×200，庫朗數(shù)Nc=0.3?？梢钥吹?，WENO-JS、WENO-Z、MDCD-HY和MDAD-HY均能很好地捕捉到流場結(jié)構(gòu)，包括馬赫桿和壁面射流。但是，MDAD-HY在第二個三波點附近的渦結(jié)構(gòu)明顯更加細致，對小尺度結(jié)構(gòu)捕捉更好。這說明相比于其他三種格式，MDAD-HY格式有著更高的分辨率。

(a) WENO-JS

4.2.2 弱激波與等熵渦相互作用問題

為了研究不同格式在二維N-S方程中的效果，可以考察帶有黏性的弱激波與等熵渦的相互作用問題[44]。該算例的計算域為[0,2L0]×[0,2L0]，其中參考長度L0=1.0。一個靜止的平面激波位于x=1.0，波前來流參數(shù)為ρ0=1.0，p0=1/γ，預(yù)設(shè)的激波前后壓差為Δp/p0=0.4，對應(yīng)的激波強度為M0=1.1588，其余的流場參數(shù)由Rankine-Hugoniot關(guān)系確定。雷諾數(shù)Re=5000，它的參考長度為L0，參考密度和參考速度均為來流參數(shù)。初始狀態(tài)下，一個孤立的等熵渦被疊加在來流中，渦中心位于(x0,y0)=(0.5,1.0)，它表現(xiàn)為速度(u,v)和溫度T=γp/ρ的擾動，而熵S=p/ργ為常數(shù)。

(a) t=0.7時中心線上的密度分布

5 結(jié) 論

本文綜述了我們在有限差分格式的色散優(yōu)化和耗散控制等方面的工作。相應(yīng)的計算格式可以計算無激波的光滑解，也可以作為非線性激波捕捉格式(如WENO格式)或者混合格式的線性部分。

在理論上，我們得到了半離散格式色散、耗散獨立的充分條件，使色散優(yōu)化和耗散控制不相互干擾，為格式性能優(yōu)化打下了基礎(chǔ)。進一步的分析表明，即使對于全離散格式，上述充分條件也可以保證色散和耗散近似獨立。

以此為基礎(chǔ)，我們得到了色散最小、耗散可控的MDCD格式，該格式的優(yōu)良性能已經(jīng)在一系列應(yīng)用中得到了驗證，并被推廣到有限體積方法，以及求解多介質(zhì)流動、不可壓流動等。MDCD格式的主要缺點在于耗散控制。無論是采用試錯方法還是基于色散、耗散條件確定通用耗散參數(shù)，得到的格式在耗散特性方面仍有很大的改進空間。

為了解決耗散控制問題，我們提出了將數(shù)值耗散與數(shù)值解的尺度相關(guān)聯(lián)的思路。為此，我們首先提出了一種基于數(shù)值通量模板點上導(dǎo)數(shù)的尺度識別器，用于判斷流場局部的等效無量綱波數(shù)。然后通過色散耗散條件[23]，設(shè)計出耗散參數(shù)與等效無量綱波數(shù)之間的關(guān)系，從而實現(xiàn)耗散參數(shù)的自適應(yīng)調(diào)節(jié)，提出了色散最小、耗散自適應(yīng)的MDAD格式。

我們通過所發(fā)展的格式與WENO格式相混合的方法，使MDCD和MDAD格式在能夠高精度、高分辨率的求解光滑流動的同時，具有魯棒的激波捕捉能力。我們發(fā)現(xiàn)，通過在混合格式的激波探測器中，引入基于尺度識別器的修正，可以進一步提高混合格式的性能。

由于MDCD格式的性能已經(jīng)得到了充分的驗證，且MDAD格式在多個標準算例中都表現(xiàn)出優(yōu)于MDCD格式的性能，可以預(yù)見，MDAD格式及其對應(yīng)的混合格式MDAD-HY在求解含激波的復(fù)雜多尺度流動中也有很好的應(yīng)用前景。

致謝:李萬愛、潘建華、曾衛(wèi)剛、黃文鋒、李星輝、張宏杰、李俊濤、張潤東等參與了計算程序開發(fā)、驗證、應(yīng)用等工作，特致感謝。