求解双曲守恒律方程的三阶修正模板WENO格式*

王亚辉

（郑州师范学院 数学与统计学院，郑州 450044）

引 言

非线性双曲守恒律的解通常包含间断和光滑的小尺度结构.在不产生虚假振荡的情况下捕捉间断，并尽可能多地求解小尺度结构，是数值格式的必要条件.在求解双曲守恒律方程的众多数值格式中，加权本质无振荡(WENO)格式已成为最常用的方法之一.WENO 是求解不连续双曲守恒律的一种常用方法，其抑制虚假振荡的成功归功于一组动态候选模板的低阶多项式的非线性组合方法.加权过程是通过根据每个候选模板的光滑度评估其局部通量的贡献来执行的，具体地说，WENO 格式中的光滑度指示子可以驱使格式向最优方向发展，并通过分配间断模板的基本零权重来避免跨间断的插值.

由Liu 等首先介绍的WENO 格式[1]使用了所有候选模板的非线性权来凸组合ENO 格式[2-5]，以提高解在光滑区域的精度，而不丢失ENO 格式在间断附近的无振荡特性.该过程通过根据候选模板上解的光滑度指示子对局部通量的贡献进行加权来执行，使得包含间断局部通量的模板权重基本上为零.Jiang 和Shu[6]设计了一个光滑度指标子，即所有低阶多项式的各阶导数的L2范数的归一化平方和，从而发展了经典的有限差分WENO(WENO-JS)格式.Henrick 等[7]提出了WENO-M 格式，用映射函数修改了权，使之满足格式最优收敛阶的充分必要条件.Borges 等[8]发展了WENO-Z 格式，该方法将全局高阶光滑度指示子纳入权重的构造中.WENO-Z 格式可以达到与WENO-M 格式相同的精度.由于该格式赋予间断模板较大的权重，因此发展了七阶到十一阶保单调WENO 格式.Gerolymos 等[9]将WENO-Z 格式推广到十一阶WENO 格式，发展了七阶到十一阶WENO-Z 格式.

与高阶WENO 格式相比，三阶WENO 格式具有许多优点.例如，它对激波问题具有很强的鲁棒性，使用较少的网格点，从而降低了边界处理和求解的难度，很容易地推广到非结构网格上.尽管三阶WENO 格式具有上述优点，但三阶WENO 格式很少受到关注.这是因为经典的三阶WENO 格式在光滑区域的精度在2 到3 阶之间，而在临界点退化为1 阶，而三阶WENO 格式由于成本较低，是理想的格式.利用WENO-Z 格式的加权过程，可以解决提高分辨率和恢复最优收敛阶的问题.通过这种方法，Wu 和Zhao[10]开发了一个改进的三阶WENO(WENO-N3)格式.他们设计了一个四阶参考光滑度指示子 τ4N，由候选模板的光滑度指示子与全局模板的光滑度指示子的线性组合得到的.然而 τ4N不能使格式在临界点附近达到最优三阶收敛.为了弥补这个缺陷，Wu 等[11]用τ4Np的p幂形式改进了参考光滑度指示子，得到了WENO-Np3 格式.Xu 和Wu[12]提出了另一个四阶参考光滑度指示子 τ4P，这样得到的WENO-P3 格式耗散性更小.然而，该算法也只能在没有临界点的光滑区域实现三阶精度，在有临界点的光滑区域降到一阶准确度.Wang 等[13]给出了三阶WENO 格式的参考光滑度指示子 τ4Rp，理论分析表明，所得到的WENO-Rp3 格式(p=1.5)在光滑区域(包括临界点)能达到三阶精度，并且WENO-R3 格式(p=1)保持了ENO 性质.关于三阶WENO 格式的其他改进格式可参见文献[14-15].

为了降低经典三阶WENO 格式的数值耗散，提高数值精度，这里我们提出了一种新的WENO 格式的改进模板近似方法.这里我们所采取的主要策略与文献[10-13]的不同，不是构造更高阶的参考光滑度指示子，而是改进经典WENO-JS3 格式中各候选模板上数值通量的一阶多项式逼近，通过加入二次项使模板逼近达到三阶精度，并且计算了相应的候选通量.它能使所得格式(WENO-MS)在包括一阶临界点在内的光滑区域内达到三阶收敛点数.为了抑制强激波和间断区域的明显数值振荡，在不影响计算精度的前提下，我们在候选模板通量的修正项中加入可调函数 φ，并给出了一系列一维和二维数值算例来验证新格式的性能.数值结果表明，与WENO-JS3、WENO-Z3 格式相比，WENO-MS 格式对精细光滑结构具有相当或更高的分辨率.

本文的主要结构如下：第1 节简要回顾了一维标量守恒律的三阶WENO 格式；第2 节提出了新的修正模板WENO 近似方法，并给出了WENO-MS 格式；第3 节给出了数值算例来说明新格式的性能；第4 节对全文进行了总结.

1 三阶WENO 格式综述

在这一节中，我们简要介绍了一维标量守恒律的三阶WENO 有限差分格式：

设Ij是 给定域的一个分区，第j个Ij:=[xj−1/2，xj+1/2].Ij的 中心表示为函数f在位置xj的 值由下标j表示，即fj=f(xj).我们假设网格点{xj+1/2}均 匀分布，单元大小Δx=xj+1/2−xj−1/2.

为了构造方程(1)的守恒有限差分格式，可以通过隐函数h(x)来定义数值通量[5]：

对方程(2)求关于x的 导数得到

因此，式(1)可由常微分方程组近似，其中空间导数由守恒有限差分近似，产生半离散形式：

图1 三阶WENO 数值通量的模板Fig.1 Stencils for the 3rd-order WENO numerical flux

1.1 三阶WENO 格式

它们结合起来定义了格式的数值通量：

这里的非线性权重 ωk定义为

其中，常系数dk被称为理想权重，因为其与的 线性组合保持最优三阶收敛到h(xj+1/2)，即

dk的具体值如下所示(参见文献[16])：

式(9)中的0 <ϵ ≪1是 为防止被零除而引入的小正参数， βk是第k个模板的光滑度指示子.Jiang 和Shu[6]提出的光滑度指示子为

其具体形式为

Borges 等[8]提出了一个WENO-Z 权：

其中 ϵ =10−40，τ3=|β1−β0|是 一个三阶参考光滑探测子，它驱使非线性权重 ωk(式(13)的WENO-Z)向最优权重dk靠近，并且比J-S 权重(9)快.当幂q=2时，WENO-Z 格式在一阶临界点处可以达到最优的三阶，但在间断附近具有更多的耗散.一个不那么耗散的方法是设计一个足够的高阶参考光滑探测子，并且总是使用幂q=1.

为了要为WENO-Z 权设计一个合适的参考光滑探测子(13)，需要一个充分条件来恢复三阶WENO 格式的最优阶精度.文献[12]给出了三阶WENO 格式达到三阶精度的一个条件：

方程(14)仅仅是一个充分条件，但不是必要条件.然而，这一条件为设计高阶参考光滑性指示子提供了有用的信息.

2 新的修正三阶WENO 格式

2.1 修正模板近似

“MS”表示“修正模板”.

注1修改后的候选通量(25)依赖于整个三点模板，实际上，通过简单的计算不难发现，它就是三阶迎风格式，因此新的格式(26)失去了ENO 性质.为了抑制强激波和不连续区域内明显的数值振荡，在这里，应用了一个可调函数 φ来限制式(25)中最后一项的影响，使新发展的格式具有ENO 性质.修正后的候选通量(25)变为

这里φ 应该满足以下条件：

(a) 如果模板S3={xj−1,xj,xj+1}是 一个间断模板，则φ 是一个很小的值；

(b) 如果模板S3={xj−1,xj,xj+1}是 一个光滑模板，则φ 是 接近1，并且φ 不会影响新格式的收敛精度.

所以φ 的具体形式定义如下：

因而在接下来的数值算例中我们取定q=2.

3 数值测试

在本节中，我们将在几个数值例子中验证所提出的WENO-MS 格式(WENO-MS-JS 和WENO-MS-Z)的性能，并与WENO-JS 和WENO-Z 格式进行比较.数值算例从简单标量对流方程求解开始，然后是一维和二维Euler 方程的数值求解.本文所有的格式和算例均采用局部Lax-Friedrichs 通量分裂.时间推进使用三阶TVD-Runge-Kutta 方法[5]：

其 中 σ是CFL(Courant-Friedrichs-Lewy)条件数.

3.1 标量测试问题

考虑具有周期边界的标量对流方程条件并用不同的初始数据测试任意初始剖面的传播，包含跳跃间断和角点.

例1 (线性方程)让我们考虑如下线性对流方程：

初值条件为

u(x,0)=u0(x).

我们首先在两组初始数据的线性方程(31)上检验了WENO-MS 格式的数值收敛阶：

时间步长为Δt=0.6Δx， 保证与空间精度兼容.误差的L1范 数是通过与t=2的精确解比较得出的，根据以下式子计算：

表1 和表2 分别显示t=2时 初值(32a)的L1、L∞误差和收敛阶.我们发现误差随着WENO-JS、WENOZ3、WENO-MS-JS3 和WENO-MS-Z3 的顺序而减小.新的WENO-MS 格式能达到最优三阶收敛精度.表3 和表4显示初值(32b)的结果.尽管有临界点的出现，但是新的WENO-MS 格式的数值收敛阶仍然能达到最优阶.相比之下，新的WENO-MS 格式的性能优于其他WENO 格式.

表1 线性对流方程(31)在初值(32a)下，不同格式在t =2.0的L1 误差和收敛阶Table 1 L1 errors and convergence rates at t =2.0 of different schemes for linear advection eq.(31) with initial data eq.(32a)

表2 线性对流方程(31)在初值(32a)下，不同格式在t =2.0 的 L∞误差和收敛阶Table 2 L ∞ errors and convergence rates at t =2.0 of different schemes for linear advection eq.(31) with initial data eq.(32a)

表3 线性对流方程(31)在初值(32b)下，不同格式在t =2.0的L1 误差和收敛阶Table 3 L1 errors and convergence rates at t =2.0 of different schemes for linear advection eq.(31) with initial data eq.(32b)

表4 线性对流方程(31)在初值(32b)下，不同格式在t =2.0的L∞误差和收敛阶Table 4 L∞ errors and convergence rates at t= 2.0 of different schemes for linear advection eq.(31) with initial data eq.(32b)

然后我们讨论WENO-MS 格式在长时间计算中的性能.具体来说，在式(31)上用

这是一个分段正弦函数，在x=0处 有一个跳跃，其中CFL 数为0.5，我们将其解到t=40处，以观察格式在跳跃间断处的行为.图2 显示了使用网格点N=400的不同WENO 格式的数值解与解析解的比较.我们可以看到，WENO-MS-JS3 和WENO-MS-Z3 格式的耗散比WENO-JS3 和WENO-Z3 格式小.

图2 线性对流方程(31)在初值(33)下数值解与解析解的比较，t =40 Fig.2 Comparison of the analytical solution with the numerical solutions of linear advection eq.(31) with initial value (33), at t =40

我们进一步考虑初始条件：

我们求解方程(31)，初始条件(34)在计算域 −1≤x≤1中 时间为t=41， 其中 Δx=0.005，CFL 数为0.5.数值结果如图3 所示.易见在间断附近，WENO-MS-JS3 和WENO-MS-Z3 的耗散比WENO-JS3 和WENO-Z3 小.

图3 线性对流方程(31)在初值(34)下不同格式的数值解与解析解的比较，t =41 ，N=400Fig.3 Comparison of the analytical solution with the numerical solutions of linear advection eq.(31) with initial value (34), at t=41, N=400

最后，我们求解线性对流方程(31)，初始条件包括Gauss 波、方波、三角波和半椭圆波：

其中G(x,β,z)=e−β(x−z)2，F(x,α,a)=a=0.5,z=−0.7，δ=0.005，α=10和 β=(ln2)/(36δ2).这种情况的解包括接触间断，角点奇异和光滑区域.我们用时间步长 Δt=CCFLΔx，CFL 数CCFL为0.5，计算到最终时间为t=6.数值结果如图4 所示.结果表明，WENO-MS 格式的性能优于WENO-JS3 和WENO-Z3 格式，特别是在间断面附近.所有改进的三阶WENO 格式都比经典的WENO-JS3 格式具有更好的分辨率，其中WENO-MS 格式的耗散最小.

图4 线性对流方程(31)在初值(35)下不同格式的数值解与解析解的比较，t =6 ，N=400Fig.4 Comparison of the analytical solution with the numerical solutions of linear advection eq.(31) with initial value (35), at t=6, N=400

3.2 一维Euler 系统

在这一小节中，我们将给出双曲守恒律方程的数值结果.考虑一维Euler 方程组：

其中

U=(ρ,ρu,E)T,F(U)=(ρu,ρu2+p,u(E+p))T.

状态方程为

ρ,u,p和E分别是密度、速度、压强和总能， γ 是比热比.相应的Jacobi 矩阵A(U)=∂F/∂U的特征值为

λ1=u−c, λ2=u, λ3=u+c,

这里c=(γp/ρ)1/2为声速.采用特征场的特征分解与局部Lax-Friedrichs 分裂，并将WENO 格式推广到一维Euler 系统.

例2 (一维 Riemann 问题)我们考虑具有如下形式的三个经典初值问题：

① Sod 激波管[17-18]

其中γ =1.4,CCFL=0.6， 网格尺度为Δx=0.005， 计算时间为t=0.2.密度场的数值结果与精确的Riemann 解进行了比较，如图5 所示.可以看出，WENO-MS-Z3 格式的性能优于WENO-MS-JS3、WENO-Z3 和WENO-JS3 格式.与其他WENO 格式相比，WENO-MS-Z3 格式在很好地捕捉激波的同时，具有最小的耗散性.

图5 Sod 激波管问题[16]的数值结果，t =0.2 ，N =200，CCFL=0.6Fig.5 Numerical results of the Sod problem[16], at t=0.2 ， N=200，CCFL=0.6

② 冲击波相互作用[17-18]

其中γ =1.4,CCFL=0.6.我们用网格尺度为 Δx=0.001 25计 算到t=0.038.密度场的数值结果如图6 所示.由于精确解是未知的，参考解是通过使用五阶WENO-JS 格式[6]在3 201 个网格点上获得的.易知, 在接触间断处的耗散的顺序是WENO-MS-Z3

图6 冲击波的相互作用的数值结果，t =0.038 ，N =800，CCFL=0.6Fig.6 Numerical results of interacting blast waves, at t =0.038 ， N =800，CCFL = 0.6

③ 一维激波等熵波作用(Shu-Osher 问题)[5]我们计算[ −5,5]区域上的近似解，这里采用周期边界条件.初始条件为

其中 γ =1.4,CCFL=0.6， ϵ 和k分 别是熵波的振幅和波数.这里，我们设置ϵ =0.2，k=5.在这个问题中，一个右移的超音速(Ma=3)激波与密度扰动中的正弦波相互作用，产生一个具有光滑结构和不连续性的流场.这股气流在波数高于初始密度变化波数k的右行激波后面诱导波迹.由于精确解是未知的，参考解是通过使用五阶WENO-JS 格式[6]在3201 个网格点上获得的.图7 将密度剖面的数值结果与参考解进行了比较.所有的WENO 格式都比WENO-JS3 格式的耗散性小.WENO-MS 格式比WENO-Z3 和WENO-JS3 格式更好地捕捉了右行激波后面的高频波，而且WENO-MS 格式的分辨率最接近WENO-JS5 格式.

图7 激波等熵波相互作用(Shu-Osher)[5]的密度分布，t =1.8 ，N =401，CCFL=0.6Fig.7 Density profiles of the shock entropy interacting of Shu-Osher[5] at t =1.8 with 401 points and CCFL=0.6

3.3 二维 Euler 系统

在本小节中，考虑二维可压缩Euler 系统：

其中

这里 ρ ,u,v,p和E分别是密度，x和y坐 标方向上的速度分量，压力和总能量.另外，U是守恒变量的向量，F(U)是x方向的通量分量，G(U)是y方向的通量分量.对二维可压缩Euler 方程进行了逐维求解.

例3 (Rayleigh-Taylor 不稳定)当加速度从较重流体定向到较轻流体时，Rayleigh-Taylor 不稳定性发生在两种密度流体之间的界面上.这一问题在文献中得到了广泛的模拟[19-20].计算域为[ 0,1/4]×[0,1]，初始条件为

比热比γ =5/3.引力效应是通过将源项S=(0,0,ρ,ρv)T添加到二维Euler 方程(37)的右侧来引入的.左右边界为反射边界条件，顶部和底部边界设置为 (ρ,u,v,p)=(1,0,0,2.5)和 (ρ,u,v,p)=(2,0,0,1)，CFL 数为0.6，时间为t=1.95.图8 显示WENO-MS-Z3 格式比WENO-JS3、WENO-Z3 和WENO-MS-JS3 格式结构更为丰富，尤其是在界面上产生的涡流较多，表明其耗散性比其他格式小.

图8 不同格式关于Rayleigh-Taylor 不稳定性问题的密度等值线，Δ x=Δy=1/240，C CFL=0.6，t=1.95Fig.8 Density contours of the Rayleigh-Taylor instability at t =1.95 with Δ x=Δy=1/240，CCFL=0.6

例4 (双Mach 反射问题)[19]激波在斜面上的二维双Mach 反射描述了平面Mach 激波在空气中撞击楔形物的反射.这种试验被广泛用于验证数值方法的性能.我们在[ 0,4]×[0,1]上计算这个测试问题，并像往常一样在[0,3]×[0,1]上 显示结果.一开始施加一个右移Mach 数为10 的激波在x=1/6处 ，激波前缘与x轴成 60◦角.在沿底部边界x=0至x=1/6的区域内，为初始激波后流动赋值，其余区域采用反射边界条件.左边界为初始激波后流动赋值.对于x=4处的右侧边界都设置为零.上边界是用来描述Mach 数为10 的激波的精确运动.有关此问题的详细描述，请参见文献[19].我们采用的网格尺度为 Δx=Δy=1/480， 比热比 γ =1.4，CCFL=0.6，计算时间t=0.2.图9 比较了WENO-MS 格式与WENO-JS3 和WENO-Z3 格式的数值结果.可见，WENO-MS-JS3 格式的分辨率与WENO-Z3 格式大致相当，WENO-MS-Z3 格式分辨率最高，WENO-MS-Z3 较好地解决了Mach杆附近的不稳定性问题.

图9 双Mach 反射问题在t =0.2时 的密度等值线，(网格点为1 920×480)Fig.9 Density contours of the double Mach reflection problem at t =0.2 with 1 920×480 grid points

4 结 论

为了减少传统的三阶加权本质无振荡格式的数值耗散以及提高经典三阶WENO 格式的数值精度，本文发展了非线性双曲守恒律数值解WENO 三阶WENO-MS 格式的一种新的修正模板逼近方法.通过改进经典WENO-JS3 格式中各候选模板上数值通量的一阶多项式逼近，并加入二次项使模板逼近达到三阶精度.计算出了相应的候选通量.它可以新发展的WENO-MS 格式在光滑区域(包括一阶临界点)实现三阶收敛.在不影响计算精度的前提下，为了使新的修正模板WENO-MS 格式具有ENO 特性，我们在候选模板通量的修正项中加入可调函数 φ，并且给出了一系列一维和二维数值算例，验证了新格式的有效性.数值结果表明，与WENOJS3、WENO-Z3 格式相比，WENO-MS 格式对精细光滑结构具有相当或更高的分辨率.