1. 引言
双曲守恒律方程是流体力学的基本控制方程之一,应用于描述流体运动和行为。由于该方程存在间断解,即使初值条件足够光滑,依然会产生激波和间断现象。因此,为了有效地捕捉这些间断现象,并减少其带来的数值不稳定性,学者们展开了众多数值格式研究。其中,加权本质无震荡(Weighted Essentially Non-Oscillatory, WENO)格式的提出,推动了对包含激波等复杂流动结构的流场进行精确数值模拟的进程。
自1987年Harden和Osher在文献[1]中首次提出了ENO格式以来,1994年Liu等在此基础上通过非线性加权改进了模板,在文献[2]中提出WENO格式。1996年,Jiang和Shu对WENO格式进行了理论分析,推广至任意阶有限差分形式,并在文献[3]中提出了光滑因子和非线性权重的计算方法,从而得到了WENO-JS格式,该格式一经提出,便成为了经典WENO格式之一,许多学者在此基础上提出了新的WENO格式。2005年Henrick等对光滑因子进行映射,在文献[4]中提出了WENO-M格式,该格式确保了在临界点附近精度不降低。然而权重计算公式复杂,计算量较大。为此,学者们在文献[5]-[7]中通过映射函数来提高收敛精度,提出了多种新格式。2008年Borges等对子模板的低阶光滑因子进行线性组合,在文献[8]中提出了全局光滑因子概念,得到了WENO-Z格式,该格式在不增加计算量的前提下,能够达到与WENO-M格式相同的精度,随后,在文献[9]-[12]中学者们提出了一系列光滑因子的改进方法,但仍存在较大的数值耗散。为了减少耗散,Fan等在文献[13]中基于拉格朗日插值多项式,提出了WENO-η格式,而Arker等在文献[14]中通过增加相对不光滑的子模板权重,提出WENO-Z + 格式。Luo等在文献[15]中进一步提出WENO-Z + I格式,该格式不仅保证获得高于WENO-Z格式的分辨率,在临界点处也能达到最佳收敛精度。
尽管上述研究者们的格式改进提升了WENO格式的性能,然而如文献[16]-[19]中对非线性权重计算方法的改进与模板优化之间关系的研究仍然较为有限。本文在五阶WENO-Z格式的子模板基础上,通过改进其三点子模板,得到四个低阶的两点子模板,并对子模板进行重新线性组合,得到新的全局光滑因子。理论分析中,新因子在一阶临界点保持了五阶精度。数值实验结果表明,改进的格式相较于其他经典格式具有更低的耗散和更高的对流场结构分辨率。
2. 基本的WENO格式
2.1. 模型方程
文献[1]中,一维标量守恒定律可以写成:
(1.1)
其中
是守恒变量,
是通量函数。考虑一维的均匀网格
,等分成
个网格,网格单元
,
,
。其中定义
,
,且
,
,将其半离散化为:
(1.2)
其中
为
的数值近似,存在通量函数
,通过隐式定义:
(1.3)
由上得到(1.2)式的守恒性质:
(1.4)
可以得到关于
的四次多项式逼近函数
,
(1.5)
2.2. WENO-Z格式
在文献[8]中,给定五点模版
,五阶WENO-Z格式三个数值通量模板分别为
,
和
,
的数值通量为
(1.6)
对应子模板上各数值通量具体形式如下
(1.7)
式(1.6)中,非线性权重计算公式为
(1.8)
其中
为理想权重,
,
,
,
为光滑因子对非线性权重影响的指数,本文中选取
,
为防止分母为零的正数,一般取
,局部光滑因子
通用计算公式为
(1.9)
对于五阶WENO-JS格式,
表示为
(1.10)
Borges给出的五阶WENO-Z格式的非线性权为
(1.11)
其中
为全局光滑因子。
3. 改进的WENO格式
3.1. 模板重组
如图1所示,
上部分为WENO-Z格式的重构模版,下部分为改进WENO-Z格式的重构模版,记作WENO-NEW格式。
Figure 1. Reconstruction stencils of fifth-order WENO-NEW
图1. 五阶WENO-NEW格式重构模板
新格式子模板四个光滑因子计算公式为
(2.1)
通过线性组合得到新的全局光滑因子如下:
(2.2)
更新后的非线性权重为
(2.3)
相较于WENO-Z格式,新格式通过计算相邻点之间的差异性,提供了对流场局部光滑性更准确的评估,更敏感的反映流动特征的变化,既降低了激波或间断处数值耗散,也提高了光滑区域解的稳定性。
3.2. 精度分析
Henrick等人在文献[5]中给出了由光滑因子判断WENO格式五阶收敛的充分条件,即每个模板上的数值通量,根据泰勒展开式可表示为
(2.4)
其中
为与
无关的常数。Borges等人在文献[8]中给出了关于全局光滑因子
收敛性的详细证明过程与分析,其将(1.10)式在
处泰勒展开至五阶,由
的定义计算光滑因子,得到
(2.5)
将新定义模板下的四个光滑因子在
处泰勒展开至五阶,由
定义进行计算,可以得到:
(2.6)
中项的权重组合表明,光滑因子计算更依赖于二阶与三阶导数的乘积,而新格式不仅在一阶临界点保持了五阶精度,并且提供了更平衡的权重分配,处理高阶导数时,在光滑区域的计算更加简化,计算更为稳定,在非光滑区域,能够更好地控制数值震荡。
3.3. 频谱特性分析
利用Pirozzoli在文献[20]中提出的近似色散关系(ADR)。可以对非线性格式进行频谱特性分析。ADR分析通过求解模型方程。对给定网格上的全部Fourier模式进行求解。对于波数为
的正弦波,从数值结果换算得到的修正波数为
,其是一个复数。其实部反映格式的色散性质,虚部则反映格式的耗散性质。
从图2可知WENO-NEW格式频谱特性明显优于WENO-Z格式,证实了改进格式数值性能优越性,同时实验结果也证实了WENO-NEW格式有更小的耗散性,则格式为非光滑模板分配了更大的权重其频谱特性更好。另一个可以注意到WENO-NEW格式的耗散性总体上几乎等同于线性迎风格式,甚至在极小的波数范围内小于线性迎风格式。
(a)
(b)
Figure 2. Comparison of the spectral properties of different WENO schemes
图2. 不同WENO格式的频谱特性曲线比较
4. 数值实验
4.1. 线性对流方程问题
(3.1)
初始条件为
,表示一个平滑的正弦波,采用周期边界条件,计算终止时间为
,时间离散方法使用三阶TVD-RK方法。表1与表2分别表示五阶WENO-JS格式、WENO-M格式、WENO-Z格式与WENO-NEW格式的
和
误差与相应收敛阶的比较。
由表1与表2数据可知,WENO-NEW格式能实现最优收敛精度阶,计算误差在相同网格数下与WENO-Z格式、WENO-M格式保持同一量级,都较WENO-JS格式小一个量级,且新格式的截断误差数值更小。
Table 1. Difference WENO scheme L1 error and convergence order of the linear transport equation at t = 2 s
表1. t = 2 s时线性对流方程不同WENO格式的L1误差与收敛阶
网格数 |
WENO-JS |
WENO-M |
WENO-Z |
WENO-NEW |
Error |
Order |
Error |
Order |
Error |
Order |
Error |
Order |
20 |
|
|
|
|
|
|
|
|
40 |
|
4.97 |
|
5.03 |
|
5.07 |
|
5.01 |
80 |
|
5.00 |
|
5.00 |
|
5.05 |
|
5.00 |
160 |
|
5.01 |
|
5.00 |
|
5.00 |
|
5.00 |
320 |
|
5.04 |
|
5.00 |
|
5.00 |
|
5.00 |
Table 2. Difference WENO scheme
error and convergence order of the linear transport equation at t = 2 s
表2. t = 2 s时线性对流方程不同WENO格式的
误差与收敛阶
网格数 |
WENO-JS |
WENO-M |
WENO-Z |
WENO-NEW |
Error |
Order |
Error |
Order |
Error |
Order |
Error |
Order |
20 |
|
|
|
|
|
|
|
|
40 |
|
4.71 |
|
4.98 |
|
5.03 |
|
4.99 |
80 |
|
4.94 |
|
4.99 |
|
5.01 |
|
5.00 |
160 |
|
5.01 |
|
5.00 |
|
5.00 |
|
5.00 |
320 |
|
5.15 |
|
5.00 |
|
5.00 |
|
5.00 |
4.2. 一维波动方程
(3.2)
表示波动的状态量,
是通量函数,
是源项,表示外部作用力或影响,此处为0。
4.2.1. 线性问题
其中
,表示波的传播速度,采用周期边界条件,初始条件为高斯波,
。
(3.3)
表示波的中心位置,此处设为0,网格数
,计算终止时间为
。
Figure 3. Numerical solutions of the linear problem
图3. 线性问题的数值解
4.2.2. Buckley-Leverett问题
采用周期边界条件,初始条件为高斯波
(3.3)
其中
是扩散系数,计算终止时间为
,网格数
,
。
Figure 4. Numerical solutions of the Buckley-Leverett problem
图4. Buckley-Leverett问题的数值解
图3和图4是上述两种问题下各类型格式的数值计算结果,上述问题用于研究波动的传播与扩散现象,如图3、图4所示,WENO-NEW格式在间断处更好的逼近精确解,在间断处能够更加有效地保持解的稳定性和准确性。
4.3. 一维欧拉方程组问题
(3.4)
其中守恒变量
,通量函数
,一维欧拉方程组问题可以研究冲击波、激波与稀疏波等基本流动现象,适合用来验证和测试数值方法的准确性与稳定性。通量分裂使用Lax-Friedrichs通量分裂方法,网格单元变量的计算使用Roe平均方法,边界条件为Neumann周期边界条件。
4.3.1. Sod激波管问题
初始条件为
(3.5)
式(3.5)可表示为Sod激波管问题。图5为各类格式的密度计算结果及激波与接触间断附近的分布放大图,其中计算终止时间为
,
,网格数
,可以观察出,WENO-NEW格式数值耗散较小,计算效果更好。
Figure 5. Numerical density solutions of Sod shock tube problem
图5. Sod激波管问题的数值密度结果
4.3.2. Lax激波管问题
初始条件为
(3.6)
此时式(3.6)可表示为Lax激波管问题。图6为各类数值格式密度计算结果及激波与接触间断附近的分布放大图,其中计算终止时间为
,
,
,网格数
,图中可以看出,各类数值格式虽都在接触间断附近产生数值震荡,但WENO-NEW格式计算结果更靠近精确解。
Figure 6. Numerical density solutions of Lax shock tube problem
图6. Lax激波管问题的数值密度结果
4.4. 二维欧拉方程组问题
(3.7)
其中守恒变量
,
,
是关于
和
的通量函数。二维欧拉方程组问题能够验证处理复杂流动特征时的能力,用于比较格式的分辨率与解的稳定性。
4.4.1. 二维Riemann问题
计算区域为
,该区域被
与
划分成四个区域。
初值条件如下
(3.8)
此时,式(3.8)可表示为Riemann问题。各类WENO格式计算Riemann问题的密度等值线分布图如图7所示,计算网格选用
的均匀网格,边界使用Neumann边界条件,
,计算终止时间为
,共绘制30条等值线。
如图7中(a)~(c)显示,相较于其他WENO格式的数值结果,WENO-NEW格式对滑移线附近的流场微小变化结构捕捉最为精确,有效解析流场中的细微结构,局部适应性更强,具有更高的分辨率和精确度。
(a)
(b)
(c)
WENO的数值结果:(a) WENO-JS;(b) WENO-Z;(c) WENO-NEW。
Figure 7. Density contours of the 2D Riemann problem computed with the WENO-JS, WENO-Z and WENO-NEW schemes
图7. 二维Riemann问题不同格式(WENO-JS、WENO-Z、WENO-NEW)密度曲线图
4.4.2. Double Mach Reflection (DMR)问题
初始条件如下
(3.9)
上述问题描述为马赫数为10的激波从以倾角为30˚向右移动
处以倾角为30˚向右移动,
的区间为
,计算网格为
的均匀网格,计算终止时间为
,
右侧使用反射边界条件,绘制30条密度等值线,数值模拟结果如图8所示。
(a)
(b)
(c)
WENO格式的数值结果:(a) WENO-JS,CPU time:6576.9221 s;(b) WENO-Z,CPU time:6659.405 s;(c) WENO-NEW,CPU time:6711.5042 s。
Figure 8. Density contours of the DMR problem computed with the WENO-JS, WENO-Z and WENO-NEW schemes
图8. DMR问题不同格式(WENO-JS、WENO-Z、WENO-NEW)密度曲线图
图中可以观察出,WENO-NEW格式在马赫杆区域描绘的漩涡结构更加丰富,滑移线卷起的涡流特征也更明显,耗散更低。这表明新格式具有更低的数值耗散和更高的空间分辨率,并且WENO-NEW格式计算时间仅比WENO-Z格式增加0.7%。
5. 结论
本文在五阶WENO-Z格式的基础上提出一种新的模板重组方法,得到改进全局光滑因子,理论分析表明,WENO-NEW格式在一阶临界点处能够达到最佳计算精度,且对导数的权重分配更加平衡,充分利用了局部光滑性信息,提高了数值解的准确性。通过线性对流方程的算例研究,新格式展现出更小的数值误差以及更优的精度稳定性。此外,基于一维和二维欧拉方程组的算例比较,数值实验结果表明,WENO-NEW格式不仅降低了间断处的数值耗散、分辨率更高,还在相同精度下,提高了对复杂流场细微变化结构的捕捉能力。
基金项目
山西省基础研究计划,项目编号202403021221030。
NOTES
*通讯作者。