数学杂志  2024, Vol. 44 Issue (5): 441-452   PDF    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
高志辉
张学莹
基于混合模板改进的五阶WENO格式在求解双曲守恒律方程中的应用
高志辉, 张学莹    
河海大学数学学院, 江苏 南京 211100
摘要:本文在WENO-MZQ3格式的基础上, 将其模板推广到五阶WENO, 即在经典五阶WENO格式的三个三点模板的基础上引入五点模板, 并引入新的改进的全局光滑因子, 对其进行加权处理得到新的五阶WENO格式. 理论分析表明新的五阶WENO格式在光滑区域比五阶WENO-JS格式具有更高阶的计算精度, 通过一系列一维和二维数值算例也验证了新的五阶WENO格式具有更高的精度和分辨率.
关键词WENO格式    五点模板    全局光滑因子    
APPLICATION OF THE IMPROVED FIFTH ORDER WENO SCHEMES BASED ON A HYBRID TEMPLATE IN SOLVING HYPERBOLIC CONSERVATION EQUATION
GAO Zhi-hui, ZHANG Xue-ying    
School of Mathematics, Hohai University, Nanjing 211100, China
Abstract: In this paper, we extend the template of WENO-MZQ3 scheme to fifth-order WENO schemes. This means that we introduce a five-point template based on the three three-point templates of the conventional fifth-order WENO scheme. We introduce a new improved global smoothness indicator and apply weighted treatment to obtain a new fifth-order WENO scheme. Theoretical analysis demonstrate that the new fifth-order WENO scheme achieves higher computational accuracy in smooth regions compared to the fifth-order WENO-JS scheme. A series of one-dimensional and two-dimensional numerical examples validate the higher accuracy and resolution of the new fifth-order WENO scheme.
Keywords: WENO scheme     five-point template     global smoothness indicator    
1 引言

随着计算机的不断发展, 数值模拟是求解计算流体力学方程的重要手段之一. 双曲守恒律方程是流体力学的基本控制方程之一. 在特定参数下, 初始条件即便为光滑也可能产生激波和接触间断等不连续解, 为了消除间断处的数值振荡, 许多学者建立了一系列理论与方法. 1959年Goduov提出迎风格式的思想[1]. 1983年Harten建立了TVD格式[2]. 1987年Harten和Osher放宽了TVD条件[3], 提出了ENO格式. 1994年, Liu等采用加权思想[4], 将ENO候选模板以凸组合形式进行非线性加权, 得到WENO格式. 1996年, Jiang和Shu将WENO格式推广到有限差分中[5], 并定义新的光滑因子, 提出了WENO-JS格式. 许多学者在WENO-JS的基础上提出新的WENO格式, 如通过映射函数来提高WENO格式的收敛精度[69, 13], 通过光滑因子的修改来提高WENO格式的精度[1012]等. Anurag Kumar等在经典三阶WENO模板中引入三点模板[14], 在此模板基础上进行重构, 并提出新的全局光滑因子, 进而得到一种改进的三阶WENO格式, 称为WENO-MZQ3. 本文在WENO-MZQ3格式的基础上, 将其推广到五阶WENO格式, 即在经典五阶WENO格式的三个三点模板的基础上引入五点模板, 通过对各模板相应的光滑因子的数值通量进行泰勒展开, 提出新的全局光滑因子, 对其进行加权处理得到新的五阶WENO格式. 通过理论分析可知新的五阶WENO格式在光滑区域相比五阶WENO-JS具有更高的精度. 通过数值算例验证, 新格式具有很好的分辨率和精度.

2 数值方法
2.1 五阶WENO-JS格式回顾[5]

在本节中, 考虑如下一维标量守恒律方程

$ \begin{equation} {u_t}(x, t) + {f_x}(u(x, t)) = 0. \end{equation} $ (2.1)

(2.1)式中, $ {u}(x, t) $表示空间$ {x} $和时间$ {t} $的守恒变量, $ {f}(u) $为通量函数. (2.1)式中$ x $计算区间为[a, b], 令$ {x_i}=a+i{\Delta x} $, 则单元格大小为$ {{\Delta x}=(b-a)/N} $. (2.1)式在$ {x_i} $处的半离散形式如下

$ \begin{equation} \frac{{d{u_i}}}{{dt}} = - \frac{1}{{\Delta x}}\left( {{{\hat f}_{i + \frac{1}{2}}} - {{\hat f}_{i - \frac{1}{2}}}} \right), \ \end{equation} $ (2.2)

其中$ {\hat f_{i + \frac{1}{2}}} = \hat f_{i + \frac{1}{2}}^ + + \hat f_{i + \frac{1}{2}}^ - $为中心节点$ {{x_{i+\frac{1}{2}}}={x_i}+\frac{\Delta x}{2}(i=1, 2, \cdots, N.)} $处的数值通量, $ {u_i} $$ {{u}(x, t)} $在节点$ {x_i} $处的近似值.单元格界面处的数值通量$ \hat f_{i + \frac{1}{2}}^ \pm $需要通过通量分裂法进行重构.本文中一律采取全局Lax-Friedrichs数值通量, $ {\hat f^ \pm }(u) = \frac{1}{2}(\hat f(u) \pm \alpha u) $, 其中$ \alpha = \max (\left| {f'(u)} \right|) $.在本文中只求$ \hat f_{i + \frac{1}{2}}^ + $, $ \hat f_{i + \frac{1}{2}}^ - $可以由对称性得到.为了方便起见, 在本文的以下部分中, 忽略上标中的$ {\pm} $.

数值通量函数$ {\hat f_{i + \frac{1}{2}}} $在五点模板上的五阶有限差分WENO格式[5]定义如下:

$ \begin{equation} {\hat f_{i + \frac{1}{2}}} = {\omega _0}{p_0} + {\omega _1}{p_1} + {\omega _2}{p_2}, \end{equation} $ (2.3)

其中

$ \begin{equation} \left\{ \begin{array}{l} {p_0} = \left( {2{f_{i - 2}} - 7{f_{i - 1}} + 11{f_i}} \right)/6, \\ {p_1} = \left( { - {f_{i - 1}} + 5{f_i} + 2{f_{i + 1}}} \right)/6, \\ {p_2} = \left( {2{f_i} + 5{f_{i + 1}} - {f_{i + 2}}} \right)/6. \end{array} \right.\ \end{equation} $ (2.4)

WENO格式的数值性能由非线性权$ {\omega_k(k=0, 1, 2.)} $决定.下面给出经典五阶WENO-JS格式[5]的非线性权重:

$ \begin{equation} \omega _k^{JS} = \frac{{\alpha _k^{JS}}}{{\sum\nolimits_{j = 0}^2 {\alpha _j^{JS}} }}, \;\;\alpha _k^{JS} = \frac{{{d_k}}}{{{{\left( {\varepsilon + {\beta _k}} \right)}^q}}}.\qquad {\rm{ k}} = 0, 1, 2. \end{equation} $ (2.5)

其中$ d_0=0.1, d_1=0.6, d_2=0.3 $是线性权值, 设置$ {\varepsilon > 0} $为了避免分母为零, $ q $为可调正整数.$ {\beta_k(k=0, 1, 2.)} $是每个子模版的光滑因子, 表达式如下:

$ \begin{equation} {\beta _0} = \frac{{13}}{{12}}{\left( {{f_{i - 2}} - 2{f_{i - 1}} + {f_i}} \right)^2} + \frac{1}{4}{\left( {{f_{i - 2}} - 4{f_{i - 1}} + 3{f_i}} \right)^2;} \end{equation} $ (2.6)
$ \begin{equation} {\beta _1} = \frac{{13}}{{12}}{\left( {{f_{i - 1}} - 2{f_{i}} + {f_{i + 1}}} \right)^2} + \frac{1}{4}{\left( {{f_{i - 1}} - {f_{i+1}}} \right)^2;} \end{equation} $ (2.7)
$ \begin{equation} {\beta _2} = \frac{{13}}{{12}}{\left( {{f_i} - 2{f_{i + 1}} + {f_{i + 2}}} \right)^2} + \frac{1}{4}{\left( {3{f_i} - 4{f_{i + 1}} + {f_{i + 2}}} \right)^2.} \end{equation} $ (2.8)
2.2 Embedded WENO-JS格式回顾[15]

Embedded WENO-JS格式是在WENO-JS的基础上提出的.模板与WENO-JS模板一致, 其各个子模版对应的光滑因子为(2.6)、(2.7)、(2.8).Embedded WENO-JS格式贡献在于提出了新的非线性权重$ \alpha_k^{E} $, 通过修改非线性权重的值来提高格式的精度.对应的非线性权重$ \alpha_k^{E} $为[15, 16]:

$ \begin{equation} \ \left\{ \begin{array}{l} \alpha _0^{E} = \frac{1}{3}{d_0}(3 - {c_2} + {c_2}\frac{{{\beta _2}}}{{{\beta _0} + \varepsilon }}), \\ \alpha _1^{E} = \frac{1}{3}{d_1}(1 + \frac{{{\beta _2}}}{{{\beta _1} + \varepsilon }} + \frac{{{\beta _0}}}{{{\beta _1} + \varepsilon }}), \\ \alpha _2^{E} = \frac{1}{3}{d_2}(3 - {c_0} + {c_0}\frac{{{\beta _0}}}{{{\beta _2} + \varepsilon }}). \end{array} \right.\ \end{equation} $ (2.9)

其中$ c_0=\frac{2}{3}, c_2=\frac{6}{7}. $

$ \begin{equation} \ \omega _k^{E} = \frac{{\alpha _k^{E}}}{{\sum\nolimits_{j = 0}^2 {\alpha _j^{E}} }}.\qquad {\rm{ k}} = 0, 1, 2. \end{equation} $ (2.10)

单元界面$ {x=x_{i+\frac{1}{2}}} $处的数值通量为:

$ \begin{equation} \ {\hat f_{i + \frac{1}{2}}} = {\omega _0^{E}}{p_0} + {\omega _1^{E}}{p_1} + {\omega _2^{E}}{p_2}. \end{equation} $ (2.11)
2.3 混合模板的五阶WENO格式

在本节中, 将介绍新型的五阶WENO格式.

在五阶WENO格式模板的基础上增添一个中心大模板, 从而构造出新的WENO格式, 即WENO-HT格式.其对应的模板为$ {{S_0}=\left\{ {{x_{i - 2}}, {x_{i - 1}}, {x_i}} \right\};} $ $ {{S_1}=\left\{ {{x_{i - 1}}, {x_{i}}, {x_{i + 1}}} \right\};} $ $ {{S_2}=\left\{ {{x_{i}}, {x_{i + 1}}, {x_{i + 2}}} \right\};} $ $ {{S_5}=\left\{ {{x_{i - 2}}, {x_{i - 1}}, {x_i}, {x_{i + 1}}, {x_{i + 2}}} \right\}.} $四个模板在$ {x_{i+\frac{1}{2}}} $对应的数值通量为:

$ \begin{equation} \left\{ \begin{array}{l} {p_0} = \left( {2{f_{i - 2}} - 7{f_{i - 1}} + 11{f_i}} \right)/6, \\ {p_1} = \left( { - {f_{i - 1}} + 5{f_i} + 2{f_{i + 1}}} \right)/6, \\ {p_2} = \left( {2{f_i} + 5{f_{i + 1}} - {f_{i + 2}}} \right)/6, \\ {p_5} = \left( {2{f_{i-2}} - 13{f_{i-1}}+ 47{f_i} + 27{f_{i + 1}} - 3{f_{i + 2}}} \right)/60. \end{array} \right.\ \end{equation} $ (2.12)

其中$ {{p_0}, {p_1}, {p_2}} $对应的光滑因子由(2.6)、(2.7)、(2.8)式给出, 接下来介绍$ {p_5} $相对应的光滑因子[12]:

$ \begin{equation} \begin{aligned} {\beta _5} &= \frac{1}{{144}}{\left( {{f_{i - 2}} - 8{f_{i - 1}} + 8{f_{i + 1}} - {f_{i + 2}}} \right)^2}\\ &+\frac{1}{{15600}}{\left( { - 11{f_{i - 2}} + 174{f_{i - 1}} - 362{f_i} + 174{f_{i + 1}} - 11{f_{i + 2}}} \right)^2}\\ &+\frac{{781}}{{2880}}{\left( { - {f_{i - 2}} + 2{f_{i - 1}} - 2{f_{i + 1}} + {f_{i + 2}}} \right)^2}\\ &+\frac{{1421461}}{{1310400}}{\left( {{f_{i - 2}} - 4{f_{i - 1}} + 6{f_i} - 4{f_{i + 1}} + {f_{i + 2}}} \right)^2}. \end{aligned} \end{equation} $ (2.13)

对光滑因子(2.6)、(2.7)、(2.8)、(2.13)式在$ {x={x_i}} $进行泰勒级数展开可得到[12, 13]:

$ \begin{equation} \begin{aligned} {\beta _0} &= {({{f_i}'})^2}{(\Delta x)^2} - \left( {\frac{2}{3}{{f_i}'}f_i^{(3)} - \frac{{13}}{{12}}{{({{f_i}''})}^2}} \right){(\Delta x)^4} + \left( {\frac{1}{2}{{f_i}'}f_i^{(4)} - \frac{{13}}{6}{{f_i}''}f_i^{(3)}} \right){(\Delta x)^5}\\ &- \left( {\frac{7}{{30}}{{f_i}'}f_i^{(5)} - \frac{{91}}{{72}}{{f_i}''}f_i^{(4)} - \frac{{43}}{{36}}{{(f_i^{(3)})}^2}} \right){(\Delta x)^6} \\ &+ \left( {\frac{1}{{12}}{{f_i}'}f_i^{(6)} - \frac{{13}}{{24}}{{f_i}''}f_i^{(5)} - \frac{{103}}{{72}}f_i^{(3)}f_i^{(4)}} \right){(\Delta x)^7} + {\rm O}({(\Delta x)^8}); \end{aligned} \end{equation} $ (2.14)
$ \begin{equation} \begin{aligned} {\beta _1} &= {({{f_i}'})^2}{(\Delta x)^2} + \left( {\frac{1}{3}{{{f_i}'}f_i^{(3)}} + \frac{{13}}{{12}}{{({{f_i}''})}^2}} \right){(\Delta x)^4} + \left( {\frac{1}{{60}}{{f_i}'}f_i^{(5)} + \frac{{13}}{{72}}{{f_i}''}f_i^{(4)}} \right){(\Delta x)^6}\\ &+\left( {\frac{{1}}{{36}}{{(f_i^{(3)})}^2}}\right){(\Delta x)^6}+ {\rm O}({(\Delta x)^8}); \end{aligned} \end{equation} $ (2.15)
$ \begin{equation} \begin{aligned} {\beta _2} &= {({{f_i}'})^2}{(\Delta x)^2} - \left( {\frac{2}{3}{{{f_i}'}f_i^{(3)}} - \frac{{13}}{{12}}{{({{f_i}''})}^2}} \right){(\Delta x)^4} - \left( {\frac{1}{2}{{f_i}'}f_i^{(4)} - \frac{{13}}{6}{{f_i}''}f_i^{(3)}} \right){(\Delta x)^5}\\ &- \left( {\frac{7}{{30}}{{f_i}'}f_i^{(5)} - \frac{{91}}{{72}}{{f_i}''}f_i^{(4)} - \frac{{43}}{{36}}{{(f_i^{(3)})}^2}} \right){(\Delta x)^6} \\ &- \left( {\frac{1}{{12}}{{f_i}'}f_i^{(6)} - \frac{{13}}{{24}}{{f_i}''}f_i^{(5)} - \frac{{103}}{{72}}f_i^{(3)}f_i^{(4)}} \right){(\Delta x)^7} + {\rm O}({(\Delta x)^8}); \end{aligned} \end{equation} $ (2.16)
$ \begin{equation} \begin{aligned} {\beta _5} &= {({f_i}')^2}{(\Delta x)^2} + \frac{{13}}{{12}}{({f_i}'')^2}{(\Delta x)^4} + \frac{1}{{5040}}\left( { - 336{{f_i}'}f_i^{(5)} + 5467{{(f_i^{(3)})}^2}} \right){(\Delta x)^6}\\ &- \frac{1}{{360}}\left( {{f_i}''}f_i^{(4)}\right){(\Delta x)^6} + {\rm O} ({(\Delta x)^8}). \end{aligned} \end{equation} $ (2.17)

由(2.14)、(2.15)、(2.16)式可得:

$ \begin{equation} \begin{aligned} \beta &= \frac{1}{6}{({{4}\beta _1+\beta_0+\beta_2})} = {({f_i}')^2}{(\Delta x)^2} + \frac{{13}}{{12}}{({f_i}'')^2}{(\Delta x)^4}\\ &+ \left( {{-}\frac{1}{{15}}{{f_i}'}f_i^{(5)} + \frac{{13}}{{24}}{{f_i}''}f_i^{(4)} + \frac{5}{{12}}{( {{f_i}^{(3)}})}^{2}} \right){(\Delta x)^6} + {\rm O} ({(\Delta x)^8}). \end{aligned} \end{equation} $ (2.18)

由(2.17)、(2.18)式可得:

$ \beta_5 - \beta = \left( {- \frac{{49}}{{90}}{{f_i}''}f_i^{(4)} + \frac{{481}}{{720}}{{(f_i^{(3)})}^2}} \right){(\Delta x)^6}+ {\rm O} ({(\Delta x)^8}). $

由上式可知:

$ \begin{equation} {\beta _5} - \beta = \left\{ \begin{array}{l} {\rm O}({(\Delta x)^6}), \;\;{{f''}_i} \ne 0, f_i^{(3)} \ne 0, f_i^{(4)} \ne 0;{\rm{ }}\\ {\rm O}({(\Delta x)^6}), \;\;{{f''}_i} = 0, f_i^{(3)} \ne 0, f_i^{(4)} \ne 0;{\rm{ }}\\ {\rm O}({(\Delta x)^8}), \;\;{{f''}_i} = 0, f_i^{(3)} = 0{\rm{ }}. \end{array} \right.{\rm{ }}\ \end{equation} $ (2.19)

此时, 令$ {\tau _{HT}} = \left| {{\beta _5} - \beta } \right|. $定义以下非线性权:

$ \begin{equation} {\omega _k} = \frac{{{{\bar \alpha }_k}}}{{{{\bar \alpha }_0} + {{\bar \alpha }_1} + {{\bar \alpha }_2} + {{\bar \alpha }_5}}}, {\bar \alpha _k} = {\eta _k}(1 + \frac{{{\tau _{HT}}}}{{{\beta _k} + \varepsilon }}).\qquad {\rm{k}} = 0, 1, 2, 5. \end{equation} $ (2.20)

在(2.20)式中, $ {\varepsilon} $是很小的正数, 只是为了避免分母为零.单元界面$ {x=x_{i+\frac{1}{2}}} $处的数值通量为:

$ \begin{equation} {\hat f_{i + \frac{1}{2}}} = {\omega _5}\left( {\frac{1}{{{\eta _5}}}{p_5} - \frac{{{\eta _0}}}{{{\eta _5}}}{p_0} - \frac{{{\eta _1}}}{{{\eta _5}}}{p_1} - \frac{{{\eta _2}}}{{{\eta _5}}}{p_2}} \right) + {\omega _0}{p_0} + {\omega _1}{p_1} + {\omega _2}{p_2}. \end{equation} $ (2.21)

在(2.21)式中: $ {\eta _5, \eta _0, \eta _1, \eta _2} $为线性权, 其和为1, 即$ {\eta _5 + \eta _0 + \eta _1 + \eta _2 = 1.} $其组合使格式在光滑区域达到更高精度.

2.4 时间离散

文中时间推进采用三阶TVD Runge-Kutta方法.如下:

$ \begin{equation} \left\{ \begin{array}{l} {u^{(1)}} = {u^n} + \Delta tL({u^n});\\ {u^{(2)}} = \frac{3}{4}{u^n} + \frac{1}{4}{u^{(1)}} + \frac{1}{4}\Delta tL({u^{(1)}});\\ {u^{n + 1}} = \frac{1}{3}{u^n} + \frac{2}{3}{u^{(2)}} + \frac{2}{3}\Delta tL({u^{(2)}}). \end{array} \right. \end{equation} $ (2.22)

其中$ L(u_i)=- \frac{1}{{\Delta x}}\left( {{{\hat f}_{i + \frac{1}{2}}} - {{\hat f}_{i - \frac{1}{2}}}} \right). $在第三节中时间步长采用变时间步长, $ {\Delta t} $在一维和二维分别为:

$ \Delta t=\sigma \frac{\Delta x}{\max\limits_i(|u_i| + a_i)}. $
$ \Delta t=\sigma \frac{\Delta t_x\Delta t_y}{\Delta t_x+\Delta t_y}, \Delta t_x=\sigma \frac{\Delta x}{\max\limits_{i, j}(|u_{i, j}| + a_{i, j})}, \Delta t_y=\sigma \frac{\Delta y}{\max\limits_{i, j}(|v_{i, j}| + a_{i, j})}. $

其中$ {\sigma} $是CFL$ (Courant-Friedrichs-Lewy) $条件数, 本文算例中CFL条件数均取$ {0.5} $.

3 数值实验

在本节中, 将五阶WENO-HT格式和WENO-JS、Embedded WENO-JS格式用于求解一维和二维欧拉方程组, 从而来验证WENO-HT格式具有更高的精度和分辨率.为了检验$ \eta $的选择是否会影响WENO-HT格式的最优精度, 在(2.22) 式$ \eta_k(k=0, 1, 2, 5) $的选取按如下几组数值进行测试:

$ (1)\eta_0=0.01, \eta_1=0.01, \eta_2=0.01, \eta_5=0.97; $

$ (2)\eta_0=\frac{1}{222}, \eta_1=\frac{10}{111}, \eta_2=\frac{1}{222}, \eta_5=\frac{100}{111}; $

$ (3)\eta_0=\frac{1}{110}, \eta_1=\frac{6}{110}, \eta_2=\frac{3}{110}, \eta_5=\frac{100}{110}; $

$ (4)\eta_0=0.05, \eta_1=0.3, \eta_2=0.15, \eta_5=0.5; $

$ (5)\eta_0=0.25, \eta_1=0.25, \eta_2=0.25, \eta_5=0.25. $

其中(1)、(2)、(5)中$ \eta_k(k=0, 1, 2, 5) $的选取见参考文献[12][14], (3)、(4)中$ \eta_k(k=0, 1, 2, 5) $的选取是将三个小模版对应的$ \eta_k(k=0, 1, 2) $之和与大模板对应的$ \eta_5 $按如下比例设置: 在WENO-HT-3格式中$ (\eta_0+\eta_1+\eta_2):\eta_5=1:10 $; 在WENO-HT-4格式中$ (\eta_0+\eta_1+\eta_2):\eta_5=1:1 $, 小模版对应的$ \eta_k(k=0, 1, 2) $按(2.5)式中的$ d_k(k=0, 1, 2) $的比例设置, 即$ \eta_0:\eta_1:\eta_2=1:6:3 $.

3.1 Sod激波管问题

初始条件如下

$ \begin{equation} \begin{aligned} \left( {\rho , u, p} \right) = \left\{ \begin{array}{rcl} \left( {0.125, 0, 0.1} \right), &&{ 0 \leq x \leq 1}, \\ \left( {1, 0, 1} \right), &&{- 1 \leq x < 0}.\\ \end{array} \right. \end{aligned} \end{equation} $ (3.1)

表 3.1给出了WENO-JS、Embedded WENO-JS与WENO-HT-1—WENO-HT-5格式的$ L^1 $误差.由表 3.1知, WENO-HT格式相比于WENO-JS, Embedded WENO-JS格式具有更小的误差.图 3.1为WENO-JS, Embedded WENO-JS, WENO-HT格式的模拟结果, 其中网格取400个均匀网格, 时间计算到t=0.14. 从图 3.1中可以看出WENO-HT-2出现震荡, WENO-HT-3与WENO-HT-4更加贴近近似Riemann解, 因此本文将WENO-JS, Embedded WENO-JS, WENO-HT-3, WENO-HT-4与近似Riemann解相比, 如图 3.1(b)所示.

表 3.1$ t=0.14$时(3.1)式WENO方案的$ L^1$误差

图 3.1 不同WENO格式对Sod激波管问题的模拟结果
3.2 Lax激波管问题

初始条件如下

$ \begin{equation} \begin{aligned} \left( {\rho , u, p} \right) = \left\{ \begin{array}{rcl} \left( {0.5, 0, 0.571} \right), &&{ 0 \leq x \leq 5}, \\ \left( {0.445, 0.698, 3.528} \right), &&{- 5 \leq x < 0}.\\ \end{array} \right. \end{aligned} \end{equation} $ (3.2)

表 3.2给出了WENO-JS、Embedded WENO-JS与WENO-HT-1—WENO-HT-5格式的$ L^1 $误差.由表 3.2知, WENO-HT格式相比于WENO-JS, Embedded WENO-JS格式具有更小的误差.图 3.2为WENO-JS, Embedded WENO-JS, WENO-HT格式对的模拟结果, 其中网格取1800个均匀网格, 时间计算到t=1.3. 从图 3.2中可以看出WENO-HT-1, WENO-HT-2出现震荡, WENO-HT-3与WENO-HT-4更加贴近近似Riemann解, 因此本文将WENO-JS, Embedded WENO-JS, WENO-HT-3, WENO-HT-4与近似Riemann解相比, 如图 3.2(b)所示.

表 3.2$ t=1.3$时(3.2)式WENO方案的$ L^1$误差

图 3.2 不同WENO格式对Lax激波管问题的模拟结果
3.3 Shu-Osher问题

初始条件如下

$ \begin{equation} \begin{aligned} \left( {\rho , u, p} \right) = \left\{ \begin{array}{rcl} \left( {1+0.2\sin{5x}, 0, 1} \right), &&{ -4 \leq x \leq 5}, \\ \left( {\frac{27}{7}, \frac{4\sqrt{35}}{9}, \frac{31}{3}} \right), &&{- 5 \leq x < -4}.\\ \end{array} \right. \end{aligned} \end{equation} $ (3.3)

图 3.3为WENO-JS, Embedded WENO-JS, WENO-HT格式对Shu-Osher问题的模拟结果, 其中网格取800个均匀网格, 时间计算到t=1.8. 其中参考解是网格取4000个均匀网格时WENO-JS的计算结果.从图 3.3可观察出WENO-HT-3与WENO-HT-4更加贴近参考解, 因此本文将WENO-JS, Embedded WENO-JS, WENO-HT-3, WENO-HT-4与参考解相比, 如图 3.3(b)所示.

图 3.3 不同WENO格式对Shu-Osher激波管问题的模拟结果
3.4 Titarev-Toro问题

初始条件如下

$ \begin{equation} \begin{aligned} \left( {\rho , u, p} \right) = \left\{ \begin{array}{rcl} \left( {1+0.1\sin{20\pi x}, 0, 1} \right), &&{ -4.5 \leq x \leq 5}, \\ \left( 1.515695, 0.523346, 1.805 \right), &&{- 5 \leq x < -4.5}.\\ \end{array} \right. \end{aligned} \end{equation} $ (3.4)

图 3.4为WENO-JS, Embedded WENO-JS, WENO-HT格式对Titarev-Toro问题的模拟结果, 其中网格取1500个均匀网格, 时间计算到t=5.其中参考解是网格取5000个均匀网格时WEWNO-JS的计算结果.从图 3.4中可以看出WENO-HT-1, WENO-HT-5出现强烈震荡, WENO-HT-3与WENO-HT-4更加贴近参考解, 因此本文将WENO-JS, Embedded WENO-JS, WENO-HT-3, WENO-HT-4与参考解相比, 如图 3.4(b)所示.

图 3.4 不同WENO格式对Titarev-Toro激波管问题的模拟结果

不同的$ {\eta _0}, {\eta _1}, {\eta _2}, {\eta _5} $的选取, 导致各模板的贡献不同.在光滑区域, 大模板$ S_5 $和中心模板$ S_1 $的贡献越大, 光滑区域的精度越高. 在间断区域, 需要各个模板选择合适的权系数来控制震荡. 因此需要调整大小模板之间的比例系数, 从而使小模版能够控制在间断处的震荡. 由表 3.1表 3.2图 3.1图 3.2图 3.3图 3.4知, $ \eta $的选取影响WENO-HT的精度和分辨率.通过实验, 比较图 3.1-3.4知, WENO-HT-3, WENO-HT-4与其它格式相比具有更高的分辨率.在图 3.2中看出WENO-HT-3出现一处略微明显的震荡.在接下来的二维欧拉方程组的数值模拟中选取WENO-HT-4.

3.5 Rayleigh-Taylor问题

初始条件如下

$ \begin{equation} \begin{aligned} \left( {\rho , u, v, p} \right) = \left\{ \begin{array}{rcl} \left( {2, 0, -0.025 c\cdot \cos{8\pi x}, 2y+1} \right), &&{ 0 \leq x \leq 0.25}, {0 \leq y < 0.5}, \\ \left( {1, 0, -0.025 c\cdot \cos{8\pi x}, y+1.5} \right), &&{ 0 \leq x \leq 0.25}, {0.5 \leq y \leq 1}.\\ \end{array} \right. \end{aligned} \end{equation} $ (3.5)

(3.5)式中, $ c=\sqrt{\frac{\gamma p}{\rho}} $, $ \gamma=\frac{5}{3} $.图 3.5(a), 图 3.5(b), 图 3.5(c)给出了网格为$ {480\times120} $的均匀网格, 时间计算到t=1.95的数值模拟结果.图 3.5(d), 图 3.5(e), 图 3.5(f)给出了网格为$ {960\times240} $的均匀网格, 时间计算到t=1.95的数值模拟结果.从图 3.5中可以看出WENO-HT的分辨率高于WENO-JS格式, 略高于EmbeddedWENO-JS格式, 表明WENO-HT的数值耗散低于WENO-JS, Embedded WENO-JS格式.

图 3.5 不同WENO格式对Rayleigh-Taylor问题的模拟结果
3.6 Double Mach Reflection问题

初始条件如下

$ \begin{equation} \begin{aligned} \left( {\rho , u, v, p} \right) = \left\{ \begin{array}{rcl} \left( {8, 4.125\sqrt{3}, -4.125, 116.5} \right), &&{x \leq \frac{1}{6} + \frac{y}{\sqrt{3}} }, \\ \left( {1.4, 0, 0, 1} \right), &&{ x > \frac{1}{6} + \frac{y}{\sqrt{3}}}.\\ \end{array} \right. \end{aligned} \end{equation} $ (3.6)

(3.6)式中, $ x, y $的区间为$ [0, 4]\times[0, 1] $, 比热比$ \gamma=1.4 $.图 3.6(a), 图 3.6(c), 图 3.6(e)给出了网格为$ {480\times120} $的均匀网格, 时间计算到t=0.2的数值模拟结果.图 3.6(b), 图 3.6(d), 图 3.6(f)给出了网格为$ {960\times240} $的均匀网格, 时间计算到t=0.2的数值模拟结果.相比于WENO-JS, Embedded WENO-JS格式WENO-HT能够观察到清晰的漩涡免卷结构, 本算例证明了WENO-HT数值耗散更低.

图 3.6 不同WENO格式对Double Mach Reflection问题的模拟结果
3.7 2D Riemann问题

初始条件如下

$ \begin{equation} \begin{aligned} \left( {\rho , u, v, p} \right) = \left\{ \begin{array}{rcl} \left( {0.13, 1.206, 1.206, 0.029} \right), &&{0 \leq x < 0.8}, {0 \leq y < 0.8}, \\ \left( {1.5, 0, 0, 1.5} \right), &&{ 0.8 \leq x \leq 1}, {0.8 \leq y \leq 1}, \\ \left( {0.5323, 1.206, 0, 0.3} \right), &&{0 \leq x < 0.8}, {0.8 \leq y \leq 1}, \\ \left( {0.5323, 0, 1.206, 0.3} \right), &&{0.8 \leq x \leq 1}, {0 \leq y < 0.8}.\\ \end{array} \right. \end{aligned} \end{equation} $ (3.7)

(3.7)式中比热比$ \gamma=1.4 $. 图 3.7(a), 图 3.7(b), 图 3.7(c)给出了网格为$ {200\times200} $的均匀网格, 时间计算到t=0.8的数值模拟结果.图 3.7(d), 图 3.7(e), 图 3.7(f)给出了网格为$ {400\times400} $的均匀网格, 时间计算到t=0.8的数值模拟结果.从图 3.7中可以看出WENO-HT的分辨率高于WENO-JS格式, 略高于Embedded WENO-JS格式, 本算例证明了WENO-HT数值耗散更低.

图 3.7 不同WENO格式对2D Riemann问题的模拟结果

以上二维算例的结果验证了WENO-HT格式在模拟复杂问题时与WENO-JS, Embedded WENO-JS格式相比具有更高的精度和分辨率.

4 结论

本文在WENO-MZQ3格式的基础上, 将其模板推广到五阶WENO, 即在经典五阶WENO格式的三个三点模板的基础上引入五点模板, 并引入新的改进的全局光滑因子, 对其进行加权处理得到新的五阶WENO格式, 即WENO-HT格式.理论分析表明, 五阶WENO-HT格式在光滑区域比五阶WENO-JS具有更高阶的计算精度. 一维欧拉方程的数值实验结果表明线性权重$ \eta_k $的选取影响WENO-HT格式的精度和分辨率, 通过比较最终选取权系数为$ \eta_0=0.05, \eta_1=0.3, \eta_2=0.15, \eta_5=0.5 $的五阶WENO-HT格式. 二维欧拉方程的数值实验结果显示在复杂问题中WENO-HT格式仍具有较高的精度和分辨率. 通过一系列的数值实验证明了WENO-HT格式相比于WENO-JS格式具有更高的精度和分辨率.

参考文献
[1] Godunov S K. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[J]. Mathematicheskii Sbornik, 1959, 47(3): 271–306.
[2] Harten A. High resolution schemes for hyperbolic conservation laws[J]. Journal of Computational Physics, 1983, 49(3): 357–393. DOI:10.1016/0021-9991(83)90136-5
[3] Harten A, Osher S. Uniformly high-order accurate nonoscillatory schemes. Ⅰ[J]. SIAM Journal on Numerical Analysis, 1987, 24(2): 279–309. DOI:10.1137/0724022
[4] Liu Xudong, Osher S, Tony Chan. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115(1): 200–212. DOI:10.1006/jcph.1994.1187
[5] Jiang Guangshan, Shu Chiwang. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126(1): 202–228. DOI:10.1006/jcph.1996.0130
[6] 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. DOI:10.1016/j.jcp.2005.01.023
[7] Feng Hui, Hu Fuxing, Wang Rong. A new mapped weighted essentially non-oscillatory scheme[J]. Journal of Scientific Computing, 2012, 51(2): 449–473. DOI:10.1007/s10915-011-9518-y
[8] 徐维铮, 孔祥韶, 吴卫国. 基于映射函数的三阶WENO改进格式及其应用[J]. 应用数学和力学, 2017, 38(10): 1120–1135.
[9] 钟巍, 贾雷明, 王澍霏, 等. 一类高效率高分辨率加映射的WENO格式及其在复杂流动问题数值模拟中的应用[J]. 力学学报, 2022, 54(11): 3010–3031. DOI:10.6052/0459-1879-22-247
[10] 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. DOI:10.1016/j.jcp.2007.11.038
[11] 刘旭亮, 武从海, 李虎, 等. WENO格式的一种新型光滑因子及其应用[J]. 空气动力学学报, 2023, 41(05): 20–34.
[12] Li Xiaogang, Li Guodong, Ge Yongbin. A new fifth-order finite difference WENO scheme for Dam-Break simulations[J]. Advances in Applied Mathematics and Mechanics, 2021, 13(01): 58–82. DOI:10.4208/aamm.OA-2019-0155
[13] Tang Shujiang. Improvements of the fifth-order WENO-JS-type scheme with normalized smoothing factor for gas dynamic Euler equations[J]. Applied Numerical Mathematics, 2023, 184: 301–324. DOI:10.1016/j.apnum.2022.10.010
[14] Anurag K, Bhavneet K, Kumar N T. A new type of improved third order WENO scheme with finite difference framework[J]. Chinese Journal of Physics, 2023, 83: 14–17. DOI:10.1016/j.cjph.2022.12.006
[15] Lith B S, Ten T B, Ijzerman W L. Embedded WENO: a design strategy to improve existing WENO schemes[J]. Journal of Computational Physics, 2017, 330(1): 529–549.
[16] 白晓雅, 郑秋亚, 梁益华. 求解欧拉方程的嵌入WENO格式[J]. 郑州大学学报(理学版), 2020, 52(03): 98–103.