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

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

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

0 引 言

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

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

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

1 控制方程

含激波流場(chǎng)采用可壓縮歐拉方程進(jìn)行描述，其三維形式［18］如下：

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

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

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

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

2 數(shù)值方法

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

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

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

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

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

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

2.2 WENO-N3格式

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

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

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

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

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

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

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

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

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

2.3 改進(jìn)的3階WENO-N3格式

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

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

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

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

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

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

因?yàn)棣翫＞βC，令

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

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

2.3.2 極值點(diǎn)處的精度

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

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

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

3 數(shù)值驗(yàn)證

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

3.1 精度測(cè)試

該算例選自文獻(xiàn)［4］，計(jì)算初始條件為

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

3.2 激波與熵波相互作用

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

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

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

3.3 爆轟波相互作用

該問(wèn)題選自文獻(xiàn)［23］，其初始條件如式（37）所示：

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

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

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

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

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

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

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

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

4.1 初始條件

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

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

4.2 典型測(cè)點(diǎn)超壓對(duì)比

為了比較改進(jìn)格式WENO-N+3與WENO-JS3，WENO-N3格式對(duì)于爆炸載荷計(jì)算的差異性，給出爆炸波演化過(guò)程（圖6），并對(duì)壁面2個(gè)典型測(cè)點(diǎn)No.1和No.2的超壓時(shí)間歷程曲線(xiàn)進(jìn)行了對(duì)比分析。

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

5 結(jié) 論

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

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

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

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