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

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

0 引 言

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

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

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

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

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

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

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

為了使格式具有激波捕捉能力，通過將MDAD格式的數(shù)值通量表示為一系列候選模板上的低階多項(xiàng)式的線性組合，并引入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)造混合格式的一個(gè)關(guān)鍵問題是合理的間斷探測(cè)器。為此，Li等[25]將傳統(tǒng)間斷探測(cè)器[13]與尺度識(shí)別器相結(jié)合，提出了一種新的間斷探測(cè)器，能夠有效減少對(duì)光滑區(qū)的極點(diǎn)的誤判。

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

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

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

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

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

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

1.1 MDCD格式的半離散特性

當(dāng)空間有限差分格式包含對(duì)稱的(2r+1)個(gè)模板點(diǎn)時(shí)，式(3)可寫為:

當(dāng)r=3時(shí)，四階的MDCD格式[22]形式如下：

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

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

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

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

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

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

(11)

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

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

(12)

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

γdisp=0.046 378 3

(13)

圖1 不同格式的色散特性

圖2 不同格式的色散誤差

對(duì)于耗散特性，為了使格式保持穩(wěn)定，所有波數(shù)下的耗散都應(yīng)非負(fù)。由式(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。因此當(dāng)a>0時(shí)，如果γdiss≥0，則MDCD格式有非負(fù)耗散，格式是穩(wěn)定的，且隨著γdiss的增大，格式的耗散增加。當(dāng)耗散參數(shù)取0時(shí)，MDCD為中心格式，此時(shí)格式分辨率最高，但是無法抑制高波數(shù)的非物理振蕩。

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

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

γdiss=0.012

(17)

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

圖3 不同格式的耗散特性

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

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

考慮線性對(duì)流方程：

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

其中k′為1.1節(jié)中所述的無量綱修正波數(shù)。對(duì)比式(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為半離散色散誤差。若采用顯式龍格庫(kù)塔格式進(jìn)行時(shí)間積分，則在t=Δt時(shí)刻有：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

第二步：校準(zhǔn)。在這一步中，采用正弦波來確定式(28)中的參數(shù)。若將式(28)定義的長(zhǎng)度尺度校準(zhǔn)為精確的正弦波無量綱波長(zhǎng)，最簡(jiǎn)單的方式是令p=q+2，此時(shí)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)

因此，當(dāng)β=2π時(shí)，無量綱尺度δ即為精確的無量綱波長(zhǎng)，

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

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

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

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

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

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

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

(a)

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

2.2 MDAD格式

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

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

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

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

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

3 混合格式

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

3.1 MDAD-WENO格式

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

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

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

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

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

代入每個(gè)候選模板得到:

(48)

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

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

3.2 MDAD-HY格式

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

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

其中，

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

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

(56)

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

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

(a) Harten

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

3.3 MDAD-HY格式的譜特性

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

圖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時(shí)耗散為0，小于MDCD-HY的耗散，之后隨著波數(shù)的增加，耗散逐漸增大，直到k>2.5后MDAD-HY的耗散特性曲線與MDCD-HY基本接近。這就說明MDAD-HY格式在低波數(shù)段相比MDCD-HY有更好的分辨率，而在高波數(shù)段能提供和MDCD-HY相似的耗散來抑制數(shù)值振蕩，具有較好的魯棒性。

(a)

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

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

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

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

MDCD格式自提出以來，得到了持續(xù)的發(fā)展和廣泛的應(yīng)用，而MDAD格式及其對(duì)應(yīng)的混合格式MDCD-HY最近才被提出。因此在本節(jié)中，對(duì)于MDCD格式我們主要介紹其發(fā)展和應(yīng)用，而對(duì)于MDAD格式，我們主要介紹其在求解一些標(biāo)準(zhǔn)算例中的表現(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多段翼模型標(biāo)準(zhǔn)算例的壓力分布云圖。圖16是41%翼展的壓力系數(shù)分布，對(duì)比同樣網(wǎng)格上MUSCL格式的計(jì)算結(jié)果，可以看到MDCD格式預(yù)測(cè)的壓力分布更接近實(shí)驗(yàn)結(jié)果。

(a) 機(jī)身壓力分布云圖

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

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

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

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

4.2 MDAD格式算例測(cè)試

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

4.2.1 雙馬赫反射問題

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

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

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

(a) WENO-JS

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

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

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

5 結(jié) 論

本文綜述了我們?cè)谟邢薏罘指袷降纳?yōu)化和耗散控制等方面的工作。相應(yīng)的計(jì)算格式可以計(jì)算無激波的光滑解，也可以作為非線性激波捕捉格式(如WENO格式)或者混合格式的線性部分。

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

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

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

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

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

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