数学物理学报, 2023, 43(2): 491-504

浅水波方程熵稳定格式的保平衡性

建芒芒,, 郑素佩,*, 封建湖,, 翟梦情,

长安大学理学院 西安 710064

Well-Balanced Preserving of Entropy Stable Schemes for Shallow Water Equations

Jian Mangmang,, Zheng Supei,*, Feng Jianhu,, Zhai Mengqing,

College of Science, Chang'an University, Xi'an 710064

通讯作者: *郑素佩, E-mail: zsp2008@chd.edu.cn

收稿日期: 2022-07-25   修回日期: 2022-10-17  

基金资助: 国家自然科学基金(11971075)
国家自然科学基金(11901057)

Received: 2022-07-25   Revised: 2022-10-17  

Fund supported: NSFC(11971075)
NSFC(11901057)

作者简介 About authors

建芒芒,E-mail:j3506687033@163.com

封建湖,E-mail:jhfeng@chd.edu.cn

翟梦情,E-mail:zhaimengqing2016@163.com

摘要

保平衡性是浅水波方程的一个重要性质, 满足保平衡性的格式理论上能够准确捕捉稳态的微小扰动. 针对带有底部地形源项的浅水波方程设计恰当的数值耗散算子并选择合适的源项离散方式以精确平衡非零通量与源项, 构造了一类保平衡的高阶熵稳定格式, 提出了高阶熵守恒格式以及熵稳定格式的保平衡性定理并给予理论证明. 利用新格式计算几个典型的数值算例, 数值结果表明新格式能够很好地处理稳态解的小扰动问题.

关键词: 保平衡性; 稳态的小扰动问题; 熵稳定格式; 浅水波方程

Abstract

Preserving well-balanced property is an important property for shallow water matrixs. The schemes with this property can capture small perturbations of steady state accurately in theory. For the shallow water matrixs with source terms, the first step is to construct an appropriate numerical dissipative operator and select a special discretization of source terms to accurately balance the non-zero flux and source terms, which is to obtain a class of high order well-balanced entropy stable schemes to keep balance. The new idea is to put forward the well-balanced preserving theorems for the high-order entropy conservative schemes and for the high-order entropy stable schemes. The detail proof process is also clearly given. Finally, several typical numerical examples show that new schemes can well deal with the small perturbation problems of steady-state solutions.

Keywords: Well-balanced; Steady-state problems with small perturbations; Entropy stable schemes; Shallow water matrixs

PDF (971KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

建芒芒, 郑素佩, 封建湖, 翟梦情. 浅水波方程熵稳定格式的保平衡性[J]. 数学物理学报, 2023, 43(2): 491-504

Jian Mangmang, Zheng Supei, Feng Jianhu, Zhai Mengqing. Well-Balanced Preserving of Entropy Stable Schemes for Shallow Water Equations[J]. Acta Mathematica Scientia, 2023, 43(2): 491-504

1 引言

浅水波方程是流体动力学中描述流体运动规律的数学模型, 在海洋、清洁能源、环境与水利工程等众多领域具有广泛应用, 该模型属于非线性双曲守恒律方程, 对应解的结构十分复杂, 很难得到解析解, 因此构造合适的数值格式成为一个研究热点. 目前已有多种数值方法可用来求解浅水波方程, 如TVD(Total Variation Diminishing)格式[1]、WENO(Weighted Essentially Non-oscillatory) 格式[2]、间断有限元(Discontinuous Galerkin)方法[3]以及旋转通量混合格式[4] 等. 一种好的数值格式应具有保平衡性, 即能够准确地保持静水稳态, 如何保持这种微妙的平衡是构造数值格式的一个巨大挑战. 在计算稳态解时通常会出现非物理波, 具有该性质的格式可在一定程度上抑制非物理波的出现.

稳态模拟是浅水波方程中十分重要的研究课题, 很多有趣的现象都会涉及计算静止湖泊的扰动, 不少学者致力于研究具有保平衡性的数值格式, Kurganov等[5]针对具有地形源项的圣维南系统提出了保稳态解的中心迎风格式, 该格式保证了水深的非负性, 可用来计算干旱地区的水流问题. 为了处理复杂计算域的问题, Liu等[6]对二维浅水波方程组构造了一个改进的保平衡中心迎风格式, 格式基于非结构化网格进行离散, 能够精确捕捉干湿界面. 中心迎风格式[7-8]基于黎曼扇进行积分, 不需要特征分解和黎曼求解器, 但是在计算过程中某些变量会出现负值, 与实际意义相违背, 需要对格式进行适当的调整. 除此之外, 早在1987年, Tadmor[9]设计了一种具有二阶精度的熵守恒格式, 并提出一个三点格式如果比熵守恒格式含有更多数值粘性则它是熵稳定的, 根据此性质得到了确定熵稳定格式的方法. 后来, Fjordholm等[10]给出了非平底地形浅水波方程的具有一阶和二阶耗散的保平衡熵稳定格式, 该格式可保持静水稳态和动水平衡态. 为了提高格式的精度, 2018年, 张海军等[11]在数值通量中嵌入限制器, 得到一类高分辨率熵稳定格式, 该格式结合了熵稳定格式和熵守恒格式的优点, 有效抑制了非物理现象的产生.

本文的主要目的是构造求解带源项浅水波方程的高阶保平衡熵稳定格式, 重点在于提出该格式的保平衡性定理并给出相关的理论证明. 首先基于文献[12]中提出的任意高阶熵稳定格式对原耗散项部分进行修正, 变量重构由ENO(Essentially Non-oscillatory) 重构替换为保号WENO-AO(WENO with Adaptive Order)重构[13], 从而得到一类新的高阶熵稳定格式. 为了证明格式的保平衡性, 先介绍了二阶熵守恒格式以及具有一阶耗散的低阶熵稳定格式的保平衡性, 再考虑对源项进行数值离散, 离散方式与所采用的偶数阶熵守恒通量的组合方式相对应, 进而证明高阶熵稳定格式的保平衡性定理, 并通过理论推导、数值实验验证新格式的保平衡性.

2 控制方程

浅水波方程可用来描述溃坝问题、水环境污染问题以及潮汐问题等, 下面给出浅水波方程的具体形式.

2.1 一维浅水波方程

考虑一维具有底部地形源项的浅水波方程

$\begin{matrix} {U}_{t}+ {F(U)}_{x}= {S}(x, {U}),\end{matrix}$

其中, 守恒变量$ {U}=(h \ hu)^{\rm T}$, 通量$ {F}=(hu \ hu^2+\frac{1}{2}gh^2)^{\rm T}$, 源项$ {S}=(0 \ - ghb_{x})^{\rm T}$, $h$ 代表水深, $u$为水流速度, 常数$g$为重力加速度, $b=b(x)$代表流体流过表面的底部地形, 方程(2.1)存在稳态解, 即$u \equiv 0, h+b \equiv $常数时, 非零通量与源项保持平衡.

2.2 二维浅水波方程

二维带源项浅水波方程的表达形式如下

$\begin{matrix} {U}_{t}+ {F(U)}_{x}+ {G(U)}_{y}= {S} (x,y, {U}),\end{matrix}$

其中守恒变量${U}=(h \ hu \ hv)^{\rm T}$, $x$方向通量为$ {F}=(hu \ hu^2+\frac{1}{2}gh^2 \ huv)^{\rm T}$, $y$方向通量为$ {G}=(hv \ huv \ hv^2+\frac{1}{2}gh^2 )^{\rm T}$, 源项$ {S}=(0 \ - ghb_{x} \ - ghb_{y})^{\rm T}$, $h$代表水深, $u,v$分别为$x,y$方向的水流速度, $g$为重力加速度常数, $b=b(x, y)$代表底部地势, 方程(2.2)同样存在稳态解, 即$u \equiv 0,$$ v \equiv 0, $$h+b \equiv $常数时, 各个方向的非零通量与源项保持平衡.

3 离散格式

下面以一维浅水波方程数值格式的构造过程为例来阐述. 浅水波方程是非线性双曲守恒律方程, 因此, 即使初始条件充分光滑, 浅水波方程的解也可能会在有限的时间内产生间断, 从而造成经典解理论不再适用, 为此, Lax[14]引入了弱解的定义, 但由于弱解并不唯一, 又提出了满足熵条件的解是问题有物理意义的解, 非平底地形浅水波方程的熵稳定条件可参见文献[10].

考虑网格均匀划分, 方程(2.1)的半离散有限体积格式为

$\begin{matrix} \frac{{\rm d}}{{\rm d}t} {U}_{i}= -\frac{1}{\triangle x}( {F}_{i+1/2}- {F}_{i-1/2})- {S}_{i}\,, \end{matrix}$

${U}_{i}=-\frac{1}{\triangle x}\int_{I_{i}} {U}(x,t){\rm d}\textit{x}$$ {U}(x,t)$在网格单元$I_{i}=[x_{i-1/2},x_{i+1/2}]$上的单元平均, $ {F}_{i \pm 1/2}$是与通量${F}$相容的数值通量.

3.1 熵守恒格式

定义3.1[10]$ {F}$相容的通量$ {F}_{i+1/2}$是熵守恒数值通量, 如果等式

$\begin{matrix} {[V]}^{\rm T}_{i+1/2}\cdot {F}_{i+1/2}=[\psi]_{i+1/2}+g[b]_{i+1/2}\overline{h}_{i+1/2}\overline{u}_{i+1/2}\end{matrix}$

成立, 那么满足该条件的形如(3.1)式的格式称为熵守恒格式, 且具有二阶精度. 类似平底地形浅水波方程的通量构造方式, 得出具有非平底地形源项的浅水波方程的数值通量和源项离散

$\begin{matrix} \begin{array}{lll}{F}^{EC}_{i+1/2}=\left[{\begin{array}{c}\overline{h}_{i+1/2}\overline{u}_{i+1/2}\\[2mm] \frac{g}{2}\overline{h}^2_{i+1/2}+ \overline{h}_{i+1/2}(\overline{u}_{i+1/2})^2\end{array}}\right], \\[7mm] {S}^{EC}_{i}=\left[{\begin{array}{c}0\\[2mm] \frac{g}{2\triangle x}(\overline{h}_{i+1/2}[b]_{i+1/2}+ \overline{h}_{i-1/2}[b]_{i-1/2})\end{array}}\right].\end{array}\end{matrix}$

3.2 熵稳定格式

熵守恒格式在光滑区域计算准确, 但在间断处会有熵耗散, 产生非物理现象, 为了解决该问题, 在熵守恒通量的基础上添加适当的耗散项以获得熵稳定数值通量.

定义3.2[10] 若数值耗散系数矩阵$ {D}_{i+1/2}\ge 0$, $ {F}^{EC}_{i+1/2}$是由(3.2) 式得到的二阶熵守恒通量, 如果数值通量的形式为

$\begin{matrix} {F}^{ES1}_{i+1/2}= {F}^{EC}_{i+1/2}-\frac{1}{2} {D}_{i+1/2} {[V]}_{i+1/2} \,,\end{matrix}$

那么称格式(3.1)为熵稳定格式.熵变量$ {V}=[g(h+b)-\frac{1}{2}u^2,u]^{\rm T},$$ {D}_{i+1/2}= {R}_{i+1/2}| {{\Lambda}}_{i+1/2}|$$\cdot {R}^{\rm T}_{i+1/2},$${R}_{i+1/2}$$ {\Lambda}_{i+1/2}$分别代表右特征向量矩阵和由特征值所构成的对角矩阵的离散形式, 具体的表达方式参见文献[11].

3.3 高阶熵稳定格式

观察熵稳定格式的定义方式,其中的${R},{{\Lambda}}, {V}$都是由熵守恒变量所表示, 耗散部分涉及熵变量差分, 格式整体仅有一阶精度. 为了提高格式精度, 从两方面出发, 一方面是将二阶熵守恒通量替换为高阶即任意偶数阶熵守恒通量;另一方面考虑在单元交界面处对熵变量进行重构, 常见的变量重构有ENO重构[15], WENO重构[16], CWENO(Central WENO) 重构[17]等, 本文采用WENO-AO重构, 该重构结合了三阶紧致CWENO重构和自适应WENO重构的思想, 具有保号性[13].

首先考虑熵守恒通量, 由(3.2)式所计算出的熵守恒通量是二阶精度的, LeFloch等[12]通过对二阶熵守恒通量进行线性组合得到任意偶数阶$(2p)$的熵通量, 形式如下

$\begin{matrix} \widetilde{{F}}^{2p}_{i+1/2}=\sum\limits_{r=1}^p\alpha^p_{r}\sum\limits_{s=0}^{r-1}\widetilde{{F}}({U}_{i-s},{U}_{i-s+r}) \,, \end{matrix} $

其中系数$\alpha^p_{1},\alpha^p_{2},\cdots,\alpha^p_{p}$满足$\sum\limits_{r=1}^p r\alpha^p_{r}=1,\ \ \sum\limits_{r=1}^p i^{2s-1}\alpha^p_{r}=0\ \ (s=2,3,\cdots,p). $

熵变量采用WENO-AO重构, 具体的构造过程详见文献[13], 记熵变量${V}$在单元交界面$x_{i+1/2}$处的左右重构值为

$\begin{aligned}\boldsymbol{V}_{i+1 / 2}^{-}= & \omega_4\left(\frac{\boldsymbol{V}_{i-2}}{15}-\frac{4 \boldsymbol{V}_{i-1}}{15}+\frac{11 \boldsymbol{V}_i}{15}+\frac{\not 17 \boldsymbol{V}_{i+1}}{30}-\frac{\boldsymbol{V}_{i+2}}{10}\right) \\& +\omega_1\left(\frac{3 \boldsymbol{V}_i}{2}-\frac{\boldsymbol{V}_{i-1}}{2}\right)+\omega_2\left(\frac{5 \boldsymbol{V}_{i+1}}{12}+\frac{2 \boldsymbol{V}_i}{3}-\frac{\boldsymbol{V}_{i-1}}{12}\right)+\omega_3\left(\frac{\boldsymbol{V}_{i+1}}{2}+\frac{\boldsymbol{V}_i}{2}\right),\end{aligned}$
$\begin{aligned}\boldsymbol{V}_{i+1 / 2}^{+}= & \widetilde{\omega}_4\left(\frac{17 \boldsymbol{V}_i}{30}-\frac{\boldsymbol{V}_{i-1}}{10}+\frac{11 \boldsymbol{V}_{i+1}}{15}-\frac{4 \boldsymbol{V}_{i+2}}{15}+\frac{\boldsymbol{V}_{i+3}}{15}\right) \\& +\widetilde{\omega}_1\left(\frac{\boldsymbol{V}_i}{2}+\frac{\boldsymbol{V}_{i+1}}{2}\right)+\widetilde{\omega}_2\left(\frac{5 \boldsymbol{V}_i}{12}+\frac{2 \boldsymbol{V}_{i+1}}{3}-\frac{\boldsymbol{V}_{i+2}}{12}\right)+\widetilde{\omega}_3\left(\frac{3 \boldsymbol{V}_{i+1}}{2}-\frac{\boldsymbol{V}_{i+2}}{2}\right).\end{aligned}$

定义高阶熵稳定数值通量

$\begin{matrix} {F}^{HES}_{i+1/2}=\widetilde{ {F}}^{2p}_{i+1/2}-\frac{1}{2} {D}_{i+1/2} {\langle V \rangle}_{i+1/2}\,,\end{matrix}$

注意, $ {\langle V \rangle}_{i+1/2}$代表重构值在单元交界面处的跳跃, 即${\langle V \rangle}_{i+1/2}={V}^{+}_{i+1/2}-{V}^{-}_{i+1/2}$.

进一步, 高阶熵稳定格式可表示为

$\begin{matrix}\frac{{\rm d}}{{\rm d}t} {U}_{i}= -\frac{1}{\triangle x}( {F}^{HES}_{i+1/2}- {F}^{HES}_{i-1/2})- \widetilde{{S}}^{2p}_{i}\,,\end{matrix}$

其中$\widetilde{{S}}^{2p}_{i}=\sum\limits_{r=1}^p\alpha^p_{r}\sum\limits_{s=0}^{r-1}\widetilde{{S}}({U}_{i-s},{U}_{i-s+r})$, 此时, 源项的离散方式与熵守恒通量的组合方式相对应.常用的四阶$(p=2)$和六阶$(p=3)$熵守恒通量[12]

$\begin{matrix}\widetilde{{F}}^{4}_{i+1/2}&=&\frac{4}{3}\widetilde{{F}}({U}_{i},{U}_{i+1})-\frac{1}{6}\left(\widetilde{{F}}({U}_{i-1},{U}_{i+1})+\widetilde{{F}}({U}_{i},{U}_{i+2})\right),\\\qquad \widetilde{{F}}^{6}_{i+1/2}&=&\frac{3}{2}\widetilde{{F}}({U}_{i},{U}_{i+1})-\frac{3}{10}\left(\widetilde{{F}}({U}_{i-1},{U}_{i+1})+\widetilde{{F}}({U}_{i},{U}_{i+2})\right)\\&&+\frac{1}{30}\left(\widetilde{{F}}({U}_{i-2},{U}_{i+1})+\widetilde{{F}}({U}_{i-1},{U}_{i+2})+\widetilde{{F}}({U}_{i},{U}_{i+3})\right).\end{matrix}$

4 格式的保平衡性理论

鉴于熵稳定格式本身在求解双曲守恒律问题时已经具备很好的性质, 格式的构造方式是本文的一小部分, 重点在于格式的保平衡性理论. 具有非平坦地形的浅水波方程存在稳态解, 而格式的保平衡性是计算稳态的必要条件, 其关键在于确定合适的离散方式使得非零通量与源项相互抵消, 精确保持该平衡是构造数值格式的一大难点. 当在粗网格中计算方程(2.1)的稳态解时, 稳态的小扰动也会被格式所放大, 所以格式必须满足保平衡性, 才能准确刻画稳态的扰动问题.

稳态解指方程(2.1)与时间变量无关的状态解, 本文主要研究静水稳态, 对一维问题即

$\begin{matrix}u \equiv 0, h+b \equiv \mbox{常数}.\end{matrix}$

如何维持这样的稳态解以及恰当处理带有小扰动的情况是构造数值格式的重点, 格式如果能够保持稳态的离散形式, 则可称为是保平衡的, 有时也表达成“和谐的”。湖泊上的波浪或者深海中的海啸波都是典型的静水扰动问题.

4.1 熵守恒格式的保平衡性

具有(3.3)式熵守恒通量和源项的熵守恒格式(3.1)可整理为

$\begin{matrix}\frac{{\rm d}}{{\rm d}t}h_{i}& =&-\frac{1}{\triangle x}(\overline{h}_{i+1/2}\overline{u}_{i+1/2}-\overline{h}_{i-1/2}\overline{u}_{i-1/2}),\\ \frac{{\rm d}}{{\rm d}t}(h_{i}u_{i})&=&-\frac{1}{\triangle x}(\overline{h}_{i+1/2}(\overline{u}_{i+1/2})^2+\frac{g}{2}\overline{h}^2_{i+1/2}-\overline{h}_{i-1/2}(\overline{u}_{i-1/2})^2-\frac{g}{2}\overline{h}^2_{i-1/2})\\&&-\frac{g}{2\triangle x}(\overline{h}_{i+1/2}[b]_{i+1/2}+(\overline{h}_{i-1/2}[b]_{i-1/2}).\end{matrix}$

引理4.1[10] 熵守恒格式(4.2)具有保平衡性, 即给定初始数据: $u_{i}=0, h_{i}+b_{i}= \mbox{常数}$, 对任意的整数$i$都成立, 那么由格式(4.2)计算出的解满足

$\begin{matrix}\frac{{\rm d}}{{\rm d}t}h_{i} \equiv 0, \ \ \ \frac{{\rm d}}{{\rm d}t}(h_{i}u_{i}) \equiv 0.\end{matrix}$

当然, 也可以通过上式说明格式能够保持静水稳态.

4.2 熵稳定格式的保平衡性

熵守恒格式可以保持静水稳态, 能够用来计算稳态的小扰动问题, 但是在激波等间断处会产生伪振荡, 因此, 考虑熵稳定格式的保平衡性.

熵稳定格式的表达形式可整理为

$\begin{matrix}\frac{{\rm d}}{{\rm d}t} {U}_{i}= -\frac{1}{\triangle x}( {F}^{ES1}_{i+1/2}- {F}^{ES1}_{i-1/2})-\frac{g}{2\triangle x} \left[{\begin{array}{c}0 \\\overline{h}_{i+1/2}[b]_{i+1/2}+\overline{h}_{i-1/2}[b]_{i-1/2}\\\end{array}}\right].\end{matrix}$

引理4.2[10] 熵稳定格式(4.4)具有保平衡性, 可保持离散的静水稳态(4.3).

引理4.1和引理4.2的具体证明过程在文献[10]中已详细给出, 这里不再赘述.

4.3 高阶熵稳定格式的保平衡性

本节将证明3.3节中构造的高阶熵稳定格式的保平衡性, 证明将分成两部分论述, 首先给出任意偶数阶熵守恒格式的保平衡性定理, 然后再将其推广至熵稳定格式.

4.3.1 任意偶数阶熵守恒格式的保平衡性

定理4.1 在恰当的源项离散形式下, 如果数值通量具有(3.5)式的形式,那么离散格式(3.1)是保平衡的.

以四阶熵守衡通量为例来证明熵守恒格式的保平衡性, 为方便起见, 记

$\begin{matrix}\widetilde{{F}}^{4}_{i+1/2}=\frac{4}{3}\widetilde{{F}}_{1,i+1/2}-\frac{1}{6}(\widetilde{{F}}_{2,i+1/2}+\widetilde{{F}}_{3,i+1/2}),\end{matrix}$

其中, $\widetilde{{F}}_{j,i+1/2}=\left[{\begin{array}{c}\overline{h}_{j,i+1/2}\overline{u}_{j,i+1/2}\\\frac{g}{2}\overline{h}^2_{j,i+1/2}+\overline{h}_{j,i+1/2}(\overline{u}_{j,i+1/2})^2\\\end{array}}\right]\ ( j=1,2,3 )$.

注意, 为了证明保平衡性, 这里需要对源项的离散方式进行修正, 与通量的构造方式相对应,

$\begin{matrix}\widetilde{{S}}_{i+1/2}&=&\frac{4}{3}\widetilde{{S}}(U_{i},U_{i+1})-\frac{1}{6}\left(\widetilde{{S}}(U_{i-1},U_{i+1})+\widetilde{{S}}(U_{i},U_{i+2})\right)\\&=&\frac{4}{3}\widetilde{{S}}_{1,i+1/2}-\frac{1}{6}(\widetilde{{S}}_{2,i+1/2}+\widetilde{{S}}_{3,i+1/2}),\end{matrix}$

其中, $\widetilde{{S}}_{j,i+1/2}=\left[{\begin{array}{c}0 \\\frac{g}{2\triangle x}(\overline{h}_{j,i+1/2}[b]_{j,i+1/2}+\overline{h}_{j,i-1/2}[b]_{j,i-1/2})\\\end{array}}\right]\ ( j=1,2,3 )$.

具有如上数值通量及源项离散形式(4.5)-(4.6)的有限体积格式为

$\begin{matrix}\frac{{\rm d}}{{\rm d}t}h_{i}&=&-\frac{1}{\triangle x}{\left\{\begin{array}{ccc} \left [\frac{4}{3}\overline{h}_{1,i+1/2}\overline{u}_{1,i+1/2}-\frac{1}{6}(\overline{h}_{2,i+1/2}\overline{u}_{2,i+1/2}+\overline{h}_{3,i+1/2}\overline{u}_{3,i+1/2})\right] \\[3mm] -\left [\frac{4}{3}\overline{h}_{1,i-1/2}\overline{u}_{1,i-1/2}-\frac{1}{6}(\overline{h}_{2,i-1/2}\overline{u}_{2,i-1/2}+\overline{h}_{3,i-1/2}\overline{u}_{3,i-1/2})\right]\end{array}\right\}},\\\frac{{\rm d}}{{\rm d}t}(h_{i}u_{i})&=&\frac{-1}{\triangle x}{\left\{\begin{array}{ccc} \left [\frac{4}{3}\left(\frac{g}{2}\overline{h}^2_{1,i+1/2}+\overline{h}_{1,i+1/2}(\overline{u}_{1,i+1/2})^2\right)\right.\\[3mm] \left. -\frac{1}{6}\left(\begin{array}{ccc} \frac{g}{2}\overline{h}^2_{2,i+1/2}+\overline{h}_{2,i+1/2}(\overline{u}_{2,i+1/2})^2\\[3mm] \frac{g}{2}\overline{h}^2_{3,i+1/2}+\overline{h}_{3,i+1/2}(\overline{u}_{3,i+1/2})^2\end{array}\right)\right]\\[8mm] -\left [\frac{4}{3}\left(\frac{g}{2}\overline{h}^2_{1,i-1/2}+\overline{h}_{1,i-1/2}(\overline{u}_{1,i-1/2})^2\right)\right.\\[3mm] \left.-\frac{1}{6}\left(\begin{array}{ccc} \frac{g}{2}\overline{h}^2_{2,i-1/2}+\overline{h}_{2,i-1/2}(\overline{u}_{2,i-1/2})^2\\[3mm] \frac{g}{2}\overline{h}^2_{3,i-1/2}+\overline{h}_{3,i-1/2}(\overline{u}_{3,i-1/2})^2\end{array}\right)\right]\end{array}\right\}}\\&&-\frac{g}{2\triangle x}{\left\{\begin{array}{ccc} \left [\frac{4}{3}(\overline{h}_{1,i+1/2}[b]_{1,i+1/2})-\frac{1}{6}(\overline{h}_{2,i+1/2}[b]_{2,i+1/2}+\overline{h}_{3,i+1/2}[b]_{3,i+1/2})\right]\\[3mm] +\left [\frac{4}{3}(\overline{h}_{1,i-1/2}[b]_{1,i-1/2})-\frac{1}{6}(\overline{h}_{2,i-1/2}[b]_{2,i-1/2}+\overline{h}_{3,i-1/2}[b]_{3,i-1/2})\right]\end{array}\right\}}.\end{matrix}$

根据引理4.1, 初始状态 $u_i \equiv 0, h_i+b_i \equiv \mbox{常数}$, 可得 $\frac{\rm d}{{\rm d}t}h_{i}\equiv 0$, 而

$\begin{matrix}\frac{{\rm d}}{{\rm d}t}(h_{i}u_{i})&=&-\frac{g}{\triangle x}\left[ \frac{2}{3}(\overline{h}^2_{1,i+1/2}-\overline{h}^2_{1,i-1/2})- \frac{1}{12}(\overline{h}^2_{2,i+1/2}-\overline{h}^2_{2,i-1/2})-\frac{1}{12}(\overline{h}^2_{3,i+1/2}-\overline{h}^2_{3,i-1/2})\right]\\&&-\frac{g}{\triangle x}\left[{\begin{array}{l} \frac{2}{3}(\overline{h}_{1,i+1/2}[b]_{1,i+1/2}+ \overline{h}_{1,i-1/2}[b]_{1,i-1/2})\\[3mm] -\frac{1}{12}(\overline{h}_{2,i+1/2}[b]_{2,i+1/2}+ \overline{h}_{2,i-1/2}[b]_{2,i-1/2})\\[3mm] -\frac{1}{12}(\overline{h}_{3,i+1/2}[b]_{3,i+1/2}+\overline{h}_{3,i-1/2}[b]_{3,i-1/2})\end{array}}\right],\end{matrix}$

根据恒等式

$\overline{h}^2_{j,i+1/2}-\overline{h}^2_{j,i-1/2}\equiv \overline{h}_{j,i+1/2}[h]_{j,i+1/2}+\overline{h}_{j,i-1/2}[h]_{j,i-1/2},$

进一步地, 有

$\begin{matrix}\mbox{上式}&=&\frac{-g}{\triangle x}\left[{\begin{array}{l} \frac{2}{3}(\overline{h}_{1,i+1/2}[h]_{1,i+1/2}+\overline{h}_{1,i-1/2}[h]_{1,i-1/2})\\[3mm] -\frac{1}{12}(\overline{h}_{2,i+1/2}[h]_{2,i+1/2}+\overline{h}_{2,i-1/2}[h]_{2,i-1/2})\\[3mm] -\frac{1}{12}(\overline{h}_{3,i+1/2}[h]_{3,i+1/2}+\overline{h}_{3,i-1/2}[h]_{3,i-1/2})\end{array}}\right]\\&&-\frac{g}{\triangle x}\left[{\begin{array}{l} \frac{2}{3}(\overline{h}_{1,i+1/2}[b]_{1,i+1/2}+ \overline{h}_{1,i-1/2}[b]_{1,i-1/2})\\[3mm] -\frac{1}{12}(\overline{h}_{2,i+1/2}[b]_{2,i+1/2}+ \overline{h}_{2,i-1/2}[b]_{2,i-1/2})\\[3mm]-\frac{1}{12}(\overline{h}_{3,i+1/2}[b]_{3,i+1/2}+ \overline{h}_{3,i-1/2}[b]_{3,i-1/2})\end{array}}\right]\\&=&-\frac{g}{\triangle x}\left[{\begin{array}{l} \frac{2}{3}(\overline{h}_{1,i+1/2}[h]_{1,i+1/2}+\overline{h}_{1,i+1/2}[b]_{1,i+1/2}\\[2mm]+\overline{h}_{1,i-1/2}[h]_{1,i-1/2}+ \overline{h}_{1,i-1/2}[b]_{1,i-1/2})\\[2mm]-\frac{1}{12}(\overline{h}_{2,i+1/2}[h]_{2,i+1/2}+\overline{h}_{2,i+1/2}[b]_{2,i+1/2}\\[2mm]+\overline{h}_{2,i-1/2}[h]_{2,i-1/2}+ \overline{h}_{2,i-1/2}[b]_{2,i-1/2})\\[2mm]-\frac{1}{12}(\overline{h}_{3,i+1/2}[h]_{3,i+1/2}+\overline{h}_{3,i+1/2}[b]_{3,i+1/2}\\[2mm]+\overline{h}_{3,i-1/2}[h]_{3,i-1/2}+ \overline{h}_{3,i-1/2}[b]_{3,i-1/2})\end{array}}\right]\\&=&-\frac{g}{\triangle x}\left[{\begin{array}{l}\frac{2}{3}(\overline{h}_{1,i+1/2}[h+b]_{1,i+1/2}+ \overline{h}_{1,i-1/2}[h+b]_{1,i-1/2})\\[3mm]-\frac{1}{12}(\overline{h}_{2,i+1/2}[h+b]_{2,i+1/2}+ \overline{h}_{2,i-1/2}[h+b]_{2,i-1/2})\\[3mm]-\frac{1}{12}(\overline{h}_{3,i+1/2}[h+b]_{3,i+1/2}+ \overline{h}_{3,i-1/2}[h+b]_{3,i-1/2})\end{array}}\right],\end{matrix}$

又由于对任意的$i \in \mathbb{Z}$, 都有$[h+b]_{j,i+1/2}=0\ (j=1,2,3)$, 因此$\frac{{\rm d}}{{\rm d}t}(h_{i}u_{i})\equiv 0$.证毕.

4.3.2 高阶熵稳定格式的保平衡性

上节证明了数值通量采用四阶熵守恒通量时熵守恒格式的保平衡性, 用类似的方法可以证明六阶乃至任意偶数阶格式有同样性质, 下面考虑高阶熵稳定格式的保平衡性.

定理4.2 高阶熵稳定格式(3.9)保持离散的静水稳态, 具有保平衡性.

初始数据若满足对任意的整数 $i$ 都有$u_{i}=0,h_{i}+b_{i}=\mbox{常数}$, 则有$[u]_{i+1/2}\equiv 0,$$[h+b]_{i+1/2}\equiv 0$. 而重构在单元交界面$x_{i+1/2}$处的重构跳跃值为

$\begin{matrix}{\langle {V} \rangle}_{i+1/2}&=&{[V]}^{+}_{i+1/2}-{[V]}^{-}_{i+1/2}\\&=&{[V]}_{i+1/2}+\left[{\begin{array}{l} \widetilde{\omega}_{4}\left(\frac{1}{15}{[V]}_{i+5/2}- \frac{1}{5}{[V]}_{i+3/2}-\frac{7}{15}{[V]}_{i+1/2}+\frac{1}{10}{[V]}_{i-1/2}\right)\\[3mm] -\frac{{[V]}_{i+1/2}}{2}\widetilde{\omega}_{1}-\widetilde{\omega}_{2}\left(\frac{1}{12}{[V]}_{i+3/2}+\frac{5}{12}{[V]}_{i+1/2}\right)-\widetilde{\omega}_{3}\left(\frac{1}{2}{[V]}_{i+3/2}\right)\end{array}}\right]\\&&-\left[{\begin{array}{l} \omega_{4}\left(-\frac{1}{10}{[V]}_{i+3/2}+ \frac{7}{15}{[V]}_{i+1/2}+\frac{1}{5}{[V]}_{i-1/2}-\frac{1}{15}{[V]}_{i-3/2}\right)\\[3mm]+\omega_{1}\left(\frac{1}{2}{[V]}_{i-1/2}\right)+\omega_{2}\left(\frac{5}{12}{[V]}_{i+1/2}+\frac{1}{12}{[V]}_{i-1/2}\right)+\omega_{3}\left(\frac{1}{2}{[V]}_{i+1/2}\right)\end{array}}\right],\end{matrix}$

通过熵变量${V}$的定义, 对任意的$i \in \mathbb{Z}$也有${[V]}_{i+1/2}=0$成立, 因此有${\langle {V} \rangle}_{i+1/2} \equiv 0$成立, 从而去掉(3.8)式中的耗散项部分, 熵稳定格式转换成熵守恒格式, 再由定理4.1, 可以得到结论 $\frac{{\rm d}}{{\rm d}t}h_{i}\equiv 0,\frac{{\rm d}}{{\rm d}t}(h_{i}u_{i})\equiv 0$, 故高阶熵稳定格式(3.9)具有保平衡性.证毕.

5 数值算例

下面用几个典型的一维和二维算例来验证格式的保平衡性. 分别用具有熵通量(3.4)的低阶熵稳定格式以及高阶熵稳定格式(3.9)来求解稳态解的小扰动问题. 低阶熵稳定格式的熵守恒部分为二阶熵守恒通量, 耗散项为一阶, 而高阶熵稳定格式具体为熵守恒部分采用四阶熵守恒通量, 耗散项部分采用五阶的WENO-AO重构, 低阶格式和高阶格式的总体精度理论上应分别为一阶和四阶, 两种格式在下文中分别用ES1以及HES 表示. 时间方向采用四阶强稳定Runge-Kutta法[18], 根据算例进一步验证格式的精度及性能.

算例1 首先用平底地形(即$b$是一个常数)来测试格式的精度, 考虑具有如下初值的溃坝问题

$\begin{matrix}h(x,0)=\left\{\begin{array}{lll}2, \quad \enspace x<0,\\1.5, \quad x \ge 0,\\\end{array}\right. \quad u(x,0)=0, \quad x \in [-1,1],\end{matrix}$

不妨取$b=0$, 计算$t=0.4$时刻的解.

该问题的精确解是由左行稀疏波和右行激波所构成, 为了避免间断处数据对精度的影响, 选取光滑区域即$(-1,-0.6)$间的数据来求解格式精度, 由表1可看出高阶熵稳定格式的精度可达到四阶, 而具有一阶耗散的熵稳定格式近似为一阶精度. 图1(a)分别利用ES1格式以及HES格式计算网格数为$N=100$时的水位$\omega (b=0,$$h=\omega)$, 可以看出高阶熵稳定格式能够更好地逼近精确解. 图1(b) 则给出了两种格式的总熵随时间的变化情况, 高阶格式比低阶的熵稳定格式产生的耗散少.

表1   浅水波光滑区域的误差及收敛阶

新窗口打开| 下载CSV


图1

图1   具有平底地形溃坝问题的数值结果


算例2 初始数据是扰动平稳解

$\begin{matrix}\omega (x,0)=1+\epsilon,\quad u(x,0)=0, \quad h(x,0)=\omega (x,0)-b(x),\end{matrix}$

这里的$\epsilon$在区间$0.1<x<0.2$内是非零常数, $x \in (0,1)$, 底部地形包括一个驼峰

$\begin{matrix}b(x)=\left\{\begin{array}{lll}0.25(\cos(10\pi (x-0.5))+1), \quad& |x-0.5| <0.1,\\0,\qquad \qquad \qquad \qquad \qquad \quad \enspace \quad & \mbox{其他}.\\\end{array}\right.\end{matrix}$

该算例选自文献[5], 重力加速度$g=1$, 计算$t=0.7$时刻解的情况.

假定$\epsilon=0.2$, 图2给出了$t=0$以及$t=0.7$时刻水位和底部地势的地形图, 可以清晰地观察到水位随时间推进的变化情况. 图3利用ES1格式以及HES格式分别计算了网格数为$N=100\mbox{和}200$时的水位$\omega $, 由图形可看出格式能准确求解驼峰附近$(0.4<x<0.6)$的解, 且网格数越大, 数值解与精确解的图像吻合度就越高, 对扰动的传播情况捕捉地更为准确. 图4给出了$\epsilon=0.01$, 网格数$N=150$时两种格式求出的$\omega $结果以及总熵变化情况, 当$\epsilon$ 由0.2减小至0.01后, 在驼峰附近仍无伪振荡产生. 由图4(b)可见两个格式的总熵变化量基本相同, 但由水位的变化图可得HES格式有效地改善了抹平现象.

图2

图2   不同时刻水位和底部的地形图


图3

图3   不同网格数下水位$\omega $的结果


图4

图4   $\epsilon=0.01, N=150$情况下的数值结果


算例3 考虑静止湖泊扰动问题, 湖泊存在一个非常小的扰动, 试图研究这种扰动如何随时间传播, 初始状态为

$\begin{matrix}u(x,0) \equiv 0,\quad h(x,0)=\left\{\begin{array}{lll}1.01-b(x), \quad &|x-6|<0.25,\\1-b(x), \quad \quad \; &\mbox{其他},\\\end{array}\right. \quad x \in [020].\end{matrix}$

底部地形是一个抛物形的凸起

$\begin{matrix}b(x)=\left\{\begin{array}{lll} \frac{4-(x-10)^2}{20}, \quad &|x-10|<2,\\0, \qquad\qquad \quad \quad \; &\mbox{其他},\\\end{array}\right.\quad x \in [020].\end{matrix}$

采用Neumann边界条件, 重力加速度$g=9.812$, 计算$t=1.5$时刻水位$\omega=h+b$的变化情况.

图5(a)给出了$t=1.5$时刻水位$\omega$以及底部地形$b$的整体地形图, 熵稳定格式能够很好地保持稳态, 进一步验证了格式的保平衡性. 图5(b)则具体利用上述两种格式描绘了$t=1.5$ 时刻水面的波动, 对照参考解, 高阶熵稳定格式能够相对精确地刻画扰动的传播过程.

图5

图5   算例3的数值结果


算例4 考虑一维带源项浅水波方程的稳态问题, 底部地形

$\begin{matrix}b(x)=\left\{\begin{array}{lll} \frac{1}{4}(\cos(10\pi(x-1.5))+1), \quad &1.4\le x\le1.6,\\0, \qquad\qquad\qquad\qquad\qquad \quad \; &\mbox{其他},\\\end{array}\right.\end{matrix}$

计算区域$[02]$, 初始状态为

$\begin{matrix} h(x,0)=\left\{\begin{array}{lll}1-b(x) +\theta, &\quad 1.1\le x\le1.6,\\1-b(x), \quad \ \ &\quad \mbox{其他},\\\end{array}\right. \qquad hu(x,0)\equiv 0,\end{matrix}$

参数$\theta=0.001$, $CFL$取为$0.45$, $g=0.98$, 计算$t=0.3$时刻该微小扰动问题的解.由两种熵稳定格式(ES1, HES)计算出的数值解见图6. 正如预期的那样, 图6(a)体现了熵稳定格式能够保持湖泊静止稳态. 图6(b)显示了$h+b$的变化情况, 两种格式均能准确捕捉光滑解, 但在激波等间断区域, 低阶格式的解会有所抹平, 而高阶熵稳定格式的数值解与参考解的拟合效果更好, 且有效抑制了非物理现象的产生.

图6

图6   算例4的数值结果


算例5 二维湖泊静止问题, 底部地势是一个椭圆形凸起

$\begin{matrix}b(x,y)=0.8e^{-5(x-0.9)^2-50(y-0.5)^2},\end{matrix}$

其中$(x,y)\in [02]\times[01]$, 在平坦的水平面上增加一个小扰动

$\begin{matrix} \omega(x,y,0)=\left\{\begin{array}{lll}1+0.01, &\quad 0.05\le x\le 0.15,\\1, \qquad \quad&\quad \mbox{其他}.\\\end{array}\right.\end{matrix}$

该算例选自文献[19],$g=9.812$,分别在$100 \times 50$以及$200 \times 100$的网格上计算$t=0.6, 0.9, 1.2,$$ 1.5, 1.8$时刻扰动的传播过程. 图7描绘了二维静止湖泊的水平面和底部地形图. 该问题的精确解涉及两个波, 一个左行波, 一个右行波, 尤其当扰动经过底部凸起时, 波的中间部分速度比其他地方慢, 从而导致初始平面扰动发生复杂变化, 图8用高阶熵稳定格式比较了两种粗细网格下水位$\omega$在五个不同时刻的扰动演化过程, 对比文献[19]的数值结果, 这里用更少的网格数达到了不错的效果, 避免了由于网格数过大而造成的计算困难.

图7

图7   水位$\omega$及底部地势$b$的地形图


图8

图8   分别在$100\times50$$200\times100$网格上计算五个时刻扰动的传播情况


6 结论

本文针对具有底部地形源项的浅水波方程设计了一类高阶熵稳定格式, 通过对耗散项部分采用WENO-AO 重构, 确保了格式在实现高精度的同时也满足了保平衡性, 也可以选用其他重构, 例如ENO、WENO重构等. 文中研究的系统主要考虑了底部地形源项, 而底部或者表面的摩擦以及河道宽度的变化等因素还有待进一步研究, 浅水波方程的源项处理方法还有很多进一步研究的内容.

稳态问题是研究浅水波方程的重要课题, 在稳态情况下, 浅水波方程存在非零通量由源项平衡的稳态解, 因此, 所构造的数值格式应该精确捕捉稳态解及其扰动, 文中由保平衡性来刻画这一特点, 证明了高阶熵稳定格式的保平衡性. 稳态包括静水稳态以及动水稳态, 本文主要研究的是静水稳态, 有关动水稳态[20]的性质将在后面继续研究. 最后利用特殊的一维和二维算例进行了数值实验, 旨在说明格式能够很好地处理稳态解的小扰动问题, 但在实际应用中, 常常会面临一些更复杂的情形, 例如, 具有干湿界面的浅水波问题(模拟岛屿周围的水面变化), 以及部分淹没情况下的稳态小扰动问题[6], 对以上特殊情形需要考虑水深$h$的非负性, 一旦$h$产生负值, 计算过程将会中断, 根据不同的情况对所构造的数值格式进行恰当处理是今后工作的重点.

参考文献

Lin G F, Lai J S, Guo W D.

Performance of high-resolution TVD schemes for $1$D dam-break simulations

J Chin Inst Eng, 2005, 28(5): 771-782

DOI:10.1080/02533839.2005.9671047      URL     [本文引用: 1]

Xing Y L, Shu C W.

High order finite difference WENO schemes with the exact conservation property for the shallow water matrixs

J Comput Phys, 2005, 208(1): 206-227

DOI:10.1016/j.jcp.2005.02.006      URL     [本文引用: 1]

Kesserwani G, Liang Q H.

A conservative high-order discontinuous Galerkin method for the shallow water matrixs with arbitrary topography

Int J Numer Meth Eng, 2011, 86(1): 47-69

DOI:10.1002/nme.v86.1      URL     [本文引用: 1]

郑素佩, 李霄, 赵青宇, .

求解二维浅水波方程的旋转混合格式

应用数学和力学, 2022, 43(2): 176-186

[本文引用: 1]

Zheng S P, Li X, Zhao Q Y, et al.

A rotated mixed scheme for solving 2D shallow water matrixs

Appl Math Mech, 2022, 43(2): 176-186

[本文引用: 1]

Kurganov A, Levy D.

Central-upwind schemes for the Saint-Venant system

Math Model Anal, 2002, 36(3): 397-425

[本文引用: 2]

Liu X, Albright J, Epshteyn Y, et al.

Well-balanced positivity preserving central-upwind scheme with a novel wet/dry reconstruction on triangular grids for the Saint-Venant system

J Comput Phys, 2018, 374: 213-236

DOI:10.1016/j.jcp.2018.07.038      URL     [本文引用: 2]

Kurganov A, Noelle S, Petrov A G.

Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi matrixs

SIAM J Sci Comput, 2001, 23(3): 707-740

DOI:10.1137/S1064827500373413      URL     [本文引用: 1]

董建.

中心格式在 Saint-Venant 方程组上的应用研究

数学物理学报, 2019, 39A(2): 372-385

[本文引用: 1]

Dong J.

Research on the application of central scheme in Saint-Venant system

Acta Math Sci, 2019, 39A(2): 372-385

[本文引用: 1]

Tadmor E.

The numerical viscosity of entropy stable schemes for systems of conservation laws

I. Math Comput, 1987, 49(179): 91-103

[本文引用: 1]

Fjordholm U S, Mishra S, Tadmor E.

Well-balanced and energy stable schemes for the shallow water matrixs with discontinuous topography

J Comput Phys, 2011, 230(14): 5587-5609

DOI:10.1016/j.jcp.2011.03.042      URL     [本文引用: 7]

张海军, 封建湖, 程晓晗, .

带源项浅水波方程的高分辨率熵稳定格式

应用数学和力学, 2018, 39(8): 935-945

[本文引用: 2]

Zhang H J, Feng J H, Cheng X H, et al.

An entropy stable scheme for shallow water matrixs with source terms

Appl Math Mech, 2018, 39(8): 935-945

[本文引用: 2]

Fjordholm U S, Mishra S, Tadmor E.

Arbitrarily high-order accurate entropy stable essentially non-oscillatory schemes for systems of conservation laws

SIAM J Numer Anal, 2012, 50(2): 544-573

DOI:10.1137/110836961      URL     [本文引用: 3]

郑素佩, 建芒芒, 封建湖, .

保号WENO-AO型中心迎风格式

计算物理, http://kns.cnki.net/kcms/detail/11.2011.O4.20220322.1639.008.html

URL     [本文引用: 3]

Zheng S P, Jian M M, Feng J H, et al.

Sign preserving WENO-AO-type central upwind schemes

Chinese J Comput Phys, http://kns.cnki.net/kcms/detail/11.2011.O4.20220322.1639.008.html

URL     [本文引用: 3]

Lax P D.

Weak solutions of nonlinear hyperbolic matrixs and their numerical computation

Commun Pur Appl Math, 1954, 7(1): 159-193

DOI:10.1002/(ISSN)1097-0312      URL     [本文引用: 1]

Fjordholm U S, Ray D.

A sign preserving WENO reconstruction method

J Sci Comput, 2016, 68(1): 42-63

DOI:10.1007/s10915-015-0128-y      URL     [本文引用: 1]

Shu C W.

Essentially non-oscillatory and weighted essentially non-oscillatory schemes

Acta Numer, 2020, 29: 701-762

DOI:10.1017/S0962492920000057      URL     [本文引用: 1]

Essentially non-oscillatory (ENO) and weighted ENO (WENO) schemes were designed for solving hyperbolic and convection–diffusion equations with possibly discontinuous solutions or solutions with sharp gradient regions. The main idea of ENO and WENO schemes is actually an approximation procedure, aimed at achieving arbitrarily high-order accuracy in smooth regions and resolving shocks or other discontinuities sharply and in an essentially non-oscillatory fashion. Both finite volume and finite difference schemes have been designed using the ENO or WENO procedure, and these schemes are very popular in applications, most noticeably in computational fluid dynamics but also in other areas of computational physics and engineering. Since the main idea of the ENO and WENO schemes is an approximation procedure not directly related to partial differential equations (PDEs), ENO and WENO schemes also have non-PDE applications. In this paper we will survey the basic ideas behind ENO and WENO schemes, discuss their properties, and present examples of their applications to different types of PDEs as well as to non-PDE problems.

郑素佩, 徐霞, 封建湖, .

高阶保号熵稳定格式

数学物理学报, 2021, 41A(5): 1296-1310

[本文引用: 1]

Zheng S P, Xu X, Feng J H, et al.

High order sign preserving stable schemes

Acta Math Sci, 2021, 41A(5): 1296-1310

[本文引用: 1]

徐霞.

求解双曲守恒律方程的任意阶保号熵稳定格式研究

西安:长安大学, 2021

[本文引用: 1]

Xu X.

The Research of Arbitrary High-Order Sign Prserving Entropy Stable Schemes for Hyperbolic Conservation Laws

Xi'an: Chang'an University, 2021

[本文引用: 1]

Leveque R J.

Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm

J Comput Phys, 1998, 146(1): 346-365

DOI:10.1006/jcph.1998.6058      URL     [本文引用: 2]

罗一鸣, 李订芳, 刘敏, .

一种维持 Saint-Venant 方程组移动稳态解的中心格式

数学物理学报, 2022, 42A(3): 891-903

[本文引用: 1]

Luo Y M, Li D F, Liu M, et al.

Moving-water equilibria preserving central scheme for the Saint-Venant system

Acta Math Sci, 2022, 42A(3): 891-903

[本文引用: 1]

/