改进的3阶WENO-N3格式及其在内爆炸载荷计算中的应用

徐维铮 ，吴卫国

1武汉理工大学高性能舰船技术教育部重点实验室，湖北 武汉 430063

2武汉理工大学交通学院，湖北 武汉 430063

0 引 言

高分辨率激波捕捉格式对含激波流场的数值模拟具有重要意义，其不但可以降低网格的规模，而且还能较好地分辨出流场中复杂的波系结构。Liu等［1］在本质无振荡（ENO）概念［2］的基础上提出WENO格式，将ENO格式改进为所有模板的加权平均，而权值则依据模板的光滑程度而定。Jiang等［3］通过引入新的光滑因子，构造了3阶和5阶WENO格式。

然而，Henrick等［4］指出传统的 5阶WENO格式在一阶极值点处精度会降为3阶。为了解决这个缺陷，采用映射函数法提出了5阶精度的WENO-M格式。Borges等［5］则通过线性组合低阶候选模板光滑因子构造全局高阶光滑因子的方式，提出5阶WENO-Z格式。之后，大量的研究［6-13］致力于改进高阶WENO格式（主要集中在5阶、7阶和9阶）在极值点处、间断处附近的精度。针对3阶WENO格式的改进，主要有以下几种：Yamaleev等［14］提出的能量稳定 ESWENO，Wu等［15-16］提出的WENO-N3格式以及WENO-NP3格式。近期，Acker等［17］指出，增大非光滑模板所对应的非线性权重能够改善格式的分辨率，并推导给出5阶WENO-Z+格式，数值试验表明，该格式具有比WENO-Z格式更小的耗散性。

本文将在文献［15，17］所做工作的基础上，提出改进格式WENO-N+3，并通过理论推导和数值算例的方式，验证改进格式WENO-N+3相较WENO-JS3，WENO-Z3，WENO-N3格式所具有的特性。最后，将改进格式WENO-N+3应用于封闭空间内爆炸载荷的数值计算。

1 控制方程

含激波流场采用可压缩欧拉方程进行描述，其三维形式［18］如下：

式中：ρ为密度；u，v，w为x，y，z方向上的速度分量；p为流体压力；E为单位体积流体的总能量；e为比内能；γ为气体的绝热指数，在本文数值计算中取为1.4。

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

例如，针对x方向，式（5）的数值离散形式为

式 中 ：分别为单元（xi-1/2，xi+1/2）的左、右界面对流项数值通量；h为x方向的均匀网格间距，在本文的数值计算中，x，y，z方向的网格间距相同。控制方程的空间离散采用下述的数值方法，时间项离散采用3阶TVD龙格库塔格式［19］。

2 数值方法

2.1 3阶WENO-JS格式数值重构过程

3阶WENO格式（WENO-JS3）的数值离散和推导过程如下［3］，为了简洁表述，仅给出了右界面通 量的重构过程，对于3阶WENO格式，的2种重构方式分别为：

利用上述2种2阶通量的加权组合计算最终具有3阶精度的数值通量即

对于光滑情形，有式（8）给出的形式既适合光滑流场也适合含间断流场，对于含激波的间断流场，式中的非线性权函数ωk需要根据下式求得：

式中，参数ε取值为10-6。光滑因子βk(k=0，1)的表达式为

文献［14，20］通过理论推导，给出3阶WENO格式达到收敛精度的充分条件为

2.2 WENO-N3格式

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

文献［15］首先指出，传统的3阶WENO-Z格式在极值点处会降阶，并提出了WENO-N3格式：

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

这里通过理论推导分析WENO-N3格式在1阶极值点处的计算精度，当存在1阶极值点时，β3在xi处泰勒级数展开为

式（10）中的光滑因子在xi处泰勒级数展开为

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

再根据加权法，则可得右界面通量所对应的权函数

同理，可得左界面通量所对应的权函数

式（21）和式（22）显示，式（23）和式（24）显示，均不满足充分条件，因此，WENO-N3格式在光滑流场极值点处的精度将降低。文献［15］中的极值点精度测试显示，3阶WENO-N3格式在1阶极值点处将降低为1阶精度。

2.3 改进的3阶WENO-N3格式

文献［17］表示，在相对稀疏的网格下，增大非光滑模板的非线性权重，相比单纯提高格式在极值点处的精度，前者能给出分辨率更好的结果，在该构造思想的启发下，直接给出了改进格式WENO-N+3：

式中，λ=h0=1，其取值依据3.2小节的算例分析。

2.3.1 构造原理证明

改进格式的构造思想是增大非光滑模板的非线性权重，从而降低格式的耗散性，提高格式的分辨率，以下给出基本证明［22］。

基本假定：3阶WENO格式有2个子模板SC，SD，其中SC为光滑模板，SD为含间断模板，意味着光滑因子βC＜βD。

证明：不考虑数值很小的参数ε，根据式（14）和式（25）可得

因为βD＞βC，令

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

根据式（27）和式（29），最终得证改进格式WENO-N+3相较WENO-N3格式增大了非光滑模板的非线性权重：

2.3.2 极值点处的精度

通过理论推导分析改进的WENO-N+3格式在1阶极值点处的计算精度，通过系列推导，可得右界面通量所对应的权函数

同理，可得左界面通量所对应的权函数

根据式（31）～式（34）可知，不能满足式（11）的充分条件，因此可以预知其在极值点处将降阶。

3 数值验证

为了进一步考察改进格式WENO-N+3的计算性能，选取线性精度测试、激波与熵波相互作用、双爆轰波碰撞、瑞利—泰勒不稳定性等经典算例进行自主编程计算，并将该格式计算结果与格式WENO-JS3，WENO-Z3和WENO-N3进行对比分析。

3.1 精度测试

该算例选自文献［4］，计算初始条件为

其包含2个1阶极值点。表1给出了WENOJS3，WENO-Z3，WENO-N3和 WENO-N+3格式的L1范数来计算误差和精度，其中N为网格数。可知，改进格式WENO-N+3，WENO-N3格式基本上具有相同的精度，在极值点处均没有达到3阶收敛精度，与2.2和2.3.1节的理论分析相一致。

3.2 激波与熵波相互作用

该算例初始条件［19］如式（36）所示，网格数为800，计算结束时间为1.8。图1给出了计算结束时刻密度曲线图及局部放大图。

表1 针对初始条件，在结束时间为2时不同数值计算格式L1误差和精度比较Table 1 A comparison study ofL1（error and order）for different schemes with initial condition at t=2

为了考察不同参数λ对改进格式WENO-N+3的影响规律，针对该算例，选用4个不同的数值：λ=h0，h1/2，h2/3，h3/3进行数值计算。图2给出相应的计算结果。根据图2发现，随着参数λ的增大，格式给出了耗散更低的计算结果，其中参数λ=h0=1给出了较优的计算结果。由于WENO格式需要在激波附近具备一定的数值耗散以抑制数值振荡，因此，过大的增加参数λ的数值将会造成一定的虚假振荡。本文的数值计算结果表明，参数λ=h0=1是改进格式WENO-N+3综合权衡耗散性和计算稳定性后较为折中的参数，这也是本文选取参数λ=h0=1的主要考虑。

3.3 爆轰波相互作用

该问题选自文献［23］，其初始条件如式（37）所示：

计算网格数为800，计算结束时间为0.038，两端边界条件取为刚性反射边界。

根据图1和图3中的一维算例结果可知，改进格式WENO-N+3相较于WENO-N3，WENO-Z3和WENO-JS3格式具有更低的耗散性。

3.4 瑞利—泰勒不稳定性

该问题主要描述重力场作用下，重流体加速进入轻流体界面失稳的过程。文献［17，24-25］也采用该算例探讨过数值方法的分辨率特性。计算域设置为［0，0.25］×［0，1］，计算初始条件如下所示：

该算例中，绝热指数γ取值为5/3。在欧拉方程y方向的动量方程和能量方程右端分别加入ρ，ρv作为源项来考虑重力效应。左、右边界设置成反射边界条件，顶部和底部边界条件分别为（ρ，u，v，p）=（1，0，0，2.5），（ρ，u，v，p）=（2，0，0，1）。网格数划分为240×960，计算结束时间为1.95。

各格式（WENO-JS3，WENO-Z3，WENO-N3和WENO-N+3）的密度等值线图如图4所示，共绘制出15条等值线，其取值区间为［0.952 269，2.145 89］。改进格式WENO-N+3在接触间断附近具有更低的耗散性，能给出分辨率更清晰的计算结果。

4 三维封闭空间内爆炸

将改进格式WENO-N+3应用于正方体封闭空间内柱形炸药爆炸载荷计算，并将计算结果与WENO-JS3，WENO-N3进行比较分析。

4.1 初始条件

封闭空间尺寸为800 mm×800 mm×800 mm，壁面上设置2个测点，分别编号为No.1和No.2，对爆炸超压时间历程进行输出，如图5（a）所示。

柱形炸药位于中间，瞬时爆轰后的高压、高密度气体参数为：半径50 mm，高度140 mm，密度1 630 kg/m3，压力 3.057 9×109Pa。计算初始条件如图5（b）所示，图中P为爆炸压力。网格数选取为40×40×40。壁面边界条件设置为刚性反射边界。

4.2 典型测点超压对比

为了比较改进格式WENO-N+3与WENO-JS3，WENO-N3格式对于爆炸载荷计算的差异性，给出爆炸波演化过程（图6），并对壁面2个典型测点No.1和No.2的超压时间历程曲线进行了对比分析。

三维封闭空间内部爆炸载荷由于爆炸波的壁面反射和叠加，呈现出不断衰减的多峰值特点。从图7的对比结果可以看出，改进格式WENO-N+3相较WENO-JS3和WENO-N3给出了较优的结果，尤其能以较低的耗散捕捉内爆炸载荷的前几个峰值。

5 结 论

通过本文的研究，主要得到以下2点结论：

1）改进格式WENO-N+3相较其他格式（WENO-JS3，WENO-Z3，WENO-N3）具有较低的耗散，提高了对复杂流场结构的分辨率。

2）改进格式WENO-N+3在相同的计算网格下能给出较高的冲击波峰值，其用于内爆炸载荷的数值计算具有一定的可行性。

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