王來，吳頌平

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191)

氣動聲學(xué)以及電磁波傳播的數(shù)值模擬在過去幾十年里大量采用高階數(shù)值格式.如今在湍流的大渦模擬以及直接數(shù)值模擬等數(shù)值模擬中，高階數(shù)值格式也得到了廣泛的應(yīng)用.湍流問題中存在大范圍的時間/空間尺度，低階格式由于耗散過大，往往無法捕捉到流動中不穩(wěn)定的、微小的尺度.高階格式精度更高，有著很好的應(yīng)用前景.然而，高階格式在求解超音速可壓縮流動問題時，僅有高精度還不能滿足實際需求.為了準確地捕捉流動中的間斷，還需要做到基本無數(shù)值振蕩，并且具有較高的間斷分辨率.總而言之，精度高、分辨率高、捕捉間斷能力強、魯棒性好的數(shù)值格式，是學(xué)者們不斷追求的目標［1-2］.

有限差分格式可以分為顯式格式和隱式格式.其中高階的顯式格式以ENO格式和WENO格式［3-4］為代表.這類格式能夠?qū)崿F(xiàn)高精度、高分辨率，并且具有良好的捕捉激波的能力.但是，在湍流的數(shù)值模擬中，ENO和WENO格式往往會表現(xiàn)出過大的耗散.WENO-Z［5-7］格式通過對傳統(tǒng)的WENO-JS格式中權(quán)重算子的改進，一定程度上減小了耗散，提高了分辨率.

緊致格式是隱式格式.Lele［8］最早系統(tǒng)地給出了具有類譜方法分辨率的非守恒型中心型緊致格式，能夠在較小的模板上輕松地獲得很高的精度.但是，由于這類緊致格式的耗散誤差為零，最早只被應(yīng)用于不可壓縮流動問題的求解.Liu等［9］對守恒型中心緊致格式進行了系統(tǒng)的介紹.爾后，學(xué)者們通過各種手段將緊致格式推廣到可壓縮流動問題的模擬.其中，Adams等［10］將 ENO格式與緊致格式結(jié)合，以實現(xiàn)捕捉間斷的能力.Pirozzoli［11］依照同樣的思路，將 WENO 格式與緊致格式結(jié)合起來.然而這些混合格式中，兩種子格式間的相互切換顯得有些“突然”，往往會在格式切換的附近造成數(shù)值振蕩.為此，Ren等［12］提出了新的權(quán)重計算方法，其中引入了自由參數(shù)rc.國內(nèi)學(xué)者也提出了不同的引入自由參數(shù)的混合格式權(quán)重算法［13-14］.它們都表現(xiàn)出了優(yōu)異的激波捕捉能力以及分辨率特性.但是，這類自由參數(shù)往往需要根據(jù)經(jīng)驗進行選取，需要進行大量的嘗試，這就造成了使用上的不便.

本文的目標就是要構(gòu)造出一種無自由參數(shù)的權(quán)重算法［15］，在保持混合格式的高精度、高分辨率特性的同時，提高格式的易用性和魯棒性.

1 混合格式的構(gòu)造

1.1 控制方程

考慮雙曲型標量守恒律方程:

將計算區(qū)域等距離劃分 xj=jΔx，j=0，1，…，N.守恒形式的半離散有限差分格式為

為了提高數(shù)值格式的魯棒性，通常將數(shù)值通量分裂為正通量和負通量2個部分:

由于這2個部分的對稱性，本文只討論正通量的計算過程.

1.2 守恒型緊致格式

本文構(gòu)造了2種守恒型緊致格式［9，11］

這2種緊致格式的截斷誤差分別為

截斷誤差的第1和第2兩項分別為格式中耗散誤差以及色散誤差的主項.不難看出，第2種緊致格式的耗散誤差更小.

1.3 WENO-Z 格式

本文的混合格式采用 WENO-Z格式［5-7］.相比WENO-JS格式［4］，WENO-Z格式在極值點附近能夠更好地保持高精度.為了論文的完整性，這里給出五階WENO-Z的基本形式.

式中，qr3為候選模板的數(shù)值通量，詳見文獻［4-5］;ωr(r=0，1，2)是模板 r的權(quán)重:

式中，cr為理想權(quán)重;ε為一個極小數(shù);βr為候選模板的光滑因子:

1.4 兩種子格式的混合

混合格式權(quán)重的理想狀態(tài)是:在間斷區(qū)域WENO-Z格式的權(quán)重為1，從而提高格式對間斷的捕捉能力;在光滑區(qū)域緊致格式的權(quán)重為1，這樣一來就能夠保持混合格式在光滑區(qū)域的低耗散特性.本文將按照這一標準構(gòu)造新的權(quán)重算子.

混合格式的模板為［xj－2，xj－1，xj，xj+1，xj+2］，對光滑因子 β0，β1，β2進行 Taylor展開分析:

本文為混合格式中緊致格式的權(quán)重設(shè)計了新算子:

式中ε是一個極小的數(shù)，可以取計算機所能存儲的最小浮點數(shù)，因此這并不能算作是自由參數(shù).在光滑區(qū)域:

也就是說，在光滑區(qū)域，緊致格式的權(quán)重幾乎就是1.當(dāng)混合格式的總模板存在間斷時，τ遠大于，如此一來，緊致格式的權(quán)重就接近0.本文將這種新型算子稱為Z型權(quán)重算子.

Pirozzoli在文獻［11］中定義了一個簡單的權(quán)重算子，這種算子造成了混合格式中兩種子格式的切換過于突然，切換點會產(chǎn)生不容忽視的數(shù)值振蕩乃至污染流場.為了克服上述問題Ren［12］和武從海等［14］提出了一個新的權(quán)重算子:

這種權(quán)重算子表現(xiàn)出了良好的分辨率以及對激波的精確捕捉能力.但是，該算子較為復(fù)雜，格式的表現(xiàn)很大程度上依賴算子中自由參數(shù)rc的選取，rc過小會造成計算無法進行.本文將這類算子稱作R型算子.

新權(quán)重算子Z避免了自由參數(shù)的引入以及邏輯判斷的使用［15］.為了測試新算子的特性，主要將新算子的計算結(jié)果與R型算子的結(jié)果進行比較.本文用HCW-U表示式(5)型緊致格式與WENO-Z的混合格式，用HCW-UL表示式(6)型緊致格式與WENO-Z的混合格式.HCW-UL-Z型格式即為采用Z型權(quán)重的式(6)型迎風(fēng)型低耗散緊致格式與WENO-Z格式的混合格式，文中提到的其他簡稱構(gòu)成與此類似.

2 數(shù)值試驗結(jié)果

本文所涉及到的一、二維算例控制方程均為歐拉方程，時間離散采用三階龍格庫塔(RK)方法，這里不再贅述，詳見文獻［4］.

2.1 激波管問題

LAX問題是1D激波管問題的典型算例之一.在區(qū)間［－5.0，5.0］之內(nèi)，以原點為分界點，左右兩側(cè)的氣體初始狀態(tài)不同.初始條件如下:

計算網(wǎng)格為200，計算終止時刻為t=1.3，CFL=0.3.本算例中R型權(quán)重算子取值為rc=0.5.

由圖1、圖2不難發(fā)現(xiàn)HCW-U以及HCW-UL采用Z型權(quán)重算子時，對激波以及接觸間斷的捕捉與采用R型權(quán)重時差不多.值得注意的是，這兩種混合格式在采用R型權(quán)重算子時在膨脹波頭造成的數(shù)值振蕩都比采用Z型權(quán)重時嚴重，詳見圖2中的局部放大圖.

圖1 HCW-U，HCW-UL，LAX 問題，密度值Fig.1 HCW-U，HCW-UL，LAX problem，density

圖2 HCW-U，HCW-UL，LAX 問題，壓強值Fig.2 HCW-U，HCW-UL，LAX problem，pressure

2.2 激波與熵波的干涉

本算例描述的是激波與熵波的干涉問題，其數(shù)值結(jié)果包含了間斷以及不同尺度的波.該算例能夠很好地測試數(shù)值格式對間斷的捕捉能力以及分辨率特性.計算區(qū)域為［0，10］，初始條件如下:

計算網(wǎng)格為301，計算終止時刻為t=1.8，CFL=0.1.混合格式采用 R 型權(quán)重時，rc=0.5.該問題沒有精確解，但是可以用五階WENO格式在網(wǎng)格數(shù)為1 601下的數(shù)值計算結(jié)果作為近似的精確解.為了便于觀察，本文給出了密度在區(qū)間［5，7.4］內(nèi)的數(shù)值結(jié)果，這一區(qū)間內(nèi)密度的數(shù)值結(jié)果能夠很好地說明格式的分辨率、耗散性質(zhì)以及對間斷的捕捉能力.

從圖3、圖4可以清楚地看到，對于HCW-U以及HCW-UL格式，Z型權(quán)重與R型權(quán)重都表現(xiàn)出了對間斷的良好捕捉能力以及對不同尺度波的分辨能力，混合格式相比于WENO-Z格式耗散大大降低，分辨率有了明顯的提高.對比圖3、圖4，HCW-U以及HCW-UL計算效果并沒有太顯著的區(qū)別，盡管后者在理論上來說耗散比前者更小.

圖3 HCW-U，Osher-Shu問題，區(qū)間［5，7.4］密度值Fig.3 HCW-U，Osher-Shu problem，density in［5，7.4］

圖4 HCW-UL，Osher-Shu問題，區(qū)間［5，7.4］密度值Fig.4 HCW-UL，Osher-Shu problem，density in［5，7.4］

2.3 2D雙馬赫反射

強激波的雙馬赫反射(DMR)問題是測試數(shù)值格式的分辨率以及對間斷的捕捉特性的標準算例之一.該問題的計算區(qū)域為［0，4］×［0，1］，初始條件為:馬赫數(shù)為10、與x軸成60°斜激波在x=1/6處與底部邊界相遇，激波上游(ρ，u，v，p)=(1.4，0，0，1)，下游參數(shù)滿足 RH 關(guān)系式.計算域的上邊界條件是激波傳播的精確解，左邊界以及底部段為入流邊界條件，右邊界為出流邊界條件，底部段為壁面邊界.計算終止時間t=0.2，CFL=0.5.計算網(wǎng)格為 801 ×201.

混合格式采用R型權(quán)重時，rc=0.5.由于主要的流場信息都在［0，3］×［0，1］區(qū)間之內(nèi)，本文僅僅給出了該區(qū)域內(nèi)的密度云圖，所有云圖的等值線均為將區(qū)間2～22均分為50等份.對比圖5、圖6，可以清楚地看到，HCW-U-Z與 HCWU-R都很好地分辨出了滑移線的卷曲.HCW-ULZ的計算結(jié)果由圖7給出，HCW-UL-R的計算結(jié)果略去.值得注意的是，采用R型權(quán)重計算時，如果rc過小(以rc=0.3為例)，計算無法進行.Z型權(quán)重由于未引入自由參數(shù)，魯棒性提高，在針對陌生問題進行數(shù)值求解時，則可以避免自由參數(shù)取得不當(dāng)造成計算無法進行的問題，與此同時，保持了很好的分辨率特性.

圖5 HCW-U-R，DMR問題，密度云圖，50條等值線Fig.5 HCW-U-R，DMR problem，density contour，50 levels

圖6 HCW-U-Z，DMR問題，密度云圖，50條等值線Fig.6 HCW-U-Z，DMR problem，density contour，50 levels

圖7 HCW-UL-Z，DMR問題，密度云圖，50條等值線Fig.7 HCW-UL-Z，DMR problem，density contour，50 levels

3 結(jié)論

1)新權(quán)重算子(Z型)對間斷的捕捉能力良好，同時相比R型權(quán)重能夠抑制間斷處數(shù)值振蕩的傳播.

2)相比于WENO-Z，采用Z型權(quán)重算子的混合格式HCW-Z耗散降低，分辨率提高.HCW-U與HCW-UL的數(shù)值結(jié)果并無明顯區(qū)別.

3)Z型權(quán)重算子的數(shù)值特性與R型權(quán)重算子(rc=0.5)相比區(qū)別不明顯，但是R型權(quán)重算子的計算效果依賴rc的選取，使用受到局限.

References)

［1］ Wang Z J，F(xiàn)idkowski K，Abgrall R，et al.High-order CFD methods:current status and perspective［J］.Interntional Journal for Numerical Methods in Fluids，2012，72(8):811-845.

［2］ 閻超，于劍，徐晶磊，等.CFD模擬方法的發(fā)展成就與展望［J］.力學(xué)進展，2011，41(5):562-589.Yan C，Yu J，Xu J L，et al.On the achievements and prospects for the methods of the computational fluid dynamics［J］.Advances in Mechanics，2011，41(5):562-589(in Chinese).

［3］ Harten A，Osher S.Uniformaly high order accurate essentially nonoscillatory scheme，Ⅲ［J］.Journal of Computational Physics，1987，71(2):231-303.

［4］ Shu C W，Osher S.Efficient implementation of essentially nonoscillatory shock-capturing schemes［J］.Journal of Computational Physics，1988，77(2):439-471.

［5］ Borges R，Carmona M，Costa B，et al.An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws［J］.Journal of Computational Physics，2008，227(6):3191-3211.

［6］ Henrick A K，Aslam T D，Powers J M.Mapped weighted essentially non-oscillatory schemes:achieving optimal order near critical points［J］.Journal of Computational Physics，2005，207(2):542-567.

［7］ Castro M，Costa B，Don W S.High order weighted essentially nonoscillatory WENO-Z schemes for hyperbolic conservation laws［J］.Journal of Computational Physics，2011，230(5):1766-1792

［8］ Lele S K.Compact finite difference schemes with spectral-like resolution［J］.Journal of Computational Physics，1992，103(1):16-42.

［9］ Liu X L，Zhang S H，Zhang H X，et al.A new class of central compact schemes with spectral-like resolutionⅠ:linear schemes［J］.Journal of Computational Physics，248(1):235-256.

［10］ Adams N A，Shariff K.A high-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems［J］.Journal of Computational Physics，1996，127(1):27-51.

［11］ Pirozzoli S.Conservative hybrid compact-WENO schemes for shock-turbulence interaction［J］.Journal of Computational Physics，2002，178(1):81-117.

［12］ Ren Y X，Liu M E，Zhang H X.A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws［J］.Journal of Computational Physics，2003，192(2):365-386.

［13］ Yu J，Yan C，Jiang Z H.A high resolution low dissipation hybrid scheme for compressible flows［J］.Chinese Journal of Aeronautics，2011，24(4):417-424.

［14］ 武從海，趙寧，田琳琳.一種改進的緊致 WENO混合格式［J］.空氣動力學(xué)報，2013，31(4):477-481.Wu C H，Zhao N，Tian L L.An improved hybrid cmpact-WENO scheme［J］.Acta Aero Dynamica Sinica，2013，31(4):417-481(in Chinese).

［15］ Shen Y Q，Zha G C.Generalized finite compact difference scheme for shock/complex flowfield interaction［J］.Journal of Computational Physics，2011，230(12):4419-4436.