徐維錚 ，吳衛(wèi)國

1武漢理工大學(xué)高性能艦船技術(shù)教育部重點實驗室，湖北 武漢 430063

2武漢理工大學(xué)交通學(xué)院，湖北 武漢 430063

0 引 言

高分辨率激波捕捉格式對含激波流場的數(shù)值模擬具有重要意義，其不但可以降低網(wǎng)格的規(guī)模，而且還能較好地分辨出流場中復(fù)雜的波系結(jié)構(gòu)。Liu等［1］在本質(zhì)無振蕩（ENO）概念［2］的基礎(chǔ)上提出WENO格式，將ENO格式改進為所有模板的加權(quán)平均，而權(quán)值則依據(jù)模板的光滑程度而定。Jiang等［3］通過引入新的光滑因子，構(gòu)造了3階和5階WENO格式。

然而，Henrick等［4］指出傳統(tǒng)的 5階WENO格式在一階極值點處精度會降為3階。為了解決這個缺陷，采用映射函數(shù)法提出了5階精度的WENO-M格式。Borges等［5］則通過線性組合低階候選模板光滑因子構(gòu)造全局高階光滑因子的方式，提出5階WENO-Z格式。之后，大量的研究［6-13］致力于改進高階WENO格式（主要集中在5階、7階和9階）在極值點處、間斷處附近的精度。針對3階WENO格式的改進，主要有以下幾種：Yamaleev等［14］提出的能量穩(wěn)定 ESWENO，Wu等［15-16］提出的WENO-N3格式以及WENO-NP3格式。近期，Acker等［17］指出，增大非光滑模板所對應(yīng)的非線性權(quán)重能夠改善格式的分辨率，并推導(dǎo)給出5階WENO-Z+格式，數(shù)值試驗表明，該格式具有比WENO-Z格式更小的耗散性。

本文將在文獻［15，17］所做工作的基礎(chǔ)上，提出改進格式WENO-N+3，并通過理論推導(dǎo)和數(shù)值算例的方式，驗證改進格式WENO-N+3相較WENO-JS3，WENO-Z3，WENO-N3格式所具有的特性。最后，將改進格式WENO-N+3應(yīng)用于封閉空間內(nèi)爆炸載荷的數(shù)值計算。

1 控制方程

含激波流場采用可壓縮歐拉方程進行描述，其三維形式［18］如下：

式中：ρ為密度；u，v，w為x，y，z方向上的速度分量；p為流體壓力；E為單位體積流體的總能量；e為比內(nèi)能；γ為氣體的絕熱指數(shù)，在本文數(shù)值計算中取為1.4。

在式（1）的每個方向上均可以看成雙曲守恒律方程：

例如，針對x方向，式（5）的數(shù)值離散形式為

式 中 ：分別為單元（xi-1/2，xi+1/2）的左、右界面對流項數(shù)值通量；h為x方向的均勻網(wǎng)格間距，在本文的數(shù)值計算中，x，y，z方向的網(wǎng)格間距相同。控制方程的空間離散采用下述的數(shù)值方法，時間項離散采用3階TVD龍格庫塔格式［19］。

2 數(shù)值方法

2.1 3階WENO-JS格式數(shù)值重構(gòu)過程

3階WENO格式（WENO-JS3）的數(shù)值離散和推導(dǎo)過程如下［3］，為了簡潔表述，僅給出了右界面通 量的重構(gòu)過程，對于3階WENO格式，的2種重構(gòu)方式分別為：

利用上述2種2階通量的加權(quán)組合計算最終具有3階精度的數(shù)值通量即

對于光滑情形，有式（8）給出的形式既適合光滑流場也適合含間斷流場，對于含激波的間斷流場，式中的非線性權(quán)函數(shù)ωk需要根據(jù)下式求得：

式中，參數(shù)ε取值為10-6。光滑因子βk(k=0，1)的表達式為

文獻［14，20］通過理論推導(dǎo)，給出3階WENO格式達到收斂精度的充分條件為

2.2 WENO-N3格式

3階 WENO-Z 格式［21］（WENO-Z3）的具體形式如下：

文獻［15］首先指出，傳統(tǒng)的3階WENO-Z格式在極值點處會降階，并提出了WENO-N3格式：

式中：τ，τN為高階全局光滑因子；β3表示3階WENO格式全局模板（xi-1，xi，xi+1）的光滑因子：

這里通過理論推導(dǎo)分析WENO-N3格式在1階極值點處的計算精度，當存在1階極值點時，β3在xi處泰勒級數(shù)展開為

式（10）中的光滑因子在xi處泰勒級數(shù)展開為

將式（17）～式（18）代入式（14）和式（15），可得

再根據(jù)加權(quán)法，則可得右界面通量所對應(yīng)的權(quán)函數(shù)

同理，可得左界面通量所對應(yīng)的權(quán)函數(shù)

式（21）和式（22）顯示，式（23）和式（24）顯示，均不滿足充分條件，因此，WENO-N3格式在光滑流場極值點處的精度將降低。文獻［15］中的極值點精度測試顯示，3階WENO-N3格式在1階極值點處將降低為1階精度。

2.3 改進的3階WENO-N3格式

文獻［17］表示，在相對稀疏的網(wǎng)格下，增大非光滑模板的非線性權(quán)重，相比單純提高格式在極值點處的精度，前者能給出分辨率更好的結(jié)果，在該構(gòu)造思想的啟發(fā)下，直接給出了改進格式WENO-N+3：

式中，λ=h0=1，其取值依據(jù)3.2小節(jié)的算例分析。

2.3.1 構(gòu)造原理證明

改進格式的構(gòu)造思想是增大非光滑模板的非線性權(quán)重，從而降低格式的耗散性，提高格式的分辨率，以下給出基本證明［22］。

基本假定：3階WENO格式有2個子模板SC，SD，其中SC為光滑模板，SD為含間斷模板，意味著光滑因子βC＜βD。

證明：不考慮數(shù)值很小的參數(shù)ε，根據(jù)式（14）和式（25）可得

因為βD＞βC，令

將式（28）代入式（27），可得

根據(jù)式（27）和式（29），最終得證改進格式WENO-N+3相較WENO-N3格式增大了非光滑模板的非線性權(quán)重：

2.3.2 極值點處的精度

通過理論推導(dǎo)分析改進的WENO-N+3格式在1階極值點處的計算精度，通過系列推導(dǎo)，可得右界面通量所對應(yīng)的權(quán)函數(shù)

同理，可得左界面通量所對應(yīng)的權(quán)函數(shù)

根據(jù)式（31）～式（34）可知，不能滿足式（11）的充分條件，因此可以預(yù)知其在極值點處將降階。

3 數(shù)值驗證

為了進一步考察改進格式WENO-N+3的計算性能，選取線性精度測試、激波與熵波相互作用、雙爆轟波碰撞、瑞利—泰勒不穩(wěn)定性等經(jīng)典算例進行自主編程計算，并將該格式計算結(jié)果與格式WENO-JS3，WENO-Z3和WENO-N3進行對比分析。

3.1 精度測試

該算例選自文獻［4］，計算初始條件為

其包含2個1階極值點。表1給出了WENOJS3，WENO-Z3，WENO-N3和 WENO-N+3格式的L1范數(shù)來計算誤差和精度，其中N為網(wǎng)格數(shù)?？芍?，改進格式WENO-N+3，WENO-N3格式基本上具有相同的精度，在極值點處均沒有達到3階收斂精度，與2.2和2.3.1節(jié)的理論分析相一致。

3.2 激波與熵波相互作用

該算例初始條件［19］如式（36）所示，網(wǎng)格數(shù)為800，計算結(jié)束時間為1.8。圖1給出了計算結(jié)束時刻密度曲線圖及局部放大圖。

表1 針對初始條件，在結(jié)束時間為2時不同數(shù)值計算格式L1誤差和精度比較Table 1 A comparison study ofL1（error and order）for different schemes with initial condition at t=2

為了考察不同參數(shù)λ對改進格式WENO-N+3的影響規(guī)律，針對該算例，選用4個不同的數(shù)值：λ=h0，h1/2，h2/3，h3/3進行數(shù)值計算。圖2給出相應(yīng)的計算結(jié)果。根據(jù)圖2發(fā)現(xiàn)，隨著參數(shù)λ的增大，格式給出了耗散更低的計算結(jié)果，其中參數(shù)λ=h0=1給出了較優(yōu)的計算結(jié)果。由于WENO格式需要在激波附近具備一定的數(shù)值耗散以抑制數(shù)值振蕩，因此，過大的增加參數(shù)λ的數(shù)值將會造成一定的虛假振蕩。本文的數(shù)值計算結(jié)果表明，參數(shù)λ=h0=1是改進格式WENO-N+3綜合權(quán)衡耗散性和計算穩(wěn)定性后較為折中的參數(shù)，這也是本文選取參數(shù)λ=h0=1的主要考慮。

3.3 爆轟波相互作用

該問題選自文獻［23］，其初始條件如式（37）所示：

計算網(wǎng)格數(shù)為800，計算結(jié)束時間為0.038，兩端邊界條件取為剛性反射邊界。

根據(jù)圖1和圖3中的一維算例結(jié)果可知，改進格式WENO-N+3相較于WENO-N3，WENO-Z3和WENO-JS3格式具有更低的耗散性。

3.4 瑞利—泰勒不穩(wěn)定性

該問題主要描述重力場作用下，重流體加速進入輕流體界面失穩(wěn)的過程。文獻［17，24-25］也采用該算例探討過數(shù)值方法的分辨率特性。計算域設(shè)置為［0，0.25］×［0，1］，計算初始條件如下所示：

該算例中，絕熱指數(shù)γ取值為5/3。在歐拉方程y方向的動量方程和能量方程右端分別加入ρ，ρv作為源項來考慮重力效應(yīng)。左、右邊界設(shè)置成反射邊界條件，頂部和底部邊界條件分別為（ρ，u，v，p）=（1，0，0，2.5），（ρ，u，v，p）=（2，0，0，1）。網(wǎng)格數(shù)劃分為240×960，計算結(jié)束時間為1.95。

各格式（WENO-JS3，WENO-Z3，WENO-N3和WENO-N+3）的密度等值線圖如圖4所示，共繪制出15條等值線，其取值區(qū)間為［0.952 269，2.145 89］。改進格式WENO-N+3在接觸間斷附近具有更低的耗散性，能給出分辨率更清晰的計算結(jié)果。

4 三維封閉空間內(nèi)爆炸

將改進格式WENO-N+3應(yīng)用于正方體封閉空間內(nèi)柱形炸藥爆炸載荷計算，并將計算結(jié)果與WENO-JS3，WENO-N3進行比較分析。

4.1 初始條件

封閉空間尺寸為800 mm×800 mm×800 mm，壁面上設(shè)置2個測點，分別編號為No.1和No.2，對爆炸超壓時間歷程進行輸出，如圖5（a）所示。

柱形炸藥位于中間，瞬時爆轟后的高壓、高密度氣體參數(shù)為：半徑50 mm，高度140 mm，密度1 630 kg/m3，壓力 3.057 9×109Pa。計算初始條件如圖5（b）所示，圖中P為爆炸壓力。網(wǎng)格數(shù)選取為40×40×40。壁面邊界條件設(shè)置為剛性反射邊界。

4.2 典型測點超壓對比

為了比較改進格式WENO-N+3與WENO-JS3，WENO-N3格式對于爆炸載荷計算的差異性，給出爆炸波演化過程（圖6），并對壁面2個典型測點No.1和No.2的超壓時間歷程曲線進行了對比分析。

三維封閉空間內(nèi)部爆炸載荷由于爆炸波的壁面反射和疊加，呈現(xiàn)出不斷衰減的多峰值特點。從圖7的對比結(jié)果可以看出，改進格式WENO-N+3相較WENO-JS3和WENO-N3給出了較優(yōu)的結(jié)果，尤其能以較低的耗散捕捉內(nèi)爆炸載荷的前幾個峰值。

5 結(jié) 論

通過本文的研究，主要得到以下2點結(jié)論：

1）改進格式WENO-N+3相較其他格式（WENO-JS3，WENO-Z3，WENO-N3）具有較低的耗散，提高了對復(fù)雜流場結(jié)構(gòu)的分辨率。

2）改進格式WENO-N+3在相同的計算網(wǎng)格下能給出較高的沖擊波峰值，其用于內(nèi)爆炸載荷的數(shù)值計算具有一定的可行性。

［1］LIU X D，OSHER S，CHAN T.Weighted essentially non-oscillatory schemes［J］.Journal of Computational Physics，1994，115（1）：200-212.

［2］HARTEN A，ENGQUIST B，OSHER S，et al.Uniformly high order accurate essentially non-oscillatory schemes，III［J］.Journal of Computational Physics，1987，71（2）：231-303.

［3］JIANG G S，SHU C W.Efficient implementation of weighted ENO schemes［J］.Journal of Computational Physics，1995，126（1）：202-228.

［4］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.

［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］YAMALEEV N K，CARPENTER M H.A systematic methodology for constructing high-order energy stable WENO schemes［J］.Journal of Computational Physics，2009，228（1）：4248-4272.

［7］FAN P.High order weighted essentially nonoscillatory WENO-ηschemes for hyperbolic conservation laws［J］.Journal of Computational Physics，2014，269：355-385.

［8］FAN P，SHEN Y Q，TIAN B L，et al.A new smoothness indicator for improving the weighted essentially non-oscillatory scheme［J］.Journal of Computational Physics，2014，269：329-354.

［9］FENG H， HUANG C， WANG R.An improved mapped weighted essentially non-oscillatory scheme［J］.Applied Mathematics and Computation，2014，232：453-468.

［10］SHEN Y Q，ZHA G H.Improvement of weighted essentially non-oscillatory schemes near discontinuities［J］.Computers&Fluids，2014，96：1-9.

［11］KIM C H，HA Y，YOON J.Modified non-linear weights for fifth-order weighted essentially non-oscillatory schemes［J］.Journal of Scientific Computing，2016，67（1）：299-323.

［12］MA Y K，YAN Z G，ZHU H J.Improvement of multistep WENO scheme and its extension to higher orders of accuracy［J］.International Journal for Numerical Methods in Fluids，2016，82（12）：818-838.

［13］WANG R，F(xiàn)ENG H，HUANG C.A New mapped weighted essentially non-oscillatory method using rational mapping function［J］.Journal of Scientific Computing，2016，67（2）：540-580.

［14］YAMALEEV N K，CARPENTER M H.Third-order energy stable WENO scheme［J］.Journal of Computational Physics，2013，228（8）：3025-3047.

［15］WU X S， ZHAO Y X.A high-resolution hybrid scheme for hyperbolic conservation laws［J］.International Journal for Numerical Methods in Fluids，2015，78（3）：162-187.

［16］WU X S，LIANG J H，ZHAO Y X.A new smoothness indicator for third-order WENO scheme［J］.International Journal for Numerical Methods in Fluids，2016，81（7）：451-459.

［17］ACKER F，DE R BORGES RB，COSTA B.An improved WENO-Z scheme［J］.Journal of Computational Physics，2016，313：726-753.

［18］TORO E F.Riemann solvers and numerical methods for fluid dynamics：a practical introduction［M］.Berlin Heidelberg：Springer，1999：87-114.

［19］SHU C W，OSHER S.Efficient implementation of essentially non-oscillatory shock-capturing schemes，II［J］.Journal of Computational Physics，1989，83（1）：32-78.

［20］GANDE N R，RATHOD Y，RATHAN S.Third-order WENO scheme with a new smoothness indicator［J］.International Journal for Numerical Methods in Fluids，2017，85（2）：90-112.

［21］DON W S，BORGES R.Accuracy of the weighted essentially non-oscillatory conservative finite difference schemes［J］.Journal of Computational Physics，2013，250：347-372.

［22］XU W Z，WU W G.An improved third-order WENO-Z scheme［J］.Journal of Scientific Computing，2018，75：1808-1841.

［23］WOODWARD P，COLELLA P.The numerical simulation of two-dimensional fluid flow with strong shocks［J］.Journal of Computational Physics，1984，54（1）：115-173.

［24］SHI J，ZHANG Y T，SHU C W.Resolution of high order WENO schemes for complicated flow structures［J］.Journal of Computational Physics，2003，186（2）：690-696.

［25］HU X Y，WANG Q，ADAMS N A.An adaptive central-upwind weighted essentially non-oscillatory scheme［J］.Journal of Computational Physics，2010，229（23）：8952-8965.