双曲守恒律方程是计算流体力学中一类重要的方程, 如海洋学中的浅水波方程、空气动力学中的Euler方程和交通流领域中的LWR模型等. 即使初始条件充分光滑, 该类方程的解也可能会在某个时刻出现间断. 解在间断附近有着丰富和复杂的结构, 这给数值解法的构造带来了诸多挑战. 为了研究间断解, Lax^[1]提出了弱解的概念, 允许间断解的存在. 但弱解并不唯一, 如何在众多弱解中寻找具有物理意义的相关解呢? 随后, Lax^[2]依据热力学第二定律提出了熵稳定条件, 满足该条件的数值格式产生的数值解能收敛于唯一且具有物理意义的解. 为了实现熵稳定, Tadmor^[3]介绍了熵变量和熵势的概念, 设计了一类具有二阶精度的熵守恒格式, 并严格证明了一个三点格式只需含有比熵守恒格式更多的数值粘性则是熵稳定的. 基于此, 通过对熵守恒格式添加合适的耗散项就可以得到熵稳定格式. Lefloch等^[4]提出了由二阶熵守恒格式构造任意偶数阶熵守恒格式的方法, 并将此高阶熵守恒通量与满足保号性的ENO重构数值耗散算子相结合, 得到TECNO格式. 为了构造高阶熵稳定格式, 熵变量在每个单元交界面上重构值的跳跃与原始值的跳跃需保持符号一致.

目前, 仅有一小类重构满足保号性. 程晓晗^[5]通过对二次多项式进行恰当限制, 提出了一种三阶保号的重构. Fjordholm^[6]提出的ENO插值和基于Minmod限制器的二阶TVD重构均满足保号性. 然而, TVD重构在解的极值处可能会降低精度, ENO重构由于不稳定模板的选择也可能会导致精度降低. 众所周知, WENO重构优于ENO重构. WENO重构是对各个候选模板进行非线性加权, 每个模板权值的选取依赖于该模板上数值解的局部光滑性, 进而克服了ENO重构的缺点. 在光滑区域, WENO格式比ENO格式精度更高、收敛性更好; 而在间断附近, 仍保留ENO格式的良好性质. 基于此, Fjordholm和Ray^[7]提出了一种保号的三阶WENO重构方法(称为SP-WENO), 并用标量方程进行测试. 随后, Biswas和Dubey^[8]基于TECNO框架, 构造了保号的TVD重构和三阶WENO重构耗散算子. 更多结果参见文献[9-22].

本文以熵稳定格式为基础, 以TECNO框架作为设计原则, 构造了一种高阶保号熵稳定格式. 从理论推导和数值实验两方面验证了三阶紧致CWENO重构熵变量满足保号特性.

考虑如下守恒律方程

(2.1) $ \begin{equation} \begin{array}{ll} {\bf u}_t+{\bf f(u)}_x = 0, & {\forall}\ (x, t)\in {{\Bbb R}} \times {{\Bbb R}} _+, \\ {\bf u}(x, 0) = {\bf u}_0(x), & {\forall}\ x\in {{\Bbb R}} , \end{array} \end{equation} $

其中向量$ {\bf u}: {{\Bbb R}} \times {{\Bbb R}} _+ \mapsto {{\Bbb R}} ^m $为守恒变量, 向量$ {\bf f(u)}: {{\Bbb R}} ^m \to {{\Bbb R}} ^m $为相应的通量函数. 众所周知, 即使初值充分光滑, 方程$ (2.1) $的解可能会含有激波等间断解. 因此, 需要寻找物理相关的唯一弱解, 其满足如下熵不等式

(2.2) $ \begin{equation} \eta({\bf u})_t+q({\bf u})_x \le 0, \end{equation} $

这里$ (\eta, q) $是熵对, 熵函数$ \eta: {{\Bbb R}} ^m \to {{\Bbb R}} ^m $是一个凸函数, 熵通量函数$ q: {{\Bbb R}} ^m \to {{\Bbb R}} ^m $满足$ \partial_{\bf u} q = {\bf v}^T \partial_{\bf u} {\bf f} $, 其中$ {\bf v} = \partial_{\bf u}\eta $是熵变量.

考虑$ {{\Bbb R}} $中的笛卡尔网格$ \{x_i\}_{i\in {\Bbb Z}} $, 设网格长度为$ h: = \Delta x = x_{i+1}-x_i $, $ x_i $为空间离散化单元$ I_i = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}] $的中点值. 则方程$ (2.1) $的半离散守恒格式为

(2.3) $ \begin{equation} \frac{\rm d}{{\rm d}t}{\bf u}_i (t) = -\frac{1}{h_i}(\hat{{\bf f}}_{i+\frac{1}{2}}-\hat{{\bf f}}_{i-\frac{1}{2}}), \end{equation} $

这里$ {\bf u}_i (t) = {\bf u} (x_i, t) $, $ \hat{{\bf f}}_{i+\frac{1}{2}} $是与$ {\bf f(u)} $相容的数值通量, 满足$ \hat{{\bf f}}({\bf u} $, $ {\bf u} $, $ \cdots $, $ {\bf u}) = {\bf f(u)} $. 相应地, 格式$ (2.3) $称为熵守恒的, 如果它满足离散熵等式

(2.4) $ \begin{equation} \frac{\rm d}{{\rm d}t}\eta({\bf u}_i(t))+\frac{1}{h_i}(\hat{q}_{i+\frac{1}{2}}-\hat{q}_{i-\frac{1}{2}}) = 0, \end{equation} $

其中数值通量$ \hat{q}_{i+\frac{1}{2}} $与熵守恒通量$ q $相容, 即

(2.5) $ \begin{equation} \hat{q}_{i+\frac{1}{2}} = \hat{q}({\bf u}_{i-p+1}, \cdots, {\bf u}_{i+p}), \ \ \hat{q}({\bf u}, {\bf u}, \cdots, {\bf u}) = q({\bf u}). \end{equation} $

定理2.1^[3]   若$ \hat{{{\bf f}}}_{i+\frac{1}{2}} $对所有的$ i(i\in{\Bbb Z}) $满足

(2.6) $ \begin{equation} [{\bf v}]_{i+\frac{1}{2}}\cdot\hat{{\bf f}}_{i+\frac{1}{2}} = [\psi]_{i+\frac{1}{2}}, \end{equation} $

那么格式$ (2.3) $是熵守恒的. 其数值熵通量满足

(2.7) $ \begin{equation} q_{i+\frac{1}{2}} = \bar{{\bf v}}_{i+\frac{1}{2}}\cdot\hat{{\bf f}}_{i+\frac{1}{2}}-\bar{\psi}_{i+\frac{1}{2}}. \end{equation} $

这里$ \psi $为熵势, $ [\ast]_{i+\frac{1}{2}}: = (\ast)_{i+1}-(\ast)_{i} $, $ (\bar{\ast})_{i+\frac{1}{2}}: = \frac{(\ast)_{i+1}+(\ast)_{i}}{2} $.

注意到, 条件$ (2.6) $为$ m $个未知数提供了一个代数方程, 除标量情形$ (m = 1) $外, 其解不唯一. 即对于方程$ (2.1) $的标量形式$ (m = 1) $, 存在唯一的三点熵守恒通量^[6]

(2.8) $ \begin{equation} \tilde{f}_{i+\frac{1}{2}}(u_i, u_{i+1}) = \left\{\begin{array}{ll} { } \frac{\psi(v_{i+1})-\psi(v_{i})}{v_{i+1}-v_{i}}, {\quad} & u_i\neq u_{i+1}, \\ f(u_i), & u_i = u_{i+1}.\end{array}\right. \end{equation} $

但对于守恒律方程组$ (m>1) $, 熵守恒通量的表示不唯一.Tadmor等^[23]用一条沿着基向量的分段线性路径为其构造了熵守恒通量.

上述定义的熵守恒通量只有二阶精度, LeFloch等^[4]通过对二阶精度熵守恒通量进行线性组合, 构造了$ 2p(p\in{\Bbb N}) $阶精度熵守恒通量

(2.9) $ \begin{equation} \tilde{{\bf f}}_{i+\frac{1}{2}}^{2p} = \sum\limits_{r = 1}^{p}\alpha_r^p\sum\limits_{s = 0}^{r-1}\tilde{{\bf f}}({\bf u}_{i-s}, \cdots , {\bf u}_{i-s+r}), \end{equation} $

式中$ \alpha_1^p, \alpha_2^p, \cdots, \alpha_p^p $满足如下$ p $个线性方程

(2.10) $ \begin{equation} \sum\limits_{r = 1}^{p}r\alpha_r^p = 1, {\quad} \sum\limits_{i = 1}^{p}i^{2s-1}\alpha_r^p = 0\ (s = 2, 3, \cdots, p). \end{equation} $

特别地, 四阶熵守恒通量$ (p = 2) $为

(2.11) $ \begin{equation} \tilde{{\bf f}}_{i+\frac{1}{2}}^4 = \frac{4}{3}\tilde{{\bf f}}({\bf u}_i, {\bf u}_{i+1})- \frac{1}{6}\big(\tilde{{\bf f}}({\bf u}_{i-1}, {\bf u}_{i+1})+\tilde{{\bf f}}({\bf u}_i, {\bf u}_{i+2})\big). \end{equation} $

这些熵守恒格式适用于熵保持的光滑问题, 但在激波处会出现伪振荡. 而这些伪振荡产生的原因是: 熵守恒格式无法在激波附近产生足够的熵增量. 故需为熵守恒通量添加适当的耗散算子, 以抑制伪振荡的产生.

熵稳定格式包含两个关键要素:

(i) 如何恰当地使用熵变量;

(ii) 是否比熵守恒格式含有更多的数值粘性.

定理2.2^[23]  假设$ {\bf D}_{i+\frac{1}{2}}\geqslant0 $, $ \tilde{{\bf f}}_{i+\frac{1}{2}} $为熵守恒通量, 若数值通量的形式为

(2.12) $ \begin{equation} \hat{{\bf f}}_{i+\frac{1}{2}} = \tilde{{\bf f}}_{i+\frac{1}{2}}-\frac{1}{2}{\bf D}_{i+\frac{1}{2}}[{\bf v}]_{i+\frac{1}{2}}, \end{equation} $

那么, 格式$ (2.3) $是熵稳定的. 即它满足离散熵不等式

(2.13) $ \begin{equation} \frac{\rm d}{{\rm d}t}\eta({\bf u}_i(t))+\frac{1}{h_i}(\hat{q}_{i+\frac{1}{2}}-\hat{q}_{i-\frac{1}{2}})\leqslant0. \end{equation} $

设$ {\bf A} = {\bf R}{\bf \Lambda}{\bf R}^{-1} $为Jacobian矩阵, $ \widetilde{{\bf \Lambda}} $是由$ {\bf A} $的特征值组成的对角矩阵, $ {\bf R} $为相应的右特征向量矩阵. 则正耗散矩阵可以写为$ {\bf D}_{i+\frac{1}{2}} = ({\bf R}\widetilde{{\bf \Lambda}}{\bf B}{\bf R}^{\rm T})_{i+\frac{1}{2}} $, 其中$ \widetilde{{\bf \Lambda}} $为正对角矩阵, $ {\bf B} $为相应的缩放比例矩阵, 且满足$ {\bf R}^{-1}{\rm d}{\bf u} = {\bf BR}^{\rm T} {\rm d}{\bf v} $. 常用的耗散算子有$ \widetilde{{\bf \Lambda}}^{Roe} = |\widetilde{{\bf \Lambda}}| $, $ \widetilde{{\bf \Lambda}}^{Rusanov} = \max|{\bf \Lambda}|{\bf I} $.

数值熵的耗散, 特别是间断附近的耗散, 若熵过多, 间断附近会出现抹平现象; 若熵不足, 间断附近会产生伪振荡. 为构造恰当的数值耗散, Ismail和Roe^[9]对Euler方程的耗散算子做出如下修正

(2.14) $ \begin{equation} \widetilde{{\bf \Lambda}}_{EC1} = |{\bf \Lambda}|+\frac{1}{6}|[{\bf \Lambda}_{u\pm a}]|, \end{equation} $

(2.15) $ \begin{equation} \widetilde{{\bf \Lambda}}_{EC2} = |{\bf \Lambda}_{EC2}|+\alpha_{EC2}|[{\bf \Lambda}_{u\pm a}]|, \end{equation} $

其中

(2.16) $ \begin{equation} {\bf \Lambda}_{EC2} = diag((1+\beta)(u-a), u, (1+\beta)(u+a)), \end{equation} $

(2.17) $ \begin{equation} \alpha_{EC2} = (\alpha_{\max}-\alpha_{\min})(\max(0, sign(dM_{\max}-[M])))+\alpha_{\min}, \end{equation} $

这里$ [M] $代表通量交界面处Mach数的改变, $ [{\bf \Lambda}_{u\pm a}] = diag([u-a], 0, [u+a]) $. 根据经验, 选择$ \beta = \alpha_{\min} = \frac{1}{6} $, $ \alpha_{\max} = 2.0 $, $ dM_{\max} = 0.5 $.

同理, 可以将耗散算子推广到浅水波方程组中, $ \widetilde{{\bf \Lambda}} $的表达式如下

(2.18) $ \begin{equation} \widetilde{{\bf \Lambda}} = |{\bf \Lambda}|+\frac{1}{6}|[{\bf \Lambda}_{u\pm \sqrt{gh}}]|. \end{equation} $

对$ 2\lfloor\frac{k+1}{2}\rfloor $阶熵守恒通量$ \tilde{{\bf f}}_{i+\frac{1}{2}} $ (见式(2.9))添加一个$ k $阶耗散项, 就可以获得一个$ k $阶精度的熵稳定通量, 这里$ \lfloor n\rfloor $表示$ \leq n $的最大整数. 高阶熵稳定TECNO通量^[6]为

(2.19) $ \begin{equation} \hat{{\bf f}}_{i+\frac{1}{2}} = \tilde{{\bf f}}_{i+\frac{1}{2}}-\frac{1}{2}{\bf D}_{i+\frac{1}{2}} \langle{{\bf v}}\rangle_{i+\frac{1}{2}}. \end{equation} $

这里$ \langle{{\bf v}}\rangle $代表$ k $阶重构值在单元交界面$ x_{i+\frac{1}{2}} $处的跳跃, 即

(2.20) $ \begin{equation} \langle{{\bf v}}\rangle_{i+\frac{1}{2}} = {\bf v}_{i+\frac{1}{2}}^+-{\bf v}_{i+\frac{1}{2}}^-, \end{equation} $

$ {\bf v}_{i+\frac{1}{2}}^+ $和$ {\bf v}_{i+\frac{1}{2}}^- $是单元交界面$ x_{i+\frac{1}{2}} $处$ {\bf v} $的重构值.

根据定理2.2, 若标量熵变量$ {\bf w} = {\bf R}^T{\bf v} $的重构满足

(2.21) $ \begin{equation} sign\left({\bf R}_{i+\frac{1}{2}}^T\langle{\bf v}\rangle_{i+\frac{1}{2}}\right) = sign\left({\bf R}_{i+\frac{1}{2}}^T\left[{\bf v}\right]_{i+\frac{1}{2}}\right), \end{equation} $

则通量$ (2.19) $是熵稳定的.

定义2.1  在每一个单元交界面$ x_{i+\frac{1}{2}} $处, 若$ v $的重构值$ v_{i+\frac{1}{2}}^{\pm} $满足

(2.22) $ \begin{equation} sign\left(\langle{v}\rangle_{i+\frac{1}{2}}\right) = sign\left(\left[{v}\right]_{i+\frac{1}{2}}\right), \end{equation} $

则该重构在$ x_{i+\frac{1}{2}} $处满足保号性.

该文基于TECNO框架提出了一种鲁棒的、耗散较小的熵稳定格式, 该格式允许使用跨越跳跃间断的保号高阶重构, 且耗散项仅在保号处起作用. 即

(2.23) $ \begin{equation} \hat{{\bf f}}_{i+\frac{1}{2}} = \tilde{{\bf f}}_{i+\frac{1}{2}}-\frac{1}{2}{\bf D}_{i+\frac{1}{2}}\cdot {\rm d}v^{\ast}, \end{equation} $

其中$ {\rm d}v^{\ast} = S_{i+\frac{1}{2}}\cdot\langle{{\bf v}}\rangle_{i+\frac{1}{2}} $.

(i) 对于标量方程

(2.24) $ \begin{equation} S_{i+\frac{1}{2}} = \left\{\begin{array}{ll}1, & sign\left(\langle{v}\rangle_{i+\frac{1}{2}}\right) = sign\left(\left[{v}\right]_{i+\frac{1}{2}}\right)\neq 0, \\0, & else.\end{array}\right. \end{equation} $

(ii) 对于守恒律方程组：$ S_{i+\frac{1}{2}} $的第$ l $个分量为

(2.25) $ \begin{equation} s_{i+\frac{1}{2}}^l = \left\{\begin{array}{ll}1, & sign\left(\langle{w}^l\rangle_{i+\frac{1}{2}}\right) = sign\left(\left[{w}^l\right]_{i+\frac{1}{2}}\right)\neq 0, \\0, & else, \end{array}\right. \end{equation} $

其中$ w^l $代表$ {\bf w} $的第$ l $个分量.

三阶紧致CWENO重构是一个基于不同模板上插值多项式的凸组合. 该重构基于定义一个二次函数, 将该函数添加到线性插值中, 会在光滑区域达到三阶精度; 对于不连续或大梯度区域, 会自动切换到单侧线性重构. 对于标量方程, 该文用三阶紧致CWENO重构熵变量$ v $. 则单元交界面$ x_{i\pm1/2} $处的重构值为

(3.1) $ \begin{equation} \begin{array}{ll} v_{i+\frac{1}{2}}^- = w_{1, i}p_{1}(x_{i+\frac{1}{2}})+w_{2, i}p_{2}(x_{i+\frac{1}{2}})+w_{3, i}p_{3}(x_{i+\frac{1}{2}}), \\ v_{i-\frac{1}{2}}^+ = \widetilde{w}_{1, i}p_{1}(x_{i-\frac{1}{2}})+\widetilde{w}_{2, i}p_{2}(x_{i-\frac{1}{2}})+\widetilde{w}_{3, i}p_{3}(x_{i-\frac{1}{2}}), \end{array} \end{equation} $

其中$ p_{1}(x), p_{2}(x), p_{3}(x) $是三阶紧致CWENO重构函数^[10], 权函数$ w_{k, i}, \widetilde{w}_{k, i}(k, i\in{\Bbb Z}) $满足

(3.2) $ \begin{equation} w_{k, i}\ge0, \ \sum\limits_{k = 1}^3 w_{k, i} = 1, \; \; \widetilde{w}_{k, i}\ge0, \ \sum\limits_{k = 1}^3 \widetilde{w}_{k, i} = 1. \end{equation} $

引理3.1   若不等式

(3.3) $ \begin{equation} |v_{i+1}-v_i|>\max(|v_i-v_{i-1}|, |v_{i+2}-v_{i+1}|) \end{equation} $

成立, 则称离散集$ \{v_i\} $在区间$ [x_i, x_{i+1}] $上有一个局部显著跳跃. 也就是说, 在区间$ [x_i, x_{i+1}] $上的跳跃比其相邻左右区间的跳跃要大.

定理3.1  三阶紧致CWENO重构$ (3.1) $, 跨越跳跃间断满足保号性且局部有界. 即

(3.4) $ \begin{equation} 0\le \frac{\langle v\rangle_{i+\frac{1}{2}}}{[v]_{i+\frac{1}{2}}}\le2. \end{equation} $

证   设$ v(x) $中的不连续点在区间$ [x_i, x_{i+1}] $上. 因此, 在离散集$ \{v_i\} $中, 区间$ [x_i, x_{i+1}] $上的跳跃对应于式$ (3.3) $中的局部显著跳跃, 即

(3.5) $ \begin{equation} \mid[v]_{i-\frac{1}{2}}\mid<\mid[v]_{i+\frac{1}{2}}\mid, \ \ \mid[v]_{i+\frac{3}{2}}\mid<\mid[v]_{i+\frac{1}{2}}\mid. \end{equation} $

令$ r_i^+ = \frac{\left[v\right]_{i-\frac{1}{2}}}{\left[v\right]_{i+\frac{1}{2}}} $, $ r_i^- = \frac{\left[v\right]_{i+\frac{1}{2}}}{\left[v\right]_{i-\frac{1}{2}}} $, 则有

(3.6) $ \begin{equation} \left|r_i^+\right|<1, \ \ \left|r_{i+1}^-\right|<1. \end{equation} $

考虑三阶紧致CWENO重构:

(3.7) $ \begin{eqnarray} v_{i+\frac{1}{2}}^-& = &w_{1, i}(v_i+\frac{v_i-v_{i-1}}{2}) +w_{3, i}(v_i+\frac{v_{i+1}-v_i}{2}){}\\ &&+w_{2, i}\bigg\{v_i-\frac{1}{12}(v_{i+1}-2v_i+v_{i-1})+\frac{1}{4}(v_{i+1}-v_{i-1}) +\frac{1}{4}(v_{i+1}-2v_i+v_{i-1})\bigg\} , \end{eqnarray} $

(3.8) $ \begin{eqnarray} v_{i+\frac{1}{2}}^+& = &w_{1, i+1}(v_{i+1}-\frac{v_{i+1}-v_i}{2}) +w_{3, i+1}(v_{i+1}-\frac{v_{i+2}-v_{i+1}}{2}){}\\&&+w_{2, i+1}\bigg\{v_{i+1}-\frac{1}{12}(v_{i+2}-2v_{i+1}+v_i)-\frac{1}{4}(v_{i+2}-v_i) +\frac{1}{4}(v_{i+2}-2v_{i+1}+v_i)\bigg\} .{\quad} \end{eqnarray} $

那么熵变量在单元交界面$ x_{i\pm\frac{1}{2}} $处的跳跃重构值为

(3.9) $ \begin{eqnarray} \langle v \rangle_{i+\frac{1}{2}} & = &v_{i+\frac{1}{2}}^+ -v_{i+\frac{1}{2}}^-{}\\ & = &v_{i+1}-v_i-\frac{1}{2}\left(\widetilde{w}_{1, i}+\frac{1}{6}\widetilde{w}_{2, i} \frac{\left[v\right]_{i+\frac{3}{2}}}{\left[v\right]_{i+\frac{1}{2}}}-\frac{5}{6}\widetilde{w}_{2, i} +\widetilde{w}_{3, i}\frac{\left[v\right]_{i+\frac{3}{2}}}{\left[v\right]_{i+\frac{1}{2}}}\right) \left[v\right]_{i+\frac{1}{2}}{}\\ &&-\frac{1}{2}\left(w_{1, i}\frac{\left[v\right]_{i-\frac{1}{2}}}{\left[v\right]_{i+\frac{1}{2}}} +\frac{1}{6}w_{2, i}\frac{\left[v\right]_{i-\frac{1}{2}}}{\left[v\right]_{i+\frac{1}{2}}} +\frac{5}{6}w_{2, i}+w_{3, i}\right)\left[v\right]_{i+\frac{1}{2}}{}\\ & = &\Bigg\{1-\frac{1}{2}\left(\widetilde{w}_{1, i}+\frac{1}{6}\widetilde{w}_{2, i} \frac{\left[v\right]_{i+\frac{3}{2}}}{\left[v\right]_{i+\frac{1}{2}}}-\frac{5}{6}\widetilde{w}_{2, i} +\widetilde{w}_{3, i}\frac{\left[v\right]_{i+\frac{3}{2}}}{\left[v\right]_{i+\frac{1}{2}}}\right){}\\& &-\frac{1}{2}\left(w_{1, i}\frac{\left[v\right]_{i-\frac{1}{2}}}{\left[v\right]_{i+\frac{1}{2}}} +\frac{1}{6}w_{2, i}\frac{\left[v\right]_{i-\frac{1}{2}}}{\left[v\right]_{i+\frac{1}{2}}} +\frac{5}{6}w_{2, i}+w_{3, i}\right)\Bigg\}\left[v\right]_{i+\frac{1}{2}} . \end{eqnarray} $

在条件$ (3.6) $下, 由式$ (3.2) $知对所有的$ i\in{\Bbb Z}, 0\leq w_{k, i}\leq 1, 0\leq\widetilde{w}_{k, i}\leq 1 $, 可得

(3.10) $ \begin{eqnarray} 0& \leqslant &\Bigg\{1-\frac{1}{2}\left(\widetilde{w}_{1, i}+\frac{1}{6}\widetilde{w}_{2, i}r_{i+1}^–\frac{5}{6}\widetilde{w}_{2, i} +\widetilde{w}_{3, i}r_{i+1}^-\right){}\\ &&-\frac{1}{2}\left(w_{1, i}r_i^+ +\frac{1}{6}w_{2, i}r_i^+ +\frac{5}{6}w_{2, i}+w_{3, i}\right)\Bigg\}\leqslant 2 . \end{eqnarray} $

故由条件$ (3.9) $和$ (3.10) $, 证明了式$ (3.4) $.

本节采用熵稳定格式$ (2.12) $(用ES表示)和高阶熵稳定格式$ (2.23) $(用CCWENO3表示)分别对一维算例进行数值测试, 并与精确解(用Exact表示)比较. 对于时间层的推进, 采用三阶强稳定龙格库塔方法, 以确保在某些$ CFL $条件下, 标量问题的总熵$ E(t): = \sum_{i}\eta(u_{i}(t)) $随时间减小.

(4.1) $ \begin{equation} \begin{array}{ll} u^{(1)} = u^n+\Delta t L(u^n), \\ { } u^{(2)} = \frac{3}{4}u^n+\frac{1}{4}u^{(1)}+\frac{1}{4}\Delta t L(u^{(1)}), \\ { }u^{n+1} = \frac{1}{3}u^n+\frac{2}{3}u^{(2)}+\frac{2}{3}\Delta t L(u^{(2)}), \end{array} \end{equation} $

这里$ (L(u))_i = -\frac{1}{h}\left(\hat{{\bf f}}_{i+\frac{1}{2}}-\hat{{\bf f}}_{i-\frac{1}{2}} \right) $, $ \hat{{\bf f}}_{i\pm\frac{1}{2}} $表示数值通量函数.

(4.2) $ \begin{equation} u_t+u_x = 0. \end{equation} $

采用周期性边界条件, 在区间$ [-1, 1] $上计算如下初值问题.

$ (1) $光滑初始条件

(4.3) $ \begin{equation} u(x, 0) = \sin(\pi x), \ \ |x|\leq1. \end{equation} $

$ (2) $间断初始条件

(4.4) $ \begin{equation} u(x, 0) = \left\{\begin{array}{ll} 1, {\quad} &{ } |x|\leqslant \frac{1}{3}, \\ 0, &{ } \frac{1}{3}\leqslant |x|\leqslant 1.\end{array}\right. \end{equation} $

线性对流方程的精确解为$ u(x, t) = u(x-at, 0) $, 本文取熵函数$ \eta = \frac{u^2}{2} $, 耗散矩阵$ D_{i+\frac{1}{2}} = 1 $表示波传播速度的绝对值. 该问题的解描述了接触间断的移动, 可以很好地检测格式的高分辨率性. 本节对连续型初值问题$ (4.3) $分别使用$ L^1 $范数和$ L^\infty $范数计算收敛阶, 从表 1可以看出CCWENO3格式达到了三阶精度. 对于间断初值问题, 从图 1可以看出CCWENO3格式的分辨率更高.

(4.5) $ \begin{equation} u_t+\bigg(\frac{u^2}{2}\bigg)_x = 0. \end{equation} $

采用周期性边界条件, 在区间$ [-1, 1] $上计算如下初值问题.

$ (1) $光滑初始条件

(4.6) $ \begin{equation} u(x, 0) = u_0-u_1\sin(\pi x), \ \ |x|\leq1. \end{equation} $

$ (2) $间断初始条件

(4.7) $ \begin{equation} u(x, 0) = \left\{\begin{array}{ll} 1, &{ } |x|\leqslant \frac{1}{3}, \\ -1, {\quad} &{ } \frac{1}{3}\leqslant |x|\leqslant 1.\end{array}\right. \end{equation} $

(4.8) $ \begin{equation} u(x, 0) = 1+\frac{1}{2}\sin(\pi x), \ \ |x|\leq1. \end{equation} $

Burgers方程是一类非线性标量问题, 取熵函数$ \eta(u) = \frac{u^2}{2} $, 则熵通量函数为$ q(u) = \frac{u^3}{3} $, 可得二阶熵守恒通量

(4.9) $ \begin{equation} \tilde{f}(u_i, u_{i+1}) = \frac{u_i^2+u_iu_{i+1}+u_{i+1}^2}{6}. \end{equation} $

取EC1耗散矩阵

(4.10) $ \begin{equation} {\bf D}_{i+\frac{1}{2}} = |\bar{u}_{i+\frac{1}{2}}|+\frac{1}{6}|[u]_{i+\frac{1}{2}}|. \end{equation} $

对于连续性初值问题$ 4.6 $, 本文取$ u_0 = 0 $, $ u_1 = 0.5 $, 分别使用$ L^1 $范数和$ L^\infty $范数计算CCWENO3格式在$ T = 0.32 $时刻的收敛率. 从表 2可以看出, 本文所构造的CCWENO3格式达到了三阶精度. 对于间断初值为式$ (4.7) $的无粘Burgers方程, 其精确解包含左稀疏波和在点$ x\leq1/3 $处的激波. 图 2表明, 与ES格式相比, CCWENO3格式更能准确地捕捉到稀疏波与激波. 对于间断初值$ (4.8) $, 精确解保持光滑, 直到$ T = 1/\pi $时破裂, 然后形成一个以$ x = 0 $为中心的激波. 图 3表明, CCWENO3格式很好地抑制了激波附近的伪振荡.

理想气体的一维Euler方程为

(4.11) $ \begin{equation} {\bf u}_t+{\bf f}({\bf u})_x = 0, \end{equation} $

其中$ {\bf u} $是守恒型向量, $ {\bf f} $是通量,

(4.12) $ \begin{equation} {\bf u} = \left( \begin{array}{c} \rho\\ \rho u \\ E \end{array} \right), \ {\bf f}({\bf u}) = \left( \begin{array}{c} \rho u\\ \rho u^2+p \\ \rho uH \end{array} \right). \end{equation} $

其中$ \rho $是密度, $ u $是速度, $ E $是总能, $ p = (\gamma-1)(E-\frac{1}{2}\rho u^2) $是压力, $ H = (E+p)/\rho $是总焓, $ \gamma = 1.4 $是理想气体常数.

取熵对$ (\eta, q) $, 令$ s = \ln p-\gamma \ln \rho $为热力学熵, 有

(4.13) $ \begin{equation} \eta({\bf u}) = \frac{-\rho s}{\gamma-1}, \ q({\bf u}) = \frac{-\rho us}{\gamma-1}, \end{equation} $

可得熵变量$ {\bf v} = \frac{\partial \eta}{\partial {\bf u}} = \left(\frac{\gamma-s}{\gamma-1}-\frac{\rho u^2}{2p}, \frac{\rho u}{p}, -\frac{\rho}{p}\right)^{\rm T} $, 熵势$ \psi({\bf v}) = \rho {\bf u} $, 且熵守恒通量满足$ {\bf v}^{\rm T}\tilde{{\bf f}} = [\psi] $.

Ismail和Roe^[9]为Euler方程构造了$ {\bf v}^{\rm T}\tilde{{\bf f}} = [\psi] $的显式解. 定义参数向量$ {\bf z} $为

(4.14) $ \begin{equation} {\bf z} = \left( \begin{array}{c} z_1\\ z_2 \\ z_3 \end{array} \right) = \sqrt{\frac{\rho}{p}}\left( \begin{array}{c} 1\\ u \\ p \end{array} \right), \end{equation} $

则单元交界面$ x_{i+\frac{1}{2}} $处的熵守恒通量函数为

(4.15) $ \begin{equation} \tilde{{\bf f}}_{i+\frac{1}{2}} = \left( \begin{array}{c} \hat{\rho}\hat{u}\\ \hat{\rho}\hat{u}^2+\hat{p}_1 \\ \hat{\rho}\hat{u}\hat{H} \end{array} \right)_{i+\frac{1}{2}}, \end{equation} $

其中$ \hat{u} = \frac{\bar{z}_2}{\bar{z}_1}, \hat{\rho} = \bar{z}_1\bar{z}_3^{\ln}, \hat{p}_1 = \frac{\bar{z}_3}{\bar{z}_1}, \hat{p}_2 = \frac{\gamma+1}{2\gamma}\frac{z_3^{\ln}}{\bar{z}_1}+\frac{\gamma-1}{2\gamma}\frac{\bar{z}_3}{\bar{z}_1}, \hat{a} = \left(\frac{\gamma\hat{p}_2}{\hat{\rho}}\right)^{\frac{1}{2}}, \hat{H} = \frac{\hat{a}^2}{\gamma-1}+\frac{\hat{u}^2}{2} $, 对数平均定义为

(4.16) $ \begin{equation} \left(\ast\right)_{i+\frac{1}{2}}^{\ln} = \frac{[\ast]_{i+\frac{1}{2}}}{[\ln(\ast)]_{i+\frac{1}{2}}} . \end{equation} $

下面给出求解一维Euler方程$ (4.11) $的数值测试, 且均使用EC2耗散矩阵.

(1) Sod激波管问题，初始条件为

(4.17) $ \begin{equation} \left\{ \begin{array}{ll} \left(\rho_l, u_l, p_l\right) = \left(1, 0, 1\right), \ &x\leq0.5, \\ \left(\rho_r, u_r, p_r\right) = \left(0.125, 0, 0.1\right), \ &x>0.5. \end{array} \right. \end{equation} $

Sod激波管问题的初始间断包含左稀疏波, 接触间断和右激波. 本文计算了$ T = 0.16 $时刻$ [0, 1] $区间上的解. 即扰动到达计算区域边界之前的解. 从图 4可以清晰地看到: 与ES格式相比, CCWENO3格式能更好地捕捉激波, 稀疏波和接触间断, 分辨率更高, 鲁棒性更好, 且不产生伪振荡.

(2) Lax激波管问题，初始条件为

(4.18) $ \begin{equation} \left\{\begin{array}{ll}\left(\rho_l, u_l, p_l\right) = \left(0.445, 0.698, 3.528\right), \ &x\leq0.5, \\\left(\rho_r, u_r, p_r\right) = \left(0.5, 0, 0.571\right), \ &x>0.5.\end{array} \right. \end{equation} $

Lax激波管问题包含了真正的非线性场(稀疏波, 激波)和线性退化场(接触间断). 其计算区间均为$ [0, 1] $, 终止时刻为$ T = 0.16, $CFL条件为0.25. 图 5表明: 与ES格式相比, CCWENO3格式可以避免过度抹平现象, 提高了数值结果的分辨率, 具有更好的鲁棒性.

(3) Shu-Osher问题, 初始条件为

(4.19) $ \begin{equation} \left\{ \begin{array}{ll} \left(\rho_l, u_l, p_l\right) = \left(3.857143, 2.629369, 10.3333\right), \ &x\leq0.1, \\ \left(\rho_r, u_r, p_r\right) = \left((1+0.2\sin(50x)), 0, 1\right), \ &x>0.1. \end{array} \right. \end{equation} $

在该测试中, Mach 3激波与声波相互作用, 既有强的和弱的激波, 又有高度振荡但平滑的波. 该测试求解复杂, 可以很好地检测所构造格式的性能. 本文在具有自由流动边界的计算域$ [0, 1] $上求解Euler方程. 当取800个网格点时, 从图 6可以看出, CCWENO3格式较清晰地捕捉了解的结构, 结果分辨率高, 格式鲁棒性好.

(4.20) $ \begin{equation} {\bf u}_t+{\bf f}({\bf u})_x = 0, \end{equation} $

其中$ {\bf u} $是守恒型向量, $ {\bf f} $是通量,

(4.21) $ \begin{equation} {\bf u} = \left(\begin{array}{c}h \\ hu\end{array}\right), \; \; \{\bf f}({\bf u}) = \left(\begin{array}{c}hu \\{ } hu^2+\frac{1}{2}gh^2\end{array}\right), \end{equation} $

其中$ h $为水深, $ u $为平均流速, $ g $为重力加速度常数. 浅水波方程是在重力作用下对水体进行建模, 熵在这里代表总能量, 熵对$ (\eta, q) $为

(4.22) $ \begin{equation} \eta({\bf u}) = \frac{hu^2+gh^2}{2}, \q({\bf u}) = \frac{hu^3}{2}+guh^2. \end{equation} $

相应地, 熵变量和熵势的表达式分别为

(4.23) $ \begin{equation} {\bf v} = \left(\begin{array}{c}{ } gh-\frac{u^2}{2} \\ u\end{array}\right), \\psi({\bf u}) = \frac{1}{2}guh^2. \end{equation} $

熵守恒通量$ \tilde{{\bf f}} $的二阶形式^[19]:

(4.24) $ \begin{equation} \tilde{{\bf f}}_{i+\frac{1}{2}} = \left(\begin{array}{c}\bar{h}\bar{u}\\{ } \bar{h}\bar{u}^2+\frac{g}{2}\bar{h^2}\end{array}\right)_{i+\frac{1}{2}}, \end{equation} $

满足$ {\bf v}^T\tilde{{\bf f}} = [\psi] $.

耗散矩阵$ {\bf D}_{i+\frac{1}{2}} $由以下平均值构造

(4.25) $ \begin{equation} {\bf R}_{i+\frac{1}{2}} = \frac{1}{\sqrt{2g}}\left(\begin{array}{cc}1{\quad} & 1 \\\bar{u}-\sqrt{g\bar{h}}{\quad} & \bar{u}+\sqrt{g\bar{h}}\end{array}\right)_{i+\frac{1}{2}}, \end{equation} $

(4.26) $ \begin{equation} {\bf \tilde{\Lambda}}_{i+\frac{1}{2}} = \left(\begin{array}{cc}{ }|\bar{u}-\sqrt{g\bar{h}}|+\frac{1}{6}[u-\sqrt{gh}]{\quad} & 0 \\0{\quad} &{ } |\bar{u}+\sqrt{g\bar{h}}|+\frac{1}{6}[u+\sqrt{gh}]\end{array}\right)_{i+\frac{1}{2}}. \end{equation} $

(1) 溃坝问题, 初始条件为

(4.27) $ \begin{equation} h(x, 0) = \left\{\begin{array}{ll}2, \ \ &x<0 \\1, \ \ &x>0\end{array} \right. , \ \ \u(x, 0)\equiv0. \end{equation} $

溃坝问题的解由左行稀疏波和右行激波组成. 在区间$ [-1, 1] $上用ES和CCWENO3格式计算了$ T = 0.4 $时刻200个网格点的解. 从图 7可以看出, CCWENO3格式可以准确地捕捉到解的结构, 分辨率更高.

(2) 大型溃坝问题, 初始条件为

(4.28) $ \begin{equation} h(x, 0) = \left\{\begin{array}{ll}15, \ \ &x<0 \\1, \ \ &x>0\end{array}\right. , \ \ \u(x, 0)\equiv0. \end{equation} $

大型溃坝问题结果包含左行稀疏波和右行激波. 计算到$ T = 0.1 $时刻, 由图 8可知, CCWENO3格式比ES格式分辨率高, 验证了新格式的性能.

本文给出了高阶熵稳定格式的构造过程. 高阶熵稳定格式的构造需要使用一种满足保号性质的方法来重构单元交界面上的熵变量. 只有少数几种方法满足这一性质, 该文提出了一种三阶保号的紧致CWENO型重构, 并基于修正的TECNO框架, 将高阶熵守恒通量与高阶数值耗散算子结合得到了高阶精度保号的熵稳定格式, 该格式具有如下特点:

(1) 对于双曲守恒律方程组的光滑解可以达到任意阶精度.

(2) 满足离散熵不等式, 是严格熵稳定的.

(3) 每个单元交界面上重构的跳跃与相应单元格值的跳跃具有相同的符号.

(4) 保留紧致CWENO格式的优点, 重构跳跃值具有对称性且计算经济, 即在相同精度要求下所用节点少.

(5) 精准捕捉解的结构, 有效避免非物理现象, 间断附近基本无振荡, 高分辨率, 且具有良好的鲁棒性等.