随着计算机的不断发展, 数值模拟是求解计算流体力学方程的重要手段之一. 双曲守恒律方程是流体力学的基本控制方程之一. 在特定参数下, 初始条件即便为光滑也可能产生激波和接触间断等不连续解, 为了消除间断处的数值振荡, 许多学者建立了一系列理论与方法. 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格式的收敛精度[6–9, 13], 通过光滑因子的修改来提高WENO格式的精度[10–12]等. Anurag Kumar等在经典三阶WENO模板中引入三点模板[14], 在此模板基础上进行重构, 并提出新的全局光滑因子, 进而得到一种改进的三阶WENO格式, 称为WENO-MZQ3. 本文在WENO-MZQ3格式的基础上, 将其推广到五阶WENO格式, 即在经典五阶WENO格式的三个三点模板的基础上引入五点模板, 通过对各模板相应的光滑因子的数值通量进行泰勒展开, 提出新的全局光滑因子, 对其进行加权处理得到新的五阶WENO格式. 通过理论分析可知新的五阶WENO格式在光滑区域相比五阶WENO-JS具有更高的精度. 通过数值算例验证, 新格式具有很好的分辨率和精度.
在本节中, 考虑如下一维标量守恒律方程
(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} $处的半离散形式如下
其中$ {\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]定义如下:
其中
WENO格式的数值性能由非线性权$ {\omega_k(k=0, 1, 2.)} $决定.下面给出经典五阶WENO-JS格式[5]的非线性权重:
其中$ d_0=0.1, d_1=0.6, d_2=0.3 $是线性权值, 设置$ {\varepsilon > 0} $为了避免分母为零, $ q $为可调正整数.$ {\beta_k(k=0, 1, 2.)} $是每个子模版的光滑因子, 表达式如下:
Embedded WENO-JS格式是在WENO-JS的基础上提出的.模板与WENO-JS模板一致, 其各个子模版对应的光滑因子为(2.6)、(2.7)、(2.8).Embedded WENO-JS格式贡献在于提出了新的非线性权重$ \alpha_k^{E} $, 通过修改非线性权重的值来提高格式的精度.对应的非线性权重$ \alpha_k^{E} $为[15, 16]:
其中$ c_0=\frac{2}{3}, c_2=\frac{6}{7}. $
单元界面$ {x=x_{i+\frac{1}{2}}} $处的数值通量为:
在本节中, 将介绍新型的五阶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}}} $对应的数值通量为:
其中$ {{p_0}, {p_1}, {p_2}} $对应的光滑因子由(2.6)、(2.7)、(2.8)式给出, 接下来介绍$ {p_5} $相对应的光滑因子[12]:
对光滑因子(2.6)、(2.7)、(2.8)、(2.13)式在$ {x={x_i}} $进行泰勒级数展开可得到[12, 13]:
由(2.14)、(2.15)、(2.16)式可得:
由(2.17)、(2.18)式可得:
由上式可知:
此时, 令$ {\tau _{HT}} = \left| {{\beta _5} - \beta } \right|. $定义以下非线性权:
在(2.20)式中, $ {\varepsilon} $是很小的正数, 只是为了避免分母为零.单元界面$ {x=x_{i+\frac{1}{2}}} $处的数值通量为:
在(2.21)式中: $ {\eta _5, \eta _0, \eta _1, \eta _2} $为线性权, 其和为1, 即$ {\eta _5 + \eta _0 + \eta _1 + \eta _2 = 1.} $其组合使格式在光滑区域达到更高精度.
文中时间推进采用三阶TVD Runge-Kutta方法.如下:
其中$ L(u_i)=- \frac{1}{{\Delta x}}\left( {{{\hat f}_{i + \frac{1}{2}}} - {{\hat f}_{i - \frac{1}{2}}}} \right). $在第三节中时间步长采用变时间步长, $ {\Delta t} $在一维和二维分别为:
其中$ {\sigma} $是CFL$ (Courant-Friedrichs-Lewy) $条件数, 本文算例中CFL条件数均取$ {0.5} $.
在本节中, 将五阶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给出了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.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.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.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)所示.
不同的$ {\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)式中, $ 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.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.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数值耗散更低.
以上二维算例的结果验证了WENO-HT格式在模拟复杂问题时与WENO-JS, Embedded WENO-JS格式相比具有更高的精度和分辨率.
本文在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格式具有更高的精度和分辨率.