由于双曲守恒律方程的解可能会在有限时间内产生间断, 因此数值求解此类型方程不仅要求数值格式具有高分辨率, 从而能够捕捉到更多的小尺度结构, 而且要求格式能够避免数值振荡.基本无振荡格式(ENO) 和加权基本无振荡格式(WENO), 既能够在间断处提供必要的耗散, 保持基本无振荡的特性, 又能在光滑区域实现高精度数值模拟, 因此成为求解含复杂流动问题的双曲守恒律方程的常见方法.ENO格式[1]将差商绝对值大小作为评估模板优劣的标准, 以此挑选出最为光滑的子模板对其进行插值重构.当最光滑子模板包含r个单元中心点时, ENO格式在光滑区域可达到r阶数值精度. 但在实际应用中, 过多的条件语句加大了计算成本, 而且舍弃了许多有用的模板信息.为了进一步提高数值精度, WENO格式对ENO格式进行改进, 将各个子模板的数值通量进行凸组合后的结果作为重构值.凸组合系数由相应的权函数决定. 若子模板被判定为光滑模板, 实际权重将趋近于理想权重, 格式将趋近于所对应的迎风格式, 从而在光滑区域实现高数值精度；而对于含间断的模板, 较大的局部光滑因子使得实际权重趋近于0, 从而实现基本无振荡.

Liu和Osher率先提出了首个WENO格式[2], 具有$ (r+1) $阶精度. 之后, Jiang和Shu[3]通过定义一个新的局部光滑因子, 发展了经典的WENO-JS格式, 新光滑因子将对应子模板上重构多项式的各阶导的$ L^2 $范数的平方和作为衡量子模板光滑程度的标准. 尽管WENO-JS格式将数值精度提升至$ (2r+1) $阶, 但是Henrick[4]等指出并证实了WENO-JS格式在极值点附近会产生精度下降的问题. 在Henrick等人的工作中, 系统地给出了5点WENO格式实现5阶收敛的充分条件和必要条件, 并通过设计一个映射函数克服了WENO-JS格式在极值点处精度下降的问题. 但映射算子的引入, 极大地增加了计算成本. Borges[5]等人独具匠心地引入了高阶全局光滑因子$ \tau_5 $来构造权函数, 并提出了WENO-Z格式. 该格式在计算成本上和WENO-JS格式相当, 但显著地提升了数值解在极值点处的精度, 减小格式的数值耗散. 然而, WENO-Z格式的数值模拟结果与指数$ q $之间存在着密切的关联性. 当$ q=1 $时, 5点WENO-Z格式在一阶极值点能够达到4阶精度. 当$ q=2 $时, 格式在一阶极值点能够达到5阶精度, 但$ q $值增大会导致格式耗散增大.

高阶全局光滑因子的引入为解决极值点精度降低的问题提供一个新的思路[5—8]. 不同于Jiang和Shu采用牛顿多项式插值的做法, Fan[6]等人根据拉格朗日多项式插值得出了一组新的局部光滑因子, 通过对插值多项式在$ x_i $处的高阶导数信息进行组合, 设计出了5阶、6阶、8阶的全局光滑因子[7]. 其中, 8阶全局光滑因子所对应的5点WENO-Z$ \eta $格式在极值点掉阶问题上有着重大的突破：在光滑区域甚至是一阶极值点处达到5阶收敛精度, 在二阶极值点处能够达到4阶收敛精度. 但在求解既含有光滑区域又含有非光滑区域的复杂问题时, 该格式会产生数值振荡. 在一些高维复杂问题中, 特别是双马赫问题, 其分辨率仍有很大地提升空间[8]. 许多研究者对WENO格式的权函数进行优化从而达到格式性能的提高[8—10]. 因此, 本文采用Fan等人提出的8阶全局光滑因子, 对5点WENO-Z$ \eta $格式的权函数进行优化: 将原先权函数中的数字1替换为一个关于$ \xi_k $的函数. 该函数在光滑区域为一个较大的常值, 从而使得格式更快地趋近于5阶迎风格式；而在间断区域, 函数值为一个小值, 从而强化了局部光滑因子对于权重的影响, 缓减了格式在含间断区域的数值振荡. 数值实验表明, 新格式(WENO-ZF)能够在不含极值点的光滑区域以及一阶极值点处达到5阶数值精度, 在二阶极值点处达到4阶数值精度. 在非光滑区域, 新格式能够更敏锐地捕捉间断, 同时保持基本无振荡的特性. 此外, 新提出的格式具有较好的延展性, 可以轻松地推广至更高阶的情况.

不失一般性, 以一维标量$ Euler $方程为例, 回顾WENO格式求解双曲守恒律方程的一般步骤：

$ \begin{equation} \frac{\partial u}{\partial t}+ \frac{\partial f(u)}{\partial x} = 0, \end{equation} $

(2.1)

其中$ u $为守恒变量, $ f(u) $为数值通量. 对$ [a, b] $进行网格离散有

$ \begin{equation} a={{x}_{\frac{1}{2}}}<{{x}_{\frac{3}{2}}}<...<{{x}_{N-\frac{1}{2}}}<{{x}_{N+\frac{1}{2}}}=b. \end{equation} $

(2.2)

网格单元长度为$ \Delta {x_i} = x_{i+\frac{1}{2}} - x_{i-\frac{1}{2}} $. 方便起见, 设网格等分. 此时网格中心点为$ x_i = a + i{\Delta x} - \frac{1}{2} {\Delta x}, i=1, 2, ..., N $. 对(2.1})式进行半离散, 得到一组常微分方程：

$ \begin{equation} \frac{d{u}_{i}}{dt}=-{\frac{\partial f}{\partial x}}|_{x=x_i}, \qquad i=1, 2, 3, ..., N. \end{equation} $

(2.3)

对数值通量进行高阶重构近似, 得到(\ref{eq3})式的近似形式：

$ \begin{equation} \frac{d{u}_{i}}{dt}\approx -\frac{\hat{f}_{i+\frac{1}{2}}-\hat{f}_{i-\frac{1}{2}}}{\Delta x}, \qquad i=1, 2, 3, ..., N. \end{equation} $

(2.4)

其中$ \hat {f}_{i\pm{\frac{1}{2}}} $表示单元边界处的通量的重构值, 为了计算稳定性, 对其进行通量分裂处理:

$ \begin{equation} \hat{f}_{i\pm\frac{1}{2}} = \hat{f}_{i\pm\frac{1}{2}}^+ + \hat{f}_{i\pm\frac{1}{2}}^-. \end{equation} $

(2.5)

本文采用的是$ Lax-Friedrichs $分裂, 有$ \hat {f}^{\pm}(u) = {\frac{1}{2}}({\hat{f}(u)} \pm au), a=max(|f'(u)|) $. 为了文章的简洁性, 本文只给出$ \hat{f}_{i+\frac{1}{2}}^+ $的重构方法, $ \hat{f}_{i+\frac{1}{2}}^- $可以通过对称性得到. 并且在接下来的讨论中, 将上标{“$ + $”}省略.

WENO格式采用低阶多项式重构值的非线性凸组合来构造边界值[3]. 对于5点WENO格式而言, 有

$ \begin{equation} \hat{f}_{i+\frac{1}{2}}=\sum\limits_{\text{m}=0}^{2}{{{\omega}_{m}}{\hat{f}_{i+\frac{1}{2}}^{m}}}, \end{equation} $

(2.6)

其中是$ \omega_m $非线性权重, $ \hat {f}^m _{i+\frac{1}{2}} $是基于第$ (m+1) $个子模板得到的低阶重构值:

$ \begin{equation} \begin{aligned} \hat{f}_{i+{\frac{1}{2}}}^{0}&={{\frac{1}{6}}(2{f}_{i-2}-7{f}_{i-1}+11{f}_{i})}, \\ \hat{f}_{i+{\frac{1}{2}}}^{1}&={{\frac{1}{6}}(-{f}_{i-1}+5{f}_{i}+2{f}_{i+1})}, \\ \hat{f}_{i+{\frac{1}{2}}}^{2}&={{\frac{1}{6}}(2{f}_{i}+5{f}_{i+1}-{f}_{i+2})}. \end{aligned} \end{equation} $

(2.7)

非线性权重$ \omega^{JS}_m $定义如下,

$ \begin{equation} \omega _m^{JS} = \frac{{\alpha _m^{JS}}}{{\sum\nolimits_{l = 0}^2 {\alpha _l^{JS}} }}, \alpha _m^{JS} = \frac{d_m}{({IS_m} + \varepsilon )^q}, \qquad { m = 0, 1, 2.} \end{equation} $

(2.8)

其中$ d_m $为理想权重, 有$ {d}_{0} = \frac{1}{10}, {d}_{1} = \frac{6}{10}, {d}_{2} = \frac{3}{10} $. 参数$ q $用来调节$ IS_m $对权重的影响程度, 在经典的WENO-JS格式中取$ q $为2. 参数$ \varepsilon $用来避免分母为0, 这里取值为$ 10^{-6} $. $ IS_m(m=0, 1, 2) $为局部光滑因子, 用来衡量各个子模板的光滑程度. 具体表达式如下：

$ \begin{equation} \begin{aligned} {IS}_{0} &= \frac{13}{12}{({f}_{i-2} - 2{f}_{i-1}+{f}_{i})^2}+{(\frac{{f}_{i-2}+4{f}_{i-1}+3{f}_{i}}{2})^2}, \\ {IS}_{1} &= \frac{13}{12}{({f}_{i-1} - 2{f}_{i}+{f}_{i+1})^2}+{(\frac{{f}_{i+1}-{f}_{i-1}}{2})^2}, \\ {IS}_{2} &= \frac{13}{12}{({f}_{i} - 2{f}_{i+1}+{f}_{i+2})^2}+{(\frac{-3{f}_{i}+4{f}_{i+1}-{f}_{i+2}}{2})^2}. \end{aligned} \end{equation} $

(2.9)

Henrick[4]等人系统地给出5点WENO格式最优阶收敛的充分条件, 有

$ \begin{equation} \omega_m = d_m + O(\Delta x^3). \end{equation} $

(2.10)

并指出$ \varepsilon $取值会影响格式的收敛精度. Henrick等人发现经典的5点WENO-JS格式, 在极值点处无法实现最优阶收敛, 且这一不足之处由于$ \varepsilon $取值较大而被掩盖. 当$ \varepsilon = 10^{-40} $时, 该格式在一阶极值点处只能达到3阶精度.

Borges等人率先引入了全局光滑因子$ \tau_5=|IS_2-IS_0| $来构造权函数, 这一创新为解决极值点处降阶问题开辟了新途径. WENO-Z格式的权函数如下：

$ \begin{equation} \omega _m^{Z} = \frac{{\alpha _m^{Z}}}{{\sum\nolimits_{l = 0}^2 {\alpha _l^{Z}} }}, \alpha _m^{Z} = {d_m}(1+(\frac{\tau_5}{IS_m + \varepsilon})^q ), \qquad { m = 0, 1, 2.} \end{equation} $

(2.11)

5点WENO-Z格式当$ q $取1时, 仍在极值点处存在掉阶问题：在一阶极值点处, 精度下降为4阶, 在二阶极值点处只有2阶数值精度. 虽然, $ q $取2能够使得格式在一阶极值点处达到5阶精度, 但是巨大的数值耗散会使得格式在间断附近产生更大的误差.

Fan等人沿用了Jiang和Shu等人的思想, 用所在子模板的重构多项式各阶导数的$ L^2 $范数的平方和来衡量子模板的光滑程度：

$ \begin{equation} IS_{m} = \sum\limits_{l=1}^{r-1} (\Delta x)^{2l-1} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \left( \hat{f}_m^{(l)}(x) \right)^2 \, dx. \end{equation} $

(2.12)

$ \hat{f}_m $是指在第$ (m+1) $个子模板上的插值多项式, $ \hat{f}^{(l)}_m $指的是$ \hat{f}_m $的第$ l $阶导数. 有所不同的是, Fan等人采用拉格朗日多项式进行插值重构[6], 由此得到的局部光滑因子$ IS^{\eta}_m $与传统的局部光滑因子$ IS_m $相比展现出更高的简洁性. 通过整合新构建的局部光滑因子中的信息, Fan等人设计了一组高阶全局光滑因子, 其中五阶全局光滑因子的表达式为

$ \begin{equation} \tau _{5}^{\eta }=|IS_{0}^{\eta }-IS_{2}^{\eta }|, \end{equation} $

(2.13)

六阶全局光滑因子的表达式为：

$ \begin{equation} \tau _{6}^{\eta }=|IS_{5}^{\eta }-\frac{IS_{0}^{\eta }+4IS_{1}^{\eta }+IS_{2}^{\eta }}{6}|, \end{equation} $

(2.14)

其中$ IS_{5}^{\eta } $为5点模板上的局部光滑因子, 由(2.12})式计算可得,

$ \begin{equation} IS_{5}^{\eta}=\frac{1}{144}[{{({{f}_{i-2}}-8{{f}_{i-1}}+8{{f}_{i+1}}-{{f}_{i+2}})}^{2}}+{{({{f}_{i-2}}-16{{f}_{i-1}}+30{{f}_{i}}-16{{f}_{i+1}}+{{f}_{i+2}})}^{2}}]. \end{equation} $

(2.15)

八阶全局光滑因子的表达式为：

$ \begin{equation} \begin{aligned} \tau _{81}^{\eta }=|(|f_{0}^{(1)}|-|f_{2}^{(1)}|)(f_{0}^{(2)}+f_{2}^{(2)}-2f_{1}^{(2)})|, \\ \tau _{82}^{\eta }={{(|f_{0}^{(1)}|-|f_{2}^{(1)}|)}^{2}}+{{(f_{0}^{(2)}+f_{2}^{(2)}-2f_{1}^{(2)})}^{2}}. \end{aligned} \end{equation} $

(2.16)

在数值模拟中, 采用Fan等人设计的8阶光滑因子所对应的5点WENO-Z$ \eta $格式, 能够在1阶极值点处达到5阶数值精度, 在2阶极值点处达到4阶数值精度, 这给极值点降阶问题带来重大突破. 但在求解同时包含光滑区域和非光滑区域的复杂问题时, 该格式会出现虚假数值振荡、分辨率低等问题. 鉴于现有研究多通过改进WENO格式的权函数来提升数值性能^{[5, 8, 9, 10]}, 本文着重对WENO-Z$ \eta $格式的权函数进行优化设计. 为了确保新格式在极值点处能够保持与8阶全局光滑因子所对应的WENO-Z$ \eta $格式相当的高数值精度, 同时有效缓解格式在间断区域附近产生的数值振荡, 提升其在复杂流域模拟中的分辨率, 我们对5点WENO-Z$ \eta $格式的加权策略进行了以下优化：

$ \begin{equation} \omega _m^{ZF} = \frac{{\alpha _m^{ZF}}}{{\sum\nolimits_{l = 0}^2 {\alpha _l^{ZF}} }}, \alpha _m^{ZF} = {d_m}(F+(\frac{\tau_{82}^\eta}{IS_m^\eta + \varepsilon})^q ), \qquad { m = 0, 1, 2.} \end{equation} $

(2.17)

其中$ q $用来控制非线性权重趋近于理想权重的速度, 通常取1. $ \varepsilon $取$ 10^{-40} $, 用来避免分母为0. $ IS_m^\eta $是局部光滑因子[7], 具体计算形式为

$ \begin{equation} IS_{m}^{\eta}={{(f_{m}^{(1)})}^{2}}+{{(f_{m}^{(2)})}^{2}}, m=0, 1, 2. \end{equation} $

(2.18)

其中$ f_{m}^{(1)} $可由第$ (m+1) $个子模板上的重构多项式在$ x_i $处的一阶导数值乘以$ \Delta x $得到, $ f_{m}^{(2)} $则由其在$ x_i $点处的二阶导数值乘以$ \Delta x^2 $得到, 具体表达式如下：

$ \begin{equation} \begin{cases} f_{0}^{(1)}&=\frac{1}{2}f_{i-2} -2f_{i-1} +\frac{3}{2}f_{i}, \\ f_{1}^{(1)}&=-{\frac{1}{2}}{f_{i-1}}+{\frac{1}{2}}{f_{i+1}}, \\ f_{2}^{(1)}&=-{\frac{3}{2}}{f_{i}}+2f_{i+1}-{\frac{1}{2}}{f_{i+2}}. \end{cases} \begin{cases} f_{0}^{(2)}&=f_{i-2} -2f_{i-1} + f_i, \\ f_{1}^{(2)}&=f_{i-1} -2f_{i} + f_{i+1}, \\ f_{2}^{(2)}&=f_{i}-2f_{i+1}+f_{i+2}.\\ \end{cases} \end{equation} $

(2.19)

计算得到局部光滑因子$ IS_{m}^{\eta} $的泰勒展开式：

$ \begin{equation} \begin{aligned} IS_{0}^{\eta }=&{{(f_{i}^{(1)}\Delta x)}^{2}}+\left( f{{_{i}^{(2)}}^{2}}-\frac{2}{3}f_{i}^{(1)}f_{i}^{(3)} \right)\Delta {{x}^{4}}+\left( -2f_{i}^{(2)}f_{i}^{(3)}+\frac{1}{2}f_{i}^{(1)}f_{i}^{(4)} \right)\Delta {{x}^{5}}\\&+\left( \frac{10}{9}f{{_{i}^{(3)}}^{2}}+\frac{7}{6}f_{i}^{(2)}f_{i}^{(4)}-\frac{7}{30}f_{i}^{(1)}f_{i}^{(5)} \right)\Delta {{x}^{6}}\\&+ \frac{1}{12}\left( -16f_{i}^{(3)}f_{i}^{(4)}-6f_{i}^{(2)}f_{i}^{(5)}+f_{i}^{(1)}f_{i}^{(6)} \right)\Delta {{x}^{7}}+O(\Delta {{x}^{8}});\\ \end{aligned} \end{equation} $

(2.20)

$ \begin{equation} \begin{aligned} IS_{1}^{\eta}=(f_{i}^{(1)}\Delta x)^2 +({f_{i}^{(2)}}^2+\frac{1}{3}f_{i}^{(1)}f_{i}^{(3)})\Delta x^4 +\frac{1}{180}(5{f_{i}^{(3)}}^2 + 30f_{i}^{(2)}f_{i}^{(4)}+ 3f_{i}^{(1)}f_{i}^{(5)})\Delta x^6 +O(\Delta {{x}^{8}}); \end{aligned} \end{equation} $

(2.21)

$ \begin{equation} \begin{aligned} IS_{2}^{\eta }=&{{(f_{i}^{(1)}\Delta x)}^{2}}+\left( f{{_{i}^{(2)}}^{2}}-\frac{2}{3}f_{i}^{(1)}f_{i}^{(3)} \right)\Delta {{x}^{4}}+\left(2f_{i}^{(2)}f_{i}^{(3)}-\frac{1}{2}f_{i}^{(1)}f_{i}^{(4)} \right)\Delta {{x}^{5}}\\&+\left( \frac{10}{9}f{{_{i}^{(3)}}^{2}}+\frac{7}{6}f_{i}^{(2)}f_{i}^{(4)}-\frac{7}{30}f_{i}^{(1)}f_{i}^{(5)} \right)\Delta {{x}^{6}}\\&+ \frac{1}{12}\left( 16f_{i}^{(3)}f_{i}^{(4)}+6f_{i}^{(2)}f_{i}^{(5)}-f_{i}^{(1)}f_{i}^{(6)} \right)\Delta {{x}^{7}}+O(\Delta {{x}^{8}}).\\ \end{aligned} \end{equation} $

(2.22)

因此有

$ \begin{equation} IS_m^\eta = \begin{cases} O(\Delta x^2) & f_i^{(1)} \neq 0, \\ O(\Delta x^4) & f_i^{(1)} = 0, f_i^{(2)} \neq 0, \\ O(\Delta x^6) & f_i^{(1)} = 0, f_i^{(2)} = 0. \end{cases} \end{equation} $

(2.23)

将局部光滑因子中的高阶信息进行组合, 得到8阶全局光滑因子[6]如下：

$ \begin{equation} \begin{aligned} \tau _{82}^{\eta}&={{(|f_{0}^{(1)}|-|f_{2}^{(1)}|)}^{2}}+{{(f_{0}^{(2)}+f_{2}^{(2)}-2f_{1}^{(2)})}^{2}}\\ & ={{(\frac{1}{2}f_{i}^{(4)}\Delta {{x}^{4}}+O(\Delta {{x}^{5}}))}^{2}}+{{(O(\Delta {{x}^{4}}))}^{2}}=O(\Delta {{x}^{8}}). \end{aligned} \end{equation} $

(2.24)

由上式易知, $ \tau_{82}^\eta $在光滑区域甚至是一阶或二阶极值点处, 都为$ O(\Delta x^8) $的小量. 对WENO-Z型权函数进行优化, 将常数1替换为函数$ F $:

$ \begin{equation} F=\lambda \frac{1+\xi _{\min }^{2}}{\sum\limits_{m=0}^{2}{{{d}_{m}}(1+\xi _{m}^{2})}}, {{\xi }_{m}}=\frac{\tau _{5}^{\eta }}{IS_{m}^{\eta }+\varepsilon }. \end{equation} $

(2.25)

其中, $ \lambda $是常数, $ \varepsilon ={{10}^{-40}} $. 接下来, 我们将从光滑区域与非光滑区域两个方面, 对$ F $的作用进行详细阐述. 当5点模板为光滑模板时, 将(2.20)式、(2.22})式带入(2.13)式有

$ \begin{equation} \tau _{5}^{\eta } = \begin{cases} O(\Delta x^5) & f_i^{(1)} \neq 0, \\ O(\Delta x^5) & f_i^{(1)} = 0, f_i^{(2)} \neq 0, \\ O(\Delta x^7) & f_i^{(1)} = 0, f_i^{(2)} = 0. \end{cases} \end{equation} $

(2.26)

将$ IS^{\eta }_{m} $的截断误差(2.23)式和$ \tau _{5}^{\eta } $的截断误差(2.26})式代入(2.25})式, 得到

$ \begin{equation} {{\xi }_{m}} = \begin{cases} O(\Delta x^3) & f_i^{(1)} \neq 0, \\ O(\Delta x) & f_i^{(1)} = 0, f_i^{(2)} \neq 0, \\ O(\Delta x) & f_i^{(1)} = 0, f_i^{(2)} = 0. \end{cases} F = \begin{cases} \lambda (1+O(\Delta {{x}^{6}})) & f_i^{(1)} \neq 0, \\ \lambda (1+O(\Delta {{x}^{2}})) & f_i^{(1)} = 0, f_i^{(2)} \neq 0, \\ \lambda (1+O(\Delta {{x}^{2}})) & f_i^{(1)} = 0, f_i^{(2)} = 0. \end{cases} \end{equation} $

(2.27)

因此在光滑区域, 函数值$ F $趋近于一个常数值$ \lambda $. 通过数值实验, 为了在保证数值格式充分展现ENO特性的基础上, 进一步最小化数值耗散, 最终选定了$ \lambda = 150 $作为最优参数. 又因为$ F $的值不会随着子模板而变化, 将(2.23)式、(2.24)式带入(2.17)式有

$ \begin{equation} \alpha_{m}^{ZF}=\begin{cases} {{d}_{m}}(F+O(\Delta {{x}^{6}})) & f_{i}^{(1)}\ne 0, \\ {{d}_{m}}(F+O(\Delta {{x}^{4}})) & f_{i}^{(1)}=0, f_{i}^{(2)}\ne 0, \\ {{d}_{m}}(F+O(\Delta {{x}^{2}}))& f_{i}^{(1)}=0, f_{i}^{(2)}=0.\\ \end{cases} \end{equation} $

(2.28)

$ \begin{equation} \omega_{m}^{ZF}=\begin{cases} {{d}_{m}}(1+O(\Delta {{x}^{6}})) & f_{i}^{(1)}\ne 0, \\ {{d}_{m}}(1+O(\Delta {{x}^{4}})) & f_{i}^{(1)}=0, f_{i}^{(2)}\ne 0, \\ {{d}_{m}}(1+O(\Delta {{x}^{2}}))& f_{i}^{(1)}=0, f_{i}^{(2)}=0. \end{cases} \end{equation} $

(2.29)

易知, WENO-ZF格式能够在光滑区域甚至是一阶极值点处满足5阶收敛的充分条件, 从而能够达到5阶收敛精度.

若整体五点模板上存在间断, 则$ \tau _{5}^{\eta }=O(1) $, $ {{(IS_{k}^{\eta })}_{\max }}=O(1) $, 因此$ {{\xi}_{\min }}=O(1) $. 对于未被间断穿过的子模版, (2.23)式仍成立, 因此有

$ \begin{equation} {{\xi }_{\max }} = \begin{cases} O(\Delta x^{-2}) & f_i^{(1)} \neq 0, \\ O(\Delta x^{-4}) & f_i^{(1)} = 0, f_i^{(2)} \neq 0, \\ O(\Delta x^{-6}) & f_i^{(1)} = 0, f_i^{(2)} = 0. \end{cases} \end{equation} $

(2.30)

故

$ \begin{equation} F = \begin{cases} O(1)/O(\Delta {{x}^{-4}})=O(\Delta x^4) & f_i^{(1)} \neq 0, \\ O(1)/O(\Delta {{x}^{-8}})=O(\Delta x^8) & f_i^{(1)} = 0, f_i^{(2)} \neq 0, \\ O(1)/O(\Delta {{x}^{-12}})=O(\Delta x^{12}) & f_i^{(1)} = 0, f_i^{(2)} = 0. \end{cases} \end{equation} $

(2.31)

因此$ F $在非光滑区域为一个很小的值. 这种做法能够增加局部光滑因子对实际权重的影响, 有效地抑制由高阶光滑因子所造成的数值振荡, 保持ENO特性. 此外, 由(2.17)式易知, 新格式具有很好的延展性, 能够轻松拓展到更高阶的情形. 新格式对WENO-Z$ \eta $格式中的权函数进行优化, 将常数1替换为函数$ F $, 因此称新格式为WENO-ZF格式.

在本节中, 我们采用WENO-JS、WENO-Z、WENO-Z$ \eta $($ \tau_{82} $)以及WENO-ZF等格式, 对一系列基准数值实验进行了模拟, 旨在验证WENO-ZF格式卓越的鲁棒性和高分辨率特性. 数值实验精心设计为三个部分：第一部分聚焦于线性对流方程的数值模拟, 旨在初步评估格式的精度与稳定性；第二部分和第三部分分别针对一维和二维欧拉方程展开数值模拟, 进一步验证格式在处理复杂流体动力学问题时的表现. 没有特殊说明, 在时间上采用3阶TVD龙格库塔方法迭代:

$ \begin{equation} \begin{aligned} & U^{(1)} = U^n + \Delta tL(U^n), \\ & U^{(2)} = \frac{1}{4}(3U^n + U^{(1)} + \Delta tL(U^{(1)})), \\ & U^{n+1} = \frac{1}{3}(U^n + 2U^{(2)} + 2\Delta tL(U^{(2)})). \end{aligned} \end{equation} $

(3.1)

首先, 测试格式WENO-ZF在光滑区域以及极值点处的数值精度. 对于线性对流方程

$ \begin{equation} \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0\quad -1<x<1 , \end{equation} $

(3.2)

给定初值条件

$ \begin{equation} u(x, t=0)= {u_0}(x), \end{equation} $

(3.3)

可计算出其精确解为

$ \begin{equation} u(x, t)=u(x-t, 0). \end{equation} $

给定一组初值条件$ u(x, 0)={{x}^{k}}{{e}^{x}} $, $ k=1, 2, 3. $在网格数为$ 20, 40, 80, 160, ..., 10240 $的均匀网格上计算$ \frac{\partial f}{\partial x} $在$ x=0 $处的误差. 当$ k=1 $时, 该函数在$ x=0 $处不存在极值点, 当$ k=2 $时, $ x=0 $为一阶极值点, 当$ k=3 $时, $ x=0 $为二阶极值点. 表 1、表 2、表 3分别展示了$ k=1, 2, 3 $时, 不同数值格式求解线性对流方程的误差和精度. 由数值结果可知, WENO-ZF与WENO-Z$ \eta $ ($ {{\tau }_{82}} $)格式均能在光滑区域甚至是一阶极值点处达到5阶数值精度, 在二阶极值点处格式能够达到4阶数值精度.

设初值条件为

$ \begin{equation} {u_0}(x)= {e^{\frac{-(x-90)^2}{400}}}(\cos(\frac{\pi}{8}(x-90))+\cos(\frac{\pi}{4}(x-90))), \quad 50<x<130. \end{equation} $

(3.5)

此算例的初值问题中含有多个极值点, 因此通常被用来测试数值格式在极值点处的数值精度. 对左右边界施加周期型边界条件, 图 3.1为各数值格式在网格数$ N $取$ 100 $, $ \Delta t=\frac{\Delta x}{2} $, 时间上采用4阶龙格库塔迭代到$ t=400 $时的数值结果, 可以看出WENO-ZF、WENO-Z$ \eta ({{\tau }_{82}}) $以及5阶迎风格式在极值点有着相似的数值模拟结果, 这三个格式与WENO-Z格式相比都能够在极值点处取得更接近于极值的数值解.

组合型激波包实验的初值条件为：

$ \begin{equation} u(x, 0)=\left\{ \begin{matrix} \frac{1}{6}(G(x, \beta , z-\delta )+G(x, \beta , z+\delta )+4G(x, \beta , z)), & 0.2\le x<0.4 \\ 1, & 0.6\le x<0.8 \\ \begin{matrix} 1-|10(x-0.1)|, & {} \\ \end{matrix} & 1\le x<1.2 \\ \frac{1}{6}(F(x, \alpha , a-\delta )+F(x, \alpha , a+\delta )+4F(x, \alpha , a)), & 1.4\le x<1.6 \\ \begin{matrix} 0 & {} \\ \end{matrix} & \\ \end{matrix} \right. \end{equation} $

(3.6)

其中$ G(x, \beta , z)={{e}^{-\beta {{(x-z)}^{2}}}} $, $ F(x, \alpha , a)=\sqrt{\max (1-{{\alpha }^{2}}{{(x-a)}^{2}}, 0)}. $$ z=-0.7 $, $ \delta=0.005 $, $ \beta=\frac{ln2}{36\delta^2} $, $ a=0.5 $, $ \alpha = 10. $

该算例初值条件涉及高斯波, 三角波, 方波以及半椭圆波, 通常被用来测试数值格式捕捉间断时的耗散与鲁棒性. Liu和Shen指出[8]高阶WENO-Z$ \eta $格式在复杂流域的数值模拟中, 会在间断附近产生数值振荡. 并且在长时间的数值模拟中, 数值振荡会加剧. 对左右边界施加周期型边界条件, 图 3.2展示了$ N=200 $, $ \Delta t=\frac{\Delta x}{2} $, $ t=10 $时, 求解线性对流方程的结果. 可以看到WENO-Z$ \eta ({{\tau }_{82}}) $在间断附近产生剧烈的数值振荡, 而WENO-ZF格式能够成功地抑制数值振荡并且能够更加敏锐地捕捉间断.

对于一维欧拉方程

$ \begin{equation} \frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}=0, \end{equation} $

(3.7)

其中

$ \begin{equation} U={{(\rho , \rho u, E)}^{T}}, F(U)={{(\rho u, \rho {{u}^{2}}+p, u(E+p))}^{T}}, \end{equation} $

(3.8)

对于理想气体有

$ \begin{equation} E=\frac{p}{\gamma -1}+\frac{1}{2}\rho {{u}^{2}}, \gamma =1.4. \end{equation} $

(3.9)

时间上仍然采用3阶龙格库塔方法进行迭代, 其中$ \Delta t=\frac{CFL\cdot \Delta x}{{{\max }_{i}}(|{{u}_{i}}|+{{c}_{i}})} $, $ CFL $这里取$ 0.5 $, $ c $表示声速, $ c=\sqrt{\frac{\gamma p}{\rho }} $.

$ Sod $问题的初值条件为[11]：

$ \begin{equation} \begin{aligned} (\rho, u, p) = \begin{cases} (0.125, 0, 0.1), \quad & -5\leq x <0, \\ (1, 0, 1), \quad & 0 \leq x\leq 5. \end{cases} \end{aligned} \end{equation} $

(3.10)

将WENO-Z格式在$ 2000 $个网格下的模拟结果作为$ \text{S}od $激波管问题的近似精确解. 左右边界设置为零梯度边界条件, 图 3.3为各数值格式在$ t=2 $, 网格数$ N=200 $的数值结果. 能够看出, WENO-ZF格式能够更加敏锐地捕捉数值间断且不产生数值振荡.

$ Lax $问题的初值条件为[12]：

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

(3.11)

将WENO-Z格式在$ 2000 $个网格下的数值结果作为$ Lax $激波管问题的近似精确解, 左右边界设置为零梯度边界条件, 图 3.4为各数值格式在$ t=1.3 $时刻, 网格数$ N=200 $的数值结果. 容易看出, WENO-ZF格式在间断处具有较低的数值耗散且不产生数值振荡.

123问题的初值条件为：

$ \begin{equation} \begin{aligned} (\rho, u, p) = \begin{cases} (1, -2, 0.4), \quad & -5\leq x <0, \\ (1, 2, 0.4), \quad & 0 \leq x\leq 5. \end{cases} \end{aligned} \end{equation} $

(3.12)

将WENO-Z格式在$ 2000 $个网格下的数值结果作为$ 123 $激波管问题的近似精确解, 左右边界设置为零梯度边界条件, 图 3.5为各数值格式在$ t=1 $, 网格数$ N=200 $的数值结果. 由于初值条件在$ x=0 $处存在强间断, 即使微小的数值振荡也可能导致密度出现负值, 从而引发整个数值计算的崩溃. WENO-Z$ \eta ({{\tau }_{82}}) $格式不能在上述条件下运行出结果, 正是基于这一原因. WENO-ZF相比于WENO-Z$ \eta ({{\tau }_{82}}) $格式, 能够在光滑区域具有更高的数值精度, 在间断区域维持ENO特性, 因此改善了WENO-Z$ \eta ({{\tau }_{82}}) $在$ 123 $问题中的不足, 不仅能够运行得出结果, 而且能够较为敏锐地捕捉到数值间断.

$ Shu-Osher $问题的初值条件为：

$ \begin{equation} \begin{aligned} (\rho, u, p) = \begin{cases} (3.857143, 2.629369, \frac{31}{3}), \quad & -5\leq x <-4, \\ (1+0.2\sin(5x), 0, 1), \quad & -4 \leq x\leq 5. \end{cases} \end{aligned} \end{equation} $

(3.13)

在此数值算例中, 一维的3马赫激波与一个扰动的密度场相互作用产生了许多小尺度结构与间断, 被用来测试数值格式在高频波模拟时的分辨率. 左右边界设置为零梯度边界条件, 将WENO-Z格式在$ 1000 $个网格下的数值结果作为该问题的近似精确解, 图 3.6展示了各数值格式在$ t=1.8 $时刻, 网格数$ N=200 $的数值结果. 能够看出WENO-ZF格式有较高的数值分辨率, 与其他格式相比, 具有更大的振幅.

爆炸波问题的初值条件为：

$ \begin{equation} \begin{aligned} (\rho, u, p) = \begin{cases} (1, 0, 1000), \quad & 0\leq x <0.1, \\ (1, 0, 0.01), \quad & 0.1 \leq x< 0.9, \\ (1, 0, 100), \quad & 0.9 \leq x \leq 1.0. \end{cases} \end{aligned} \end{equation} $

(3.14)

初始条件在$ x=0.1 $和$ x=0.9 $处含有强间断, 此案例用来测试格式的在间断处的耗散和基本无振荡特性. 将WENO-Z格式在$ 1000 $个网格下的数值结果作为该问题的近似精确解, 左右边界设置反射型边界条件, 图 3.7为各数值格式在$ t=0.038 $时刻, 网格数$ N=200 $的数值结果. 易见, WENO-ZF格式能够在极值点处取得更接近极值的值.

双马赫问题的初值条件为[13]：

$ \begin{equation} \begin{aligned} (\rho, u, v, p) = \begin{cases} (1.4, 0, 0, 1), \quad & y<1.732(x-0.1667), \\ (8, 7.145, -4.125, 116.50), \quad &otherwise. \end{cases} \end{aligned} \end{equation} $

(3.15)

初值条件被分成两个部分, 分割线从$ (0.6667, 0) $出发和$ x $轴形成$ {{60}^{\circ}} $的倾斜角. 下边界条件在点$ (0.6667, 0) $的左侧设置“后激波条件”, 在点$ (0.6667, 0) $的右侧设置为“对称型边界条件”. 上边界的分界点由以马赫10速度向右移动的激波决定, 左右边界采用通用的“流入流出边界条件”. 计算区域为$ [0, 4]\times [0, 1] $, 采用$ 1024\times 256 $个网格, 计算至$ t=0.2 $. 计算结果如图 3.8所示, WENO-Z$ \eta $格式在求解复杂流域的算例中, 并没有充分发挥高阶光滑因子的优势, 捕捉到较少的小尺度结构. 而WENO-ZF格式很好地弥补了这个不足, 能够捕捉到更精细的结构和更清晰的马赫杆.

R-T不稳定性算例的初值条件为[14]:

$ \begin{equation} \begin{aligned} (\rho, u, v, p) = \begin{cases} (2, 0, -0.025c\cos(8\pi x), 1+2y), \quad &if0 \leq y < \frac{1}{2}, \\ (1, 0, -0.025c\cos(8\pi x), y+\frac{3}{2}), \quad &if \frac{1}{2} \leq y < 1. \end{cases} \end{aligned} \end{equation} $

(3.16)

其中$ c $表示声速, $ c=\sqrt{\frac{\gamma p}{\rho }} $, $ \gamma $为比热比, $ \gamma =\frac{5}{3} $.

由初值条件可知, $ t=0 $时, $ y>0.5 $处流体密度较小, $ y<0.5 $处流体密度较大, 在$ y=0.5 $处形成一个间断. 上下边界设置为零梯度边界条件, 分别为：

$ (\rho , u, v, p)=(2, 0, 0, 1), (\rho , u, v, p)=(1, 0, 0, 2.5), $

左右采用对称型边界条件. 计算区域为$ [0, 0.25]\times [0, 1] $, 网格分辨率为$ 240\times 960 $, 计算至$ t=1.95 $. 图 3.9为各数值方法的数值模拟结果, 可以看出WENO-ZF格式在数值模拟中具有较高的分辨率.

本文所提出的WENO-ZF格式对WENO-Z$ \eta ({{\tau }_{82}}) $格式的权函数进行优化, 将原先权函数中的数字1替换为一个关于$ {\xi_{k}} $的函数$ F $. 函数$ F $的值与子模板$ k $无关, 且$ F $在光滑区域近似于一个较大的值$ \lambda $. 为了在数值格式的高鲁棒性和低数值耗散之间取得最佳平衡, 通过系统的数值实验, 最终确定最优参数$ \lambda $的取值为$ 150 $. 在光滑区域, WENO-ZF格式削弱了局部光滑因子对于实际权重的影响, 使得实际权重更快地趋近理想权重. 而在间断区域, 函数值$ F $为接近0的小值, 从而加强了局部光滑因子对于实际权重的影响, 基本消除了由高阶光滑因子所造成的数值振荡, 保持基本无振荡的特性. 数值实验表明, 5点WENO-ZF格式能够在光滑区域甚至是一阶极值点处达到5阶数值精度, 在二阶极值点处能够达到4阶数值精度. 在含间断的复杂流域问题求解时, WENO-ZF格式不仅能够保持基本无振荡的性质, 而且能够捕捉到更多的小尺度结构.