武從海，趙 寧，田琳琳

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016）

0 引 言

自1992年Lele提出了一類緊致差分格式的構(gòu)造方法后［1］，緊致格式受到研究者的廣泛關(guān)注［2－4］。緊致格式不僅能在較少的網(wǎng)格點(diǎn)上取得較高精度，而且分辨率接近譜方法，同時(shí)對計(jì)算域和邊界的要求不像譜方法那樣苛刻。但是，直接應(yīng)用緊致格式計(jì)算含有間斷的流場時(shí)會(huì)產(chǎn)生非物理振蕩。為了消除緊致格式在間斷附近的非物理振蕩，必須采用非線性的處理。

Ravichandran［5］對緊致格式進(jìn)行通量限制得到了具有TVD性質(zhì)的高階緊致格式。鄧小剛［6］等采用加權(quán)技術(shù)構(gòu)造了一系列非線性緊致格式（加權(quán)緊致非線性格式）。Wang等［7］在Pade格式中加入ENO算法，得到一種新的高階ENO－Pade格式，此格式能消除間斷附近的非物理振蕩。另外一類處理方法是將緊致格式與其他捕捉間斷格式混合，這類方法兼顧了高分辨率，高計(jì)算效率以及良好的魯棒性。Adams等［8］首先構(gòu)造了一種高階緊致格式和ENO格式的混合數(shù)值格式，在流場光滑的區(qū)域采用緊致格式，在間斷附近采用ENO格式。Pirozzoli［9］構(gòu)造了一類守恒型的緊致格式與WENO格式混合，由于采用了守恒型的格式以及穩(wěn)定性比ENO更好的WENO，其結(jié)果比Adams的結(jié)果更好。Ren等［10］將Pirozzoli的格式進(jìn)一步改進(jìn)，在格式的混合處理上采用了加權(quán)平均，避免了兩種不同格式之間的突然切換；另外，對于Euler方程采用了特征型的格式，使得對復(fù)雜流動(dòng)的分辨率有了明顯的提高，但該方法需要求解塊三對角線性方程組，相比僅需求解三對角線性方程組的Pirozzoli的混合格式，增加了許多計(jì)算量。

本文基于Ren的混合格式，使用了間斷分辨率更好的WENO子格式［7］。另外，為加強(qiáng)格式的穩(wěn)定性，本文改進(jìn)了子格式權(quán)系數(shù)的計(jì)算公式。對于Euler方程，本文采用全局分裂，且僅對 WENO子格式進(jìn)行特征投影處理，無需求解塊三對角方程組，相比Ren的混合格式，計(jì)算量大大減少，且依然得到了相近的結(jié)果。

1 混合緊致WENO格式

對于一維標(biāo)量守恒律方程：

取均勻網(wǎng)格劃分，xi＝iΔx，i＝0，1，…，N，Ii＝［xi－1／2，xi＋1／2］，其中單元邊界xi±1／2＝xi±Δx／2。式（1）的半離散化守恒差分格式為：

采用全局的Lax－Friedrichs分裂方法，f可分解為f＋與f－之和。假設(shè)f′＞0，對于f′＜0的情況，可由f′＞0時(shí)所用的格式作對稱變換得到。

本文采用的緊致格式為Pirozzoli提出的5階迎風(fēng)緊致格式的歸一化形式［5］，即

捕捉間斷格式采用改進(jìn)的5階WENO格式，計(jì)算公式為：

權(quán)系數(shù)ωk的計(jì)算過程請參見文獻(xiàn)［11］。

本文使用加權(quán)平均的方法混合兩個(gè)子格式，式（3）和式（4）的權(quán)系數(shù)分別為σi＋1／2和（1－σi＋1／2）。當(dāng)權(quán)系數(shù)σi＋1／2分別取1和0時(shí)，格式也分別化簡為5階迎風(fēng)緊致格式和5階WENO格式。根據(jù)建立混合格式的原則，權(quán)系數(shù)σi＋1／2由流函數(shù)的光滑程度確定。使用ri＋1／2來度量xi＋1／2附近的流場光滑程度，

Ren等［10］給出的權(quán)系數(shù)的計(jì)算公式為σi＋1／2＝其中閾值rc的范圍在區(qū)間（0，1）內(nèi)。這種權(quán)值的選取方式可以使計(jì)算量大的WENO格式僅在間斷附近啟用，從而提高計(jì)算效率。本文采用權(quán)系數(shù)

其中p為大于1的整數(shù)，文獻(xiàn)［6］中的權(quán)系數(shù)對應(yīng)p＝1的情況。當(dāng)p增大時(shí)，在間斷附近權(quán)系數(shù)會(huì)變小，WENO格式所占權(quán)重會(huì)變大，那么格式的穩(wěn)定性會(huì)有所增強(qiáng)。以下一節(jié)中的激波熵波相互作用問題為例，當(dāng)閾值rc過小時(shí)，計(jì)算過程會(huì)出現(xiàn)負(fù)壓強(qiáng)。若保證不出現(xiàn)負(fù)壓強(qiáng)，則不同p值對應(yīng)相應(yīng)的rc的下限。表1列出了一些p值及對應(yīng)的rc的取值下限。由該表可知，rc的取值下限在p＝1時(shí)比較大，在p取其他值時(shí)差別較小，并且隨著p增大而緩慢減小。當(dāng)p＝100時(shí)，rc的取值下限為0.33，反而比p＝2、3、4、5時(shí)小。實(shí)際上，當(dāng)p足夠大時(shí)，權(quán)系數(shù)的取值幾乎全部在0和1附近，這樣導(dǎo)致了間斷附近兩個(gè)子格式之間的突然切換，降低了穩(wěn)定性。那么，p取2、3、4、5時(shí)比p＝1時(shí)格式的穩(wěn)定性好一些。另外，在實(shí)際計(jì)算中，當(dāng)p分別取2、3、4、5時(shí)，計(jì)算結(jié)果之間的差別很小，因此本文選取p＝2。

表1 激波熵波問題中，N＝200，不同p值對應(yīng)的rc的取值下限Table 1 Lower limit of rcfor different p，Shu－Osher problem，N＝200

在Euler方程的計(jì)算中，若有間斷存在，使用特征投影處理可以使計(jì)算結(jié)果更加準(zhǔn)確；若流場較為光滑，沒有必要使用特征投影處理?；旌细袷皆诹鲌龅墓饣瑓^(qū)域中采用了緊致格式，對緊致格式特征投影處理會(huì)產(chǎn)生塊三對角方程組，而求解塊三對角方程組會(huì)引起計(jì)算量的大大增加。本文僅對混合格式中的WENO格式進(jìn)行特征投影處理。對于一維Euler方程＝0，u為守恒變量組成的向量，f為通量。本文格式的具體計(jì)算過程為：

（1）令λj＊＝其中λj為特征速度，xi取所有的節(jié)點(diǎn)，j＝1、2、3；

（2）令α＝1、2、3，對于每個(gè)節(jié)點(diǎn)進(jìn)行通量分裂

（3）根據(jù)密度場的光滑程度，由式（5）、（6）計(jì)算兩種格式的權(quán)系數(shù)，對于需要WENO格式計(jì)算的區(qū)域，分別對f＋、f－進(jìn)行特征投影后進(jìn)行 WENO插值，然后逆向投影；

（4）在連續(xù)區(qū)域采用迎風(fēng)緊致格式，在間斷區(qū)域由兩種格式加權(quán)得到混合格式，結(jié)合時(shí)間格式即為全部數(shù)值求解過程。

除了WENO權(quán)因子不同外，本文格式中的WENO格式與文獻(xiàn)［12］中 WENO格式處理有所不同。由于本文首先進(jìn)行全局通量分裂，所以WENO計(jì)算過程是對分裂后的通量進(jìn)行的。由于組元型WENO格式在接觸間斷處會(huì)出現(xiàn)振蕩，Pirozzoli的混合格式在Lax問題中存在此問題。因此，本文格式WENO部分對全局分裂后的通量進(jìn)行了特征投影處理。盡管這樣的處理會(huì)稍稍加大計(jì)算量，但是混合格式里需要WENO計(jì)算的部分所占比例一般較小。例如，對于雙馬赫反射問題，當(dāng)rc＝0.4時(shí)，Ren的混合格式需要WENO計(jì)算的比例根據(jù)時(shí)間不同為0.9%到7.4%不等，所以WENO格式的處理對總計(jì)算量影響較小。在Euler方程的數(shù)值試驗(yàn)中，該混合格式的計(jì)算量甚至低于Pirozzoli的低耗散LF分裂的混合格式，比Ren的混合格式計(jì)算量小很多。

2 數(shù)值算例與分析

本文給出了一種用于激波捕捉計(jì)算的守恒型混合緊致WENO格式。本節(jié)討論該混合格式在求解一維和二維Euler方程的應(yīng)用，并將計(jì)算結(jié)果與相關(guān)的幾種格式進(jìn)行比較。這里時(shí)間格式均采用3階TVD龍格庫塔格式［8］。

2.1 激波熵波相互作用問題

激波熵波相互作用問題是一個(gè)一維Euler方程問題，其初值條件為：

圖1 幾種格式在激波熵波問題的計(jì)算結(jié)果的比較，t＝1.8Fig.1 Comparison of different schemes，Shu－Osher problem，t＝1.8

這里取rc＝0.35，網(wǎng)格點(diǎn)數(shù)為200，CFL數(shù)取0.5。該問題中Mach數(shù)為3的激波與正弦波相互干擾，誘發(fā)高頻振動(dòng)。同時(shí)，采用Pirozzoli的低耗散LF分裂的混合格式（β＝0.3）和Ren的混合格式（rc＝0.35）對該算例進(jìn)行了計(jì)算，并與本文的提出的格式 比較。圖1（a）為t＝1.8時(shí)的密度圖，圖1（b）為圖1（a）的局部放大圖。圖中實(shí)線為5階WENO格式采用4000個(gè)網(wǎng)格點(diǎn)計(jì)算的結(jié)果，作為近似的準(zhǔn)確解。從圖中可以看出，在網(wǎng)格相對較粗的情況下，三種格式均能較好地捕捉流場的細(xì)節(jié)。其中Ren的格式表現(xiàn)最佳，與“準(zhǔn)確解”最為接近；Pirozzoli的格式在高頻振動(dòng)區(qū)域耗散最大；本文所用格式與Ren的格式的結(jié)果相近，在高頻振動(dòng)區(qū)域耗散小，大大好于Pirozzoli的格式的結(jié)果。表2列出了本次計(jì)算幾種格式完成該算例在不同網(wǎng)格點(diǎn)數(shù)時(shí)所花費(fèi)的時(shí)間。從時(shí)間效率來看，本文格式比Pirozzoli的格式還要高，大大優(yōu)于Ren的格式。

表2 幾種格式在計(jì)算激波熵波問題時(shí)花費(fèi)的時(shí)間Table2 Elapsed time for different schemes，Shu－Osher problem

2.2 雙馬赫反射問題

該問題描述了一個(gè)向右運(yùn)動(dòng)的Mach數(shù)為10的斜激波，與反射壁面相互作用并產(chǎn)生復(fù)雜的馬赫結(jié)構(gòu)。激波與反射壁面的夾角為60°，在波前未受影響區(qū)域激波密度為1.4，壓強(qiáng)為1［13］。求解區(qū)域以壁面為底部，壁面起始位置為x＝1／6。這里計(jì)算區(qū)域?yàn)椋?，4］×［0，1］。該問題包含強(qiáng)激波和馬赫反射，對格式的精度、分辨率、特別是穩(wěn)定性等都有很高的要求。

這里取rc＝0.5，CFL數(shù)為0.5，計(jì)算時(shí)間取t＝0.2，網(wǎng)格為800×200。在該算例的計(jì)算中，本文格式需要開啟WENO計(jì)算的部分所占比例隨計(jì)算時(shí)間不同為從0.50%到11.39%不等。圖2為本文格式、Ren的混合格式（rc＝0.5）以及5階 WENO格式的計(jì)算結(jié)果在區(qū)域［0，3］×［0，1］的密度等值線圖，圖3和圖4分別為圖2（a）和圖2（b）的馬赫桿附近的放大圖。從圖中可以看出，相比5階WENO格式，兩種混合格式能更好地模擬流場的細(xì)節(jié)，尤其是馬赫桿附近的區(qū)域。與Ren的混合格式相比，本文格式在壁面前端上方的稀疏波區(qū)域的等值線更平滑，數(shù)值噪聲更小。這是由于本文格式采用的全局通量分裂，其耗散性要略大，并且因此無需加入Ren的格式中用來處理紅斑現(xiàn)象的H－correction過程［6］。表3列出了幾種格式在本次計(jì)算中所花費(fèi)的時(shí)間。從該表中可以看到，本文格式和Ren的格式所用時(shí)間分別為1313.281s和13190.359s，本文格式所需時(shí)間僅為Ren的格式的1／10，在時(shí)間效率上表現(xiàn)出了極大的優(yōu)勢。

圖2 雙馬赫反射問題密度等值線圖，計(jì)算網(wǎng)格800×200，50條密度等值線Fig.2 Double Mach reflection problem，grid：800×200，50contours，density

圖3 圖2（a）中馬赫桿附近的放大圖Fig.3 The enlarged portion of Fig.2 （a）near the Mach stems

圖4 圖2（b）中馬赫桿附近的放大圖Fig.4 The enlarged portion of Fig.2 （b）near the Mach stems

表3 幾種格式在雙馬赫反射問題中花費(fèi)的時(shí)間Table3 Elapsed time for different schemes，Shu－Osher problem

3 結(jié) 論

本文構(gòu)造了一種混合緊致WENO格式，其中的子格式WENO采用了Borges改進(jìn)的WENO權(quán)因子［11］，增強(qiáng)了對間斷的分辨能力。另外，本文改進(jìn)了Ren的混合格式中子格式權(quán)系數(shù)的計(jì)算公式，加強(qiáng)了格式的穩(wěn)定性［10］。對于Euler方程，本文首先進(jìn)行全局通量分裂，然后對分裂后的通量進(jìn)行計(jì)算。特征投影僅對格式中的WENO部分進(jìn)行，避免了Ren的混合格式中對緊致格式特征投影引起的塊三對角方程組求解，大大降低了計(jì)算量。文中計(jì)算了兩個(gè)Euler方程的算例，與Ren的混合格式相比，本文格式數(shù)值粘性稍大，但是時(shí)間效率大大提高。

［1］LELE S K.Compact finite difference schemes with spectral－like resolution［J］.J.Comput.Phys，1992，103：16－42.

［2］沈孟育，蔣莉.滿足熵增原則的高精度高分辨率格式［J］.清華大學(xué)學(xué)報(bào)（自然科學(xué)版），1998，16：56－62.

［3］SHANG J S.High－order compact－differenced schemes for time－dependent Maxwell equations［J］.J.Comput.Phys，1999，153：312－333.

［4］傅德薰，馬延文.高精度差分格式及多尺度流場特性的數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1998，16（1）：24－34.

［5］RAVICHANDRAN K S.Higher order KFVS algorithms using compact upwind differenced operators［J］.J.Comput.Phys，1997，130：161－173.

［6］DENG X G，ZHANG H X.Developing high－order accurate nonlinear schemes［J］.J.Comput.Phys，2000，165：22－44.

［7］WANG Z，HUANG G P.An essentially nonoscillatory high－order Pade－type（ENO－Pade）scheme［J］.J.Comput.Phys，2002，177：37－58.

［8］ADAMS N A，SHARIFF K.A high－resolution hybrid compact－ENO scheme for shock－turbulence interaction problems［J］.J.Comput.Phys，1996，127：27－51.

［9］PIROZZOLI S.Conservative hybrid compact－WENO schemes for shock－turbulence interaction［J］.J.Comput.Phys，2001，178：81－117.

［10］REN Y X，LIU M，ZHANG H.A characteristic－wise hybrid compact－WENO scheme for solving hyperbolic conservation laws［J］.J.Comput.Phys，2003，192：365－386.

［11］BORGES R，CARMONA M，COSTA B，et al.An improved weighted essentially non－oscillatory scheme for hyperbolic conservation laws［J］.J.Comput.Phys，2008，227：3192－3211.

［12］JIANG G S，SHU C W.Efficient implementation of weighted ENO schemes［J］.J.Comput.Phys，1996，126：202－228.

［13］WOODWARD P，COLLELA P.The numerical simulation of two dimensional fluid flow with strong shocks［J］.J.Comput.Phys，1984，54：115－173.