众所周知, 双曲守恒律方程在众多领域中发挥着重要的作用. 而在过去的几十年里, 求解双曲守恒律方程的高阶方法的研究吸引了大量学者, 其主要挑战来自于处理激波和激波湍流的数值方法之间的矛盾特性, 这导致即使初始条件足够光滑, 解也可能发展为不连续解. 而在这些高阶数值方法中, WENO方法[1,2,3,4,5]引起了广泛的关注. 近年来, 王镇明、朱君等人[6]提出了一种五阶多分辨率WENO空间重构方法, 该方案的重构模板分别含由一个点, 三个点和五个点组成的中心模板, 而不是采用三个等大小的三点中心或偏离模板. 该方法构造简单, 易于得到高精度数值格式和推广到高维问题. 另外, WENO方案采用自适应模板的思想, 对候选模板分配非线性权并进行凸组合, 其中非线性权重取决于局部模版的光滑程度. 与通常的分量加权方法不同, 申义庆等人开发了一种共权WENO格式[7]. 共权是指在一个全局模板上, 一维空间中欧拉方程的分裂通量向量通常共享一组权重. 该方法计算效率较高, 具有良好的鲁棒性和较低的数值耗散性, 并有助于抑制相位误差. 本文基于多分辨率WENO有限差分格式和共权思想提出了一种共权多分辨率WENO有限差分格式, 用于求解一维和二维欧拉方程.
考虑以下一维欧拉方程:
其中守恒变量$ \mathbf {U} $, 无粘性通量$ \mathbf {F} $的表达式如下:
并且满足状态方程
其中$ \rho $, $ u $, $ E $和$ p $依次是密度, 速度, 总能量和压强, 而$ \gamma $是比热容.
通过欧拉方程的双曲性质, 我们可得
其中$ \Lambda = diag({\lambda _1}, {\lambda _2}, {\lambda _3}) = diag(u, u - c, u + c) $是雅克比矩阵A的特征值的对角矩阵, $ R $和$ L=R^{-1} $分别是该矩阵的右和左特征向量构成的矩阵, $ c = \sqrt {\frac{{\gamma p}}{\rho }} $代表声速.
通量$ \mathbf {F} $可以被分裂为正部和负部:
其中$ {\Lambda ^ + } = diag\left( {\lambda _1^ + , \lambda _2^ + , \lambda _3^ + } \right), {\Lambda ^ - } = diag\left( {\lambda _1^ - , \lambda _2^ - , \lambda _3^ - } \right) $分别是正部和负部的特征矩阵, $ \lambda _j^ + \ge 0, \lambda _j^ - \le 0, j = 1, 2, 3 $. 并且分裂通量$ \mathbf {F}^\pm $的广义公式可以写成如下形式:
分裂方程(2.5)的方法有很多, 例如Steger-Warming(SW)[8]通量分裂方法使用以下公式即可得到$ \mathbf {F}^+ $及$ \mathbf {F}^- $:
而本文选择使用全局Lax-Friedrichs(GLF)分裂方法[9]:
其中$ \alpha= {\max } \left\{ {{\lambda _j}} \right\} $是整个计算域上最大的特征值. 尽管采用GLF分裂的耗散更大些, 但得到的解的光滑性会更好.
注 对于二维情形, 通量分裂情况可参考[10], 这里不再赘述.
在本节中, 为了便于表达, 我们考虑将(2.1)缩写成以下一维标量形式下的双曲守恒方程:
为简单起见, 本文采用均匀网格, 网格单元$ {I_i} = \left[ {{x_{i - \frac{1}{2}}}, {x_{i + \frac{1}{2}}}} \right] $, 均匀单元长度表示为$ \Delta x = {x_{i + \frac{1}{2}}} - {x_{i - \frac{1}{2}}} $, 单元中心定义为$ {x_i} = \frac{{{x_{i + \frac{1}{2}}} + {x_{i - \frac{1}{2}}}}}{2} $. 在此均匀网格上, (3.1)的半离散格式为:
其中, $ u_i $是一个对点值$ u(x_i, t) $的数值近似, 而$ {{\widehat f}_{i + \frac{1}{2}}} $是一个数值通量. 自然地, 如果考虑了数值通量$ {{\widehat f}_{i + \frac{1}{2}}} $的五阶近似, 那么$ \frac{1}{h}\left( {{{\widehat f}_{i + \frac{1}{2}}} - {{\widehat f}_{i - \frac{1}{2}}}} \right) $是在$ x = x_i $处对$ f_x (u) $的五阶近似. 与经典的WENO5-JS方案相比, 王镇明、朱君等人提出了一种新的有限差分多分辨率WENO方法[6], 该方法的特点是利用一系列不等大小的分层中心空间模板设计高阶精度.
考虑数值通量$ f^+_{i+\frac{1}{2}} $在图 1中的三种模板($ {T_1} = \left\{ {{I_i}} \right\} $, $ {T_2} = \left\{ {{I_{i - 1}}, {I_i}, {I_{i + 1}}} \right\} $, $ {T_3} = \left\{ {{I_{i - 2}}, {I_{i - 1}}, {I_i}, {I_{i + 1}}, {I_{i + 2}}} \right\} $)下的重构多项式, 我们可得一个零次多项式$ q_1(x) $, 一个二次多项式$ q_2(x) $和一个四次多项式$ q_3(x) $满足对应点信息依次如下:
上述多项式的具体表达如下:
其中$ a_n(n=0, ..., 2) $和$ b_n(n=0, ..., 4) $分别是$ q_2(x) $和$ q_3(x) $的系数:
基于上述不同精度的多项式$ q_n $, 将它们重写成另一种等价的形式$ Q_n $, 以便以不同的精度顺序继续嵌套多分辨率WENO方案:
需要强调的是, 任意$ Q_n $中的线性权值之和必须为1(例如$ {\gamma _{1, 2}} + {\gamma _{2, 2}} = 1 $). 为了平衡不连续区域上的尖锐、本质无振荡的激波跃迁和平滑区域上的精度, 取$ {\gamma _{1, 2}} = \frac{1}{{11}}, {\gamma _{2, 2}} = \frac{{10}}{{11}};{\gamma _{1, 3}} = \frac{1}{{111}}, {\gamma _{2, 3}} = \frac{{10}}{{111}}, {\gamma _{3, 3}} = \frac{{100}}{{111}} $.
接下来需要计算$ Q_n $的光滑因子$ \beta_n, n = 1, 2, 3 $. 其中$ \beta_2 $和$ \beta_3 $按如下公式计算[5]:
其中$ k_n $是多项式$ Q_n(x) $的次数. 具体表达如下:
余下的光滑因子$ \beta_1 $的具体表达参考[11]:
其中
得到光滑因子后非线性权重的计算过程如下:
最后我们得到了$ {{\widehat f}_{i + \frac{1}{2}}} $最终的五阶重构:
值得一提的是, $ \hat f_{{i + \frac{1}{2}}}^ - $的重构相对于$ x_{i+1/2} $在空间上是镜像对称的.
申义庆等人[7]发现密度、压力和分裂后的能量方程的乘积可以作为变量来计算共同权值:(1)在激波和接触间断处, 密度发生跳跃;(2)分裂的能量通量包含速度的三次幂项(例如$ u^3 $)时能够具有迎风特性;(3)压力总是在激波处跳跃, 它可以帮助提高高速流动的稳定性, 其中动能比内能大得多.
与通常的通量分量加权方法不同, 共权是指在一个模板上, 一维空间中欧拉方程的分裂通量向量通常共享一组权重. 而对于二维欧拉方程而言, 需要将$ x $方向与$ y $方向分开考虑, 相同方向上的分裂通量向量共享一组权重. 共权WENO方案有两个显著的优点:其一是只需要计算一组权值, 具有较高的计算效率;其二是共权WENO方案对双曲方程组中数值通量的每个分量都保持相同的贡献.
因此, 本文基于上述思想结合多分辨率WENO格式提出了一个针对于欧拉方程的共权多分辨率WENO方案, 为简单起见, 以下算法仅给出(2.6)中通量分量的正部情形, 具体如下:
1. 由分裂通量(2.6)计算得到新的特殊计算变量$ \Gamma _i^ + $, 其中$ f_E $是与总能量相关的分量, 在一维情形下取$ f_E=f_3 $:
2. 计算系数$ a^E_n(n=0, ..., 2) $和$ b^E_n(n=0, ..., 4) $:
3. 计算新光滑因子.
3.1采用(4.2)中新系数并参照方程(3.12-3.13)计算$ {\beta^E_2} $和$ {\beta^E_3} $.
3.2首先由(3.16)和(3.17)式计算得到$ {\sigma _0}, {\sigma _1} $, 其中$ {\zeta _0} $, $ {\zeta _1} $的计算更改如下:
然后再计算$ {\beta^E_1} $:
4. 参照(3.18)计算新权重$ {\omega^E_n}, n=1, 2, 3 $.
5. 计算新重构$ \hat f_l^ + \left( x_{i + \frac{1}{2}} \right), l=1, 2, 3 $.
首先参照(3.6-3.10)计算$ Q_{l, n}, l=1, 2, 3 $:
然后计算$ \hat f_l^ + \left( x_{i + \frac{1}{2}} \right), l=1, 2, 3 $:
在本节中应用一些经典的无粘数值算例, 通过与原五阶有限差分多分辨率WENO格式的比较, 验证了所提出的针对欧拉方程的新的五阶有限差分共权多分辨率WENO格式的性能. 为了让结果更加有说服力, 本文对一维算例采用均方根误差(RMSE)和平均绝对误差(MAE)来衡量不同算法的效果:
其中$ y_i $为数值解, $ \mathop {{y_i}}\limits^ \wedge $为参考解. 而对于二维算例则是给出了另外两种五阶方法(WENO5-JS[5]和WENO5-Z[12])的收敛解来进行对比. 另外, 本文在时间上采用三阶龙格-库塔方法来求解欧拉方程:
其中, $ L $代表空间算子, $ U^n $是$ t^n $时刻的解. 对于一维问题而言, 时间步长采用:
二维算例则是采用:
CFL数使用$ \delta =0.5 $. $ \gamma $本文选为常数1.4.
例1 一维sod激波管问题, 该问题的初值条件如下:
这个问题被模拟到$ t=1.4 $, 分别有200个和300个网格点两种情形在图 2和图 3中被展示, 给出了参考解与多分辨率WENO格式和共权多分辨率WENO格式计算得到的密度分布. 其中参考解是基于3000个网格点由WENO5-JS格式计算的收敛解. 图 2和图 3之间的对比描述了随着网格点数的增加, 数值结果的收敛情况变好. 表 1给出了两种格式的均方根误差和平均绝对误差比较, 虽然所提出的新格式不能避免接触不连续点附近的数值振荡, 但它作为传统的分量加权方法可用于求解欧拉方程.
例2 一维Shu-Osher问题, 该问题的初值条件如下:
在这个数值算例中我们给出$ t=0.038 $时刻分别在网格点数为200和300时的两种情形, 分别如图 4和图 5所示, 其中图 4(b)和图 5(b)是复杂数值结果处的局部放大图. 通过图 4(b)和图 5(b)之间密度分布的对比, 我们能够更加直观地对比得出共权多分辨率WENO格式的数值解更贴近参考解, 具有更小的数值耗散. 表 2则是量化地指出了所提方法的优势. 值得一提的是, 本算例的参考解仍然是基于3000个网格点, 由WENO5-JS格式计算的收敛解.
例3 二维黎曼问题, 该问题的初值条件如下:
$ t=0.8 $时两种方案采用$ 400 \times 400 $网格点求解的密度等值线分别如图 6(c)和图 6(d)所示. 并且在图 6(a)和图 6(b)中, 本文还给出了WENO5-JS格式和WENO5-Z格式的收敛解. 相比图 6(a)和图 6(b)中的数值结果, 我们可以明显地看到, 共权多分辨率WENO格式在同一网格级上的表现优于一般的多分辨率WENO格式, 前者捕捉到了更微妙、更复杂的结构.
例4 作为一个检验高阶格式数值耗散的经典算例, 二维瑞利-泰勒不稳定问题描述了当加速度从重流体转向轻流体时, 不同密度的流体之间的界面不稳定性. 具体初值条件如下:
其中$ \alpha= \sqrt {\frac{{5 p}}{3 \rho }} $, 计算域为$ \left( {x, y} \right) \in \left[ {0, 0.25} \right] \times \left[ {0, 1} \right] $. 左右边界采用反射边界条件, 在顶部边界, 流量值设置为$ \left( {\rho , u, v, p} \right) = \left( {1, 0, 0, 2.5} \right) $, 而底部边界条件为$ \left( {\rho , u, v, p} \right) = \left( {2, 0, 0, 1} \right) $. 网格采用$ 240 \times 960 $, 图 7(c)和图 7(d)分别给出了两种多分辨率格式在$ t=1.95 $时的数值模拟结果. 考虑图 7(a)和图 7(b)中的收敛解, 可知由于耗散较少, 相比于一般的多分辨率WENO格式, 共权多分辨率WENO格式产生了更复杂的不稳定结构.
本文在多分辨率WENO格式的基础上引入共权思想对欧拉方程进行空间重构, 并给出了该格式求解欧拉方程的具体过程. 通过数值实验与原多分辨率WENO格式等进行对比, 得出了所提方法具有一定的鲁棒性和收敛精度.