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

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

罗一鸣^,¹, 李订芳^,¹, 刘敏^,¹, 董建^,²

¹ 武汉大学数学与统计学院武汉 430072

² 国防科技大学文理学院长沙 410073

Moving-Water Equilibria Preserving Central Scheme for the Saint-Venant System

Luo Yiming^,¹, Li Dingfang^,¹, Liu Min^,¹, Dong Jian^,²

¹ School of Mathematics and Statistics, Wuhan University, Wuhan 430072

² College of Liberal Arts and Sciences, National University of Defense Technology, Changsha 410073

通讯作者: 董建, E-mail: j.dong@whu.edu.cn

收稿日期: 2021-01-11

基金资助:

国家重点研发计划项目. 2017YFC0405901

Received: 2021-01-11

Fund supported:

the National Key Research and Development Project. 2017YFC0405901

作者简介 About authors

罗一鸣,E-mail:luoyiming@whu.edu.cn , E-mail：luoyiming@whu.edu.cn

李订芳,E-mail:dfli@whu.edu.cn , E-mail：dfli@whu.edu.cn

刘敏,E-mail:liumin@whu.edu.cn , E-mail：liumin@whu.edu.cn

Abstract

In this paper, we propose a second-order unstaggered central finite volume scheme for the Saint-Venant system. Classical central scheme can preserve still-water steady state solution by reconstructing conservative variables and the water level, but generates enormous numerical oscillation when considering moving-water steady state. We choose to reconstruct conservative variables and the energy, and design a new discretization of the source term to preserve moving-water equilibria and capture its small perturbations. In the end, several classical numerical experiments are performed to verify the proposed scheme which is convergent, well-balanced and robust.

Keywords： Saint-Venant system ; Unstaggered central scheme ; Moving-water equilibria

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

本文引用格式

罗一鸣, 李订芳, 刘敏, 董建. 一种维持Saint-Venant方程组移动稳态解的中心格式. 数学物理学报[J], 2022, 42(3): 891-903 doi:

Luo Yiming, Li Dingfang, Liu Min, Dong Jian. Moving-Water Equilibria Preserving Central Scheme for the Saint-Venant System. Acta Mathematica Scientia[J], 2022, 42(3): 891-903 doi:

1 引言

Saint-Venant方程组由Saint-Venant于1871年提出^[1], 至今仍在浅水流、涌浪传播、水库和明渠流等问题的数值求解中发挥着重要作用. 在一维情形下, 其向量形式可表示为

$\begin{equation} {{{U}}_t} + {{F}}{({{U}})_x} = {{S}}({{U}}), \end{equation}$

(1.1)

其中

$\begin{equation} {{U}} = \left( \begin{array}{c} h\\ q \end{array} \right), {\quad} {{F}}\left( {{U}} \right) = \left( \begin{array}{c} q\\ { } h{u^2} + \frac{g}{2}{h^2} \end{array} \right), {\quad} {{S}}\left( {{U}} \right) = \left( \begin{array}{c} 0\\ - gh{B_x} \end{array} \right).{\nonumber} \end{equation}$

(1.2)

$h\left( {x, t} \right)$ 表示水的深度, $u\left( {x, t} \right)$ 表示速度, $q\left( {x, t} \right) = h\left( {x, t} \right)u\left( {x, t} \right)$ 表示流量, $g$ 是重力加速度常数, $B\left( x \right)$ 表示底部地形. 初值条件 $h\left( {x, 0} \right)$ , $q\left( {x, 0} \right)$ 已知, 若求解区域并非整个实数域, 还需添加适当的边值条件.

该方程组的一个重要特征是其满足稳态解

$\begin{equation} q \equiv const, {\quad} E: = \frac{{{u^2}}}{2} + g\left( {h + B} \right) \equiv const.{\nonumber} \end{equation}$

(1.3)

其中 $E\left( {x, t} \right)$ 表示能量, $h\left( {x, t} \right) + B\left( x \right)$ 表示水位. 当上式满足 $q \equiv 0$ , $h + B \equiv const$ 时, 称之为“静稳态解”. 如果所设计的数值格式能够从离散角度维持这种稳态, 就称格式具有“和谐性”.

Nessyahu和Tadmor于1990年提出了用于求解双曲守恒律问题的二阶交错中心格式^[2], 相较于迎风格式, 其优势在于避免了求解黎曼问题, 但于相邻时间层, 数值解会在交错单元与非交错单元之间来回更替, 导致了复杂的边界条件设置. 为此, Touma提出了稳健的二阶非交错中心格式^[3], 并将其应用于Saint-Venant方程组, 提出了保正的能维持静稳态解的格式^[4], 但在干湿界面处, “保正性”与“和谐性”无法同时得到满足. Dong和Li^[5]通过重构守恒变量及水位值, 且在干湿界面处构造关于水位的映射, 使格式在干湿界面处经过“保正性”处理后仍可维持静稳态解. 然而在用中心格式求解移动稳态问题时, 在底部地形变化处( ${B_x} \ne 0$ ) 流量与能量会产生较为明显的虚假振荡.

由于能量与水深之间的非线性关系, 设计能维持移动稳态解的数值格式往往难度更大. 为此, Kurganov等人提出了中心-迎风格式^{[6, 7]}, Liu等人提出了部分松弛格式^[8], Shu等人提出了加权本质非振荡(WENO) 格式及间断有限元(DG)格式^[9-11]. 这些格式很好地避免了中心格式存在的问题. 更多结果参见文献[12-17].

本文提出了一种具有二阶精度的能够精确维持移动稳态解的非交错中心格式. 一方面, 基于移动稳态解的特征, 选择重构守恒变量以及能量值, 并通过求解非线性方程得到基于稳态时的水深; 另一方面, 在对守恒变量进行时间层的更新时, 基于稳态时的能量, 提出了一种新的源项离散方法, 使其能够平衡非零通量的梯度, 以保证格式的“和谐性”.

文章的结构如下, 第2节, 对格式做详细叙述, 包括经典格式的回顾, 新格式的构造以及“和谐性”的证明; 第3节, 通过一些经典的数值算例验证了格式的收敛性, 和谐性及稳健性; 第4节, 对格式做出总结.

2 维持移动稳态解的非交错中心格式

2.1 经典的非交错中心格式

针对Saint-Venant方程组的初值问题, 先简要介绍一下经典的非交错中心格式. 将计算区域 $\Omega = [a, b] \subset \mathbb{R}$ 做均匀网格剖分, 并定义非交错单元 ${C_i}: = [{x_{i - \frac{1}{2}}}, {x_{i + \frac{1}{2}}}]$ 以及交错单元 ${C_{i + \frac{1}{2}}}: = [{x_i}, {x_{i + 1}}]$ , 其中, ${x_{i + \frac{1}{2}}} = {x_i} + \frac{{\Delta x}}{2}$ , ${x_i} = a + \left( {i - \frac{1}{2}} \right)\Delta x$ , $\forall i \in Z$ . 不失一般性, 假设在 ${t^n} = n\Delta t$ 时刻, 数值解

$\overline {{U}} _i^n \approx \frac{1}{{\Delta x}}\int_{{x_{i - \frac{1}{2}}}}^{{x_{i + \frac{1}{2}}}} {{{U}}\left( {x, {t^n}} \right)} {\rm d}x{\nonumber}$

是已知的.

(1) 向前投影

为了获得二阶精度且避免虚假震荡, 对 ${{U}}$ 做分片线性重构

${{U}}\left( {x, {t^n}} \right) \approx {{{P}}_i}\left( x \right): = \overline {{U}} _i^n + \left( {x - {x_i}} \right){\left( {{{\overline {{U}} }_i}} \right)^\prime }, x \in {C_i},$

其中 ${\left( {{{\overline {{U}} }_i}} \right)^\prime }$ 表示对 ${{U}}$ 在点 $x = {x_i}$ 的导数的逼近, 用minmod函数^[2]来计算

${\left( {{{\overline {{U}} }_i}} \right)^\prime } = {\rm minmod} \left( {\theta \frac{{{{\overline {{U}} }_i} - {{\overline {{U}} }_{i - 1}}}}{{\Delta x}}, \frac{{{{\overline {{U}} }_{i + 1}} - {{\overline {{U}} }_{i - 1}}}}{{2\Delta x}}, \theta \frac{{{{\overline {{U}} }_{i + 1}} - {{\overline {{U}} }_i}}}{{\Delta x}}} \right),$

其中 $\theta \in [1, 2]$ 为已知参数

${\rm minmod} \left( {a, b, c} \right) = \left\{ \begin{array}{ll} \max \left( {a, b, c} \right), &\max \left( {a, b, c} \right) < 0, \\ \min \left( {a, b, c} \right), &\min \left( {a, b, c} \right) > 0, \\ 0, &\mbox{其他}. \end{array} \right.{\nonumber}$

利用上述分片线性函数 ${{{P}}_i}\left( x \right)$ , 去求其在交错单元 ${C_{i + \frac{1}{2}}}$ 的积分平均值, 并定义为 $\overline {{U}} _{i + \frac{1}{2}}^n$ :

$\begin{eqnarray*} \overline {{U}} _{i + \frac{1}{2}}^n& =& \frac{1}{{\Delta x}}\bigg( {\int_{{x_i}}^{{x_{i + \frac{1}{2}}}} {{{{P}}_i} \left( x \right){\rm d}x + \int_{{x_{i + \frac{1}{2}}}}^{{x_{i + 1}}} {{{{P}}_{i + 1}}\left( x \right){\rm d}x} } } \bigg) \\ &=& \frac{{\overline {{U}} _i^n + \overline {{U}} _{i + 1}^n}}{2} - \frac{{\Delta x}}{8}\left( {{{\left( {{{\overline {{U}} }_{i + 1}}} \right)}^\prime } - {{\left( {{{\overline {{U}} }_i}} \right)}^\prime }} \right).{\nonumber} \end{eqnarray*}$

(2) 预测步

将方程(1.1)在矩形 $[{x_i}, {x_{i + 1}}] \times [{t^n}, {t^{n + 1}}]$ 上积分并利用高斯定理得

$\begin{eqnarray} && \int_{{x_i}}^{{x_{i + 1}}} {{{U}}\left( {x, {t^{n + 1}}} \right)} {\rm d}x - \int_{{x_i}}^{{x_{i + 1}}} {{{U}}\left( {x, {t^n}} \right)} {\rm d}x + \int_{{t^n}}^{{t^{n + 1}}} {{{F}}\left( {{{U}}\left( {{x_{i + 1}}, t} \right)} \right)} {\rm d}t {}\\ & &- \int_{{t^n}}^{{t^{n + 1}}} {{{F}}\left( {{{U}}\left( {{x_i}, t} \right)} \right)} {\rm d}t = \int_{{x_i}}^{{x_{i + 1}}} {\int_{{t^n}}^{{t^{n + 1}}} {{{S}}\left( {{U}} \right)} } {\rm d}t{\rm d}x. \end{eqnarray}$

(2.1)

在 ${t^{n + 1}}$ 时刻 ${{U}}\left( {x, {t^{n + 1}}} \right)$ 也通过分片线性函数来近似

${{U}}\left( {x, {t^{n + 1}}} \right) \approx {{{P}}_{i + \frac{1}{2}}}\left( x \right): = \overline {{U}} _{i + \frac{1}{2}}^{n + 1} + \left( {x - {x_i}} \right){\left( {{{\overline {{U}} }_{i + \frac{1}{2}}}} \right)^\prime }, \; \; x \in {C_{i + \frac{1}{2}}}, {\nonumber}$

其中 ${\left( {{{\overline {{U}} }_{i + \frac{1}{2}}}} \right)^\prime }$ 表示对 ${{U}}$ 在点 $x = {x_{i + \frac{1}{2}}}$ 的导数的逼近. 对(2.1)式利用中点积分公式求通量的时间积分, 得到其离散形式

$\begin{equation} \overline {{U}} _{i + \frac{1}{2}}^{n + 1} = \overline {{U}} _{i + \frac{1}{2}}^n - \frac{{\Delta t}}{{\Delta x}}\left( {{{F}}\left( {\overline {{U}} _{i + 1}^{n + \frac{1}{2}}} \right) - {{F}}\left( {\overline {{U}} _i^{n + \frac{1}{2}}} \right)} \right) + \frac{1}{{\Delta x}}\int_{{x_i}}^{{x_{i + 1}}} {\int_{{t^n}}^{{t^{n + 1}}} {{{S}}\left( {{U}} \right)} } {\rm d}t{\rm d}x. \end{equation}$

(2.2)

由(2.2)式即可得 ${t^{n + 1}}$ 时刻非定常项在交错单元的积分平均值, 其中预测值 $\overline {{U}} _i^{n + \frac{1}{2}}$ 通过对时间的一阶泰勒展开得到

$\begin{equation} \overline {{U}} _i^{n + \frac{1}{2}} = \overline {{U}} _i^n + \frac{{\Delta t}}{2}\left( { - {{\left( {{{{F}}_i}} \right)}^\prime } + {{S}}_i^n} \right). \end{equation}$

(2.3)

${\left( {{{{F}}_i}} \right)^\prime }$ 表示对通量 ${{F}}\left( {{U}} \right)$ 在点 $x = {x_i}$ 的导数的逼近, ${{S}}_i^n \approx {{S}}\left( {{x_i}, \overline {{U}} _i^n} \right)$ 表示源项的逼近. 为使格式具有稳定性, 在时间步长的选取上应满足CFL条件

$\frac{{\Delta t}}{{\Delta x}} \times \mathop {\max }\limits_i \left| {u_i^n + \sqrt {gh_i^n} } \right| \leqslant \beta ,$

其中, $\beta$ 为CFL数, 通常取 $\beta \leqslant 1$ . 具体证明可参考文献[2].

(3) 向后投影

同理于向前投影, 可求得 ${t^{n + 1}}$ 时刻非定常项在非交错单元的积分平均值

$\overline {{U}} _i^{n + 1} = \frac{{\overline {{U}} _{i - \frac{1}{2}}^{n + 1} + \overline {{U}} _{i + \frac{1}{2}}^{n + 1}}}{2} - \frac{{\Delta x}}{8}\left( {{{\left( {{{\overline {{U}} }_{i + \frac{1}{2}}}} \right)}^\prime } - {{\left( {{{\overline {{U}} }_{i - \frac{1}{2}}}} \right)}^\prime }} \right).{\nonumber}$

2.2 新的非交错中心格式

(1) 底部函数处理

将底部地形函数 $B\left( x \right)$ 用分片线性函数 $\widetilde B\left( x \right)$ 近似^[12]:

$\widetilde B\left( x \right) = {B_{i - \frac{1}{2}}} + \left( {{B_{i + \frac{1}{2}}} - {B_{i - \frac{1}{2}}}} \right)\frac{{x - {x_{i - \frac{1}{2}}}}}{{\Delta x}}, x \in {C_i}.{\nonumber}$

其中

${B_{i + \frac{1}{2}}}: = \frac{{B\left( {{x_{i + \frac{1}{2}}} + 0} \right) + B\left( {{x_{i + \frac{1}{2}}} - 0} \right)}}{2}.{\nonumber}$

若 $B\left( x \right)$ 在 $x = {x_{i + \frac{1}{2}}}$ 处连续, 则有 ${B_{i + \frac{1}{2}}} = B\left( {{x_{i + \frac{1}{2}}}} \right)$ . 再将 $\widetilde B\left( x \right)$ 在非交错单元 ${C_i}$ 的积分平均值定义为 ${B_i}$

${B_i}: = \frac{1}{{\Delta x}}\int_{{x_{i - \frac{1}{2}}}}^{{x_{i + \frac{1}{2}}}} {\widetilde B\left( x \right){\rm d}x = \widetilde B\left( {{x_i}} \right)} = \frac{{{B_{i - \frac{1}{2}}} + {B_{i + \frac{1}{2}}}}}{2}.{\nonumber}$

(2) 向前投影

为维持移动稳态解, 首先计算每个非交错单元 ${C_i}$ 的能量 $E$ :

${E_i} = \frac{{{u_i}^2}}{2} + g\left( {{h_i} + {B_i}} \right).{\nonumber}$

对 ${{V}} = {\left( {h, q, E} \right)^{\rm T}}$ , 而非仅对 ${{U}}$ 做分片线性重构, 得到 ${h_{i + \frac{1}{2}}}$ , ${q_{i + \frac{1}{2}}}$ 及 ${E_{i + \frac{1}{2}}}$ , 并定义 ${\widehat E_{i + \frac{1}{2}}}$ 如下

$\begin{equation} {\widehat E_{i + \frac{1}{2}}}: = \frac{{{u_{i + \frac{1}{2}}}^2}}{2} + g\left( {{h_{i + \frac{1}{2}}} + {{\overline B }_{i + \frac{1}{2}}}} \right), \end{equation}$

(2.4)

其中, ${\overline B _{i + \frac{1}{2}}}$ 为 $\widetilde B\left( x \right)$ 在交错单元 ${C_{i + \frac{1}{2}}}$ 的积分平均值

${\overline B _{i + \frac{1}{2}}} = \frac{1}{{\Delta x}}\int_{{x_i}}^{{x_{i + 1}}} {\widetilde B\left( x \right){\rm d}x = \frac{{{B_i} + 2{B_{i + \frac{1}{2}}} + {B_{i + 1}}}}{4}}.{\nonumber}$

如果 ${\widehat E_{i + \frac{1}{2}}} < \min \left( {{E_i}, {E_{i + 1}}} \right)$ 或者 ${\widehat E_{i + \frac{1}{2}}} > \max \left( {{E_i}, {E_{i + 1}}} \right)$ , 就通过求解以下方程来替代之前重构得到的 ${h_{i + \frac{1}{2}}}$ :

$\begin{equation} \frac{{{q_{i + \frac{1}{2}}}^2}}{{2{h^2}}} + g\left( {h + {{\overline B }_{i + \frac{1}{2}}}} \right) - {E_{i + \frac{1}{2}}} = 0. \end{equation}$

(2.5)

由此得到 $\overline {{U}} _{i + \frac{1}{2}}^n$ .

(3) 向后投影

类似于向前投影, 首先计算每个交错单元 ${C_{i + \frac{1}{2}}}$ 的能量 $E$ :

${E_{i + \frac{1}{2}}} = \frac{{{u_{i + \frac{1}{2}}}^2}}{2} + g\left( {{h_{i + \frac{1}{2}}} + {{\overline B }_{i + \frac{1}{2}}}} \right).{\nonumber}$

再对 ${{V}}$ 做分片线性重构, 得到 ${h_i}$ , ${q_i}$ 及 ${E_i}$ , 并定义 ${\widehat E_i}$ :

$\begin{equation} {\widehat E_i}: = \frac{{{u_i}^2}}{2} + g\left( {{h_i} + {B_i}} \right). \end{equation}$

(2.6)

如果 ${\widehat E_i} < \min \left( {{E_{i - \frac{1}{2}}}, {E_{i + \frac{1}{2}}}} \right)$ 或者 ${\widehat E_i} > \max \left( {{E_{i - \frac{1}{2}}}, {E_{i + \frac{1}{2}}}} \right)$ , 就通过求解以下方程来替代之前重构的到的 ${h_i}$ :

$\begin{equation} \frac{{{q_i}^2}}{{2{h^2}}} + g\left( {h + {B_i}} \right) - {E_i} = 0. \end{equation}$

(2.7)

(4) 非线性方程求解

考虑上述过程所牵涉到的非线性代数方程

$\varphi \left( h \right): = \frac{{{q^2}}}{{2{h^2}}} + g\left( {h + B} \right) - E = 0,$

其解不唯一, 解法如下:

(Ⅰ) 如果 $q = 0$ , 方程退化为线性方程, $h = \frac{E}{g} - B$ .

(Ⅱ) 如果 $q \ne 0$ , 注意到 $\varphi \left( h \right)$ 为凸函数且在 ${h_0} = \sqrt[3]{{\frac{{{q^2}}}{g}}}$ 处达到极小值. 如果 $\varphi \left( {{h_0}} \right) > 0$ , 方程无正根, 此时无需改变重构值(此情形在数值算例中很少出现); 如果 $\varphi \left( {{h_0}} \right) = 0$ , 方程只有一个正根: $h = {h_0}$ ; 如果 $\varphi \left( {{h_0}} \right) < 0$ , 通过公式法得到方程的三个根

$h = 2\sqrt[3]{{\widehat P}}\cos \left( {\widehat Q + \frac{{2k\pi }}{3}} \right) - \frac{1}{3}\left( {B - \frac{E}{g}} \right), \; \; k = 0, 1, 2.{\nonumber}$

其中

$\widehat P = \sqrt { - {{\left( {\frac{P}{3}} \right)}^3}} , \widehat Q = \frac{1}{3}\arccos \left( {\frac{{ - Q}}{{2\widehat P}}} \right),$

$P = - \frac{1}{3}{\left( {B - \frac{E}{g}} \right)^2}, Q = \frac{2}{{27}}{\left( {B - \frac{E}{g}} \right)^3} + \frac{{{q^2}}}{{2g}}.$

这里面有一个负根和两个正根, 舍弃负根, 通过以下方法来选取符合物理现象的根:

引入弗劳德数 $Fr: = \frac{{\left| u \right|}}{{\sqrt {gh} }} = \frac{{\left| q \right|}}{{\sqrt {g{h^3}} }}$ , 显然当 $h = {h_0} = \sqrt[3]{{\frac{{{q^2}}}{g}}}$ 时, $Fr = 1$ . 而当 $h < {h_0}\left( {h > {h_0}} \right)$ 时, $Fr > 1\left( {Fr < 1} \right)$ , 称此时为超临界流(次临界流). 通过重构值来判断所求单元的状态, 如有 $Fr > 1\left( {Fr < 1} \right)$ , 就选取较小(较大)的解, 使得流体能保持状态不变.

2.3 源项的"和谐性"离散

(1) 源项离散

首先考虑(2.3)式, 注意到

${\left( {h{u^2} + \frac{1}{2}g{h^2}} \right)_x} = h{e_x} + u{q_x},$

其中, $e = \frac{{{u^2}}}{2} + gh$ , 故可将 ${\left( {{{{F}}_i}} \right)^\prime }$ 定义为

$\begin{equation} {\left( {{{{F}}_i}} \right)^\prime }: = \left( \begin{array}{cc} {\left( {{q_i}} \right)^\prime } \\ {h_i}{\left( {{e_i}} \right)^\prime } + {u_i}{\left( {{q_i}} \right)^\prime } \\ \end{array} \right), \end{equation}$

(2.8)

其中 ${\left( {{e_i}} \right)^\prime }$ 也可通过minmod函数来计算, 再借助sensor函数^[4]使 ${{S}}_i^n$ 达到源项于 ${t^n}$ 时刻在非交错单元的二阶近似:

$\begin{equation} {{{S}}_i}: = {{{S}}_{i, L}} + {{{S}}_{i, R}} + {{{S}}_{i, C}}. \end{equation}$

(2.9)

其中

$\left\{ \begin{array}{ll} { } {{\mbox{S}}_{i, L}} = \frac{{{s_i}^2\left( {1 - {s_i}} \right)\left( {2 - {s_i}} \right)}}{6}\left( \begin{array}{cc} 0 \\ { } - g{h_i}\theta \frac{{{B_i} - {B_{i - 1}}}}{{\Delta x}} \\ \end{array} \right), \hfill \\ { } {{\mbox{S}}_{i, R}} = \frac{{{s_i}^2\left( {1 + {s_i}} \right)\left( {2 - {s_i}} \right)}}{2}\left( \begin{array}{cc} 0 \\ { } - g{h_i}\theta \frac{{{B_{i + 1}} - {B_i}}}{{\Delta x}} \\ \end{array} \right), \hfill \\ { } {{\mbox{S}}_{i, C}} = \frac{{{s_i}\left( {{s_i} + 1} \right)\left( {{s_i} - 1} \right)}}{6}\left( \begin{array}{cc} 0 \\ { } - g{h_i}\frac{{{B_{i + 1}} - {B_{i - 1}}}}{{2\Delta x}} \\ \end{array} \right), \hfill \\ \end{array} \right.\quad {s_i} = \left\{ \begin{array}{ll} - 1, &{ } {\left( {{e_i}} \right)^\prime } = \theta \frac{{{e_i} - {e_{i - 1}}}}{{\Delta x}}, \\ 1, &{ } {\left( {{e_i}} \right)^\prime } = \theta \frac{{{e_{i + 1}} - {e_i}}}{{\Delta x}}, \\ 2, &{ } {\left( {{e_i}} \right)^\prime } = \frac{{{e_{i + 1}} - {e_{i - 1}}}}{{2\Delta x}}, \\ 0, &{ } {\left( {{e_i}} \right)^\prime } = 0. \\ \end{array} \right.{\nonumber}$

再考虑(2.2)式关于源项积分的逼近:

$\int_{{x_i}}^{{x_{i + 1}}} {\int_{{t^n}}^{{t^{n + 1}}} {{{S}}\left( {{U}} \right)} } {\rm d}t{\rm d}x \approx \Delta x\Delta t{{S}}\left( {\overline {{U}} _i^{n + \frac{1}{2}}, \overline {{U}} _{i + 1}^{n + \frac{1}{2}}} \right).{\nonumber}$

其中

$\begin{equation} {{S}}\left( {{{\overline {{U}} }_i}, {{\overline {{U}} }_{i + 1}}} \right) = \left( \begin{array}{cc} 0 \\ { } - g \cdot \frac{{{h_i} + {h_{i + 1}}}}{2} \cdot \frac{{{B_{i + 1}} - {B_i}}}{{\Delta x}} + \frac{{{h_{i + 1}} - {h_i}}}{{4\Delta x}} \cdot {\left( {{u_{i + 1}} - {u_i}} \right)^2} \\ \end{array} \right). \end{equation}$

(2.10)

注意到(2.10)式第一项用的是中心差分及中点积分公式, 且对于光滑解, 第二项有

${\left( {{u_{i + 1}} - {u_i}} \right)^2} = O\left( {{{\left( {\Delta x} \right)}^2}} \right),$

故此数值积分也有二阶精度.

(2) "和谐性"的证明

定理2.1 假设在 ${t^n}$ 时刻, 数值解 $\overline {{U}} _i^n$ 于离散水平满足移动稳态条件约束, 即

$\begin{equation} {q_i} \equiv \widehat q, {E_i} = \widehat E, \forall i. \end{equation}$

(2.11)

则(Ⅰ)由(2.3)式得到的预测值 $\overline {{U}} _i^{n + \frac{1}{2}} = \overline {{U}} _i^n$ ; (Ⅱ) 由(2.2)式得到的更新值 $\overline {{U}} _{i + \frac{1}{2}}^{n + 1} = \overline {{U}} _{i + \frac{1}{2}}^n$ .

证首先来证明(Ⅰ), 不妨设 ${\left( {{e_i}} \right)^\prime } = \theta \frac{{{e_{i + 1}} - {e_i}}}{{\Delta x}}$ , 由条件(2.11), ${\left( {{q_i}} \right)^\prime } = 0$ , ${\left( {{E_i}} \right)^\prime } = 0$ , 代入(2.8)和(2.9)式

$\begin{equation} {\left( {{{{F}}_i}} \right)^\prime } = \left( \begin{array}{cc} 0 \\ { } {h_i}\theta \frac{{{e_{i + 1}} - {e_i}}}{{\Delta x}} \\ \end{array} \right), {{{S}}_i} = {{{S}}_{i, R}} = \left( \begin{array}{cc} 0 \\ { } - g{h_i}\theta \frac{{{B_{i + 1}} - {B_i}}}{{\Delta x}} \\ \end{array} \right). \end{equation}$

(2.12)

将(2.12)式代入(2.3)式

$\overline {{U}} _i^{n + \frac{1}{2}} = \overline {{U}} _i^n + \frac{{\Delta t}}{{2\Delta x}}\left( \begin{array}{cc} 0 \\ { } - {h_i}\theta \left( {{E_{i + 1}} - {E_i}} \right) \\ \end{array} \right) = \overline {{U}} _i^n.{\nonumber}$

当 ${\left( {{e_i}} \right)^\prime }$ 取其他值时同理可证.

接下来证明(Ⅱ), 由(Ⅰ), (2.2)式变成

$\overline {{U}} _{i + \frac{1}{2}}^{n + 1} = \overline {{U}} _{i + \frac{1}{2}}^n - \frac{{\Delta t}}{{\Delta x}}\left( {{{F}}\left( {\overline {{U}} _{i + 1}^n} \right) - {{F}}\left( {\overline {{U}} _i^n} \right)} \right) + \Delta t{{S}}\left( {\overline {{U}} _i^n, \overline {{U}} _{i + 1}^n} \right),$

其中

$\begin{equation} {{F}}\left( {{{\overline {{U}} }_{i + 1}}} \right) - {{F}}\left( {{{\overline {{U}} }_i}} \right) = \left( \begin{array}{cc} 0 \\ { } {\widehat q^2}\left( {\frac{1}{{{h_{i + 1}}}} - \frac{1}{{{h_i}}}} \right) + \frac{1}{2}g\left( {{h_{i + 1}} - {h_i}} \right)\left( {{h_{i + 1}} + {h_i}} \right) \\ \end{array} \right). \end{equation}$

(2.13)

由于 ${E_{i + 1}} = {E_i}$ , 移项即有

$\begin{equation} g\left( {{h_{i + 1}} - {h_i}} \right) = \frac{{{{\widehat q}^2}}}{2}\left( {\frac{1}{{{h_i}^2}} - \frac{1}{{{h_{i + 1}}^2}}} \right) - g\left( {{B_{i + 1}} - {B_i}} \right). \end{equation}$

(2.14)

将(2.14)式代入(2.13)式, 有

${{F}}\left( {{{\overline {{U}} }_{i + 1}}} \right) - {{F}}\left( {{{\overline {{U}} }_i}} \right) = \left( \begin{array}{cc} 0 \\ { } - g \cdot \frac{{{h_i} + {h_{i + 1}}}}{2} \cdot \left( {{B_{i + 1}} - {B_i}} \right) + \frac{{{h_{i + 1}} - {h_i}}}{4} \cdot {\left( {\frac{{\widehat q}}{{{h_{i + 1}}}} - \frac{{\widehat q}}{{{h_i}}}} \right)^2} \\ \end{array} \right).{\nonumber}$

由(2.10)式, 有

$- \frac{{\Delta t}}{{\Delta x}}\left( {{{F}}\left( {\overline {{U}} _{i + 1}^n} \right) - {{F}}\left( {\overline {{U}} _i^n} \right)} \right) + \Delta t{{S}}\left( {\overline {{U}} _i^n, \overline {{U}} _{i + 1}^n} \right) = \left( \begin{array}{ll} 0 \hfill \\ 0 \hfill \\ \end{array} \right).{\nonumber}$

故 $\overline {{U}} _{i + \frac{1}{2}}^{n + 1} = \overline {{U}} _{i + \frac{1}{2}}^n$ , 结论得证.

定理2.2 在定理2.1的条件下, 利用向前投影(2.5)以及向后投影(2.7), 本文提出的非交错中心格式是具有“和谐性”的, 即

$\overline {{U}} _i^{n + 1} = \overline {{U}} _i^n.{\nonumber}$

证不妨设相邻三个非交错单元 ${C_{i - 1}}, {C_i}, {C_{i + 1}}$ 在 ${t^n}$ 时刻均为超临界流状态, 即 $Fr > 1$ , 则在向前投影时, 有

$q_{i - \frac{1}{2}}^n = q_{i + \frac{1}{2}}^n = \widehat q, E_{i - \frac{1}{2}}^n = E_{i + \frac{1}{2}}^n = \widehat E.{\nonumber}$

又由minmod参数 $\theta \in [1, 2]$ , 通过重构得到的水深

$h_{i \pm \frac{1}{2}}^n \in \left[ {\min \left( {h_{i - 1}^n, h_i^n, h_{i + 1}^n} \right), \max \left( {h_{i - 1}^n, h_i^n, h_{i + 1}^n} \right)} \right].{\nonumber}$

也就是说对于相邻两个交错单元 ${C_{i - \frac{1}{2}}}, {C_{i + \frac{1}{2}}}$ , 也有 $Fr > 1$ . 由定理2.1, 时间层的更新并不改变守恒变量和能量. 同理在向后投影时, 有

$q_i^{n + 1} = \widehat q, E_i^{n + 1} = \widehat E.{\nonumber}$

且对于非交错单元 ${C_i}$ , 也有 $Fr > 1$ , 因此 $h_i^{n + 1} = h_i^n$ . 次临界流的情形同理可证.

注2.1 为了避免在(临近)干湿界面处水深趋近于0而导致速度变得异常的大, 格式采取文献[6]中的去奇异性策略:

$u = \left\{ \begin{array}{ll} { } \frac{q}{h}, {\quad} &h > \varepsilon , \\ 0, &\mbox{其他}. \\ \end{array} \right.{\nonumber}$

为保持一致性, 重新计算 $q = h \cdot u$ . 在数值算例中, 通常取 $\varepsilon = {10^{ - 8}}$ .

3 数值算例

本节主要展示所提出的维持移动稳态的非交错中心格式(新格式)在一些经典算例上的数值结果, 并将其与文献[5]中的格式(旧格式)进行比较, 以展现新格式在求解移动稳态问题时的优势. 取重力加速度 $g = 9.812$ , minmod参数 $\theta = 1.3$ , CFL数 $\beta = 0.5$ .

首先检验格式在数值上是否二阶收敛, 取底部函数为

$B\left( x \right) = {\sin ^2}\left( {\pi x} \right),$

初值条件为

$h\left( {x, 0} \right) = 5 - B\left( x \right), q\left( {x, 0} \right) = 0.{\nonumber}$

计算区域为 $[0, 1]$ , 输出时间 $t = 10$ . 该算例描述的是光滑底部函数下的静水问题^[4], 即初始值就是精确解. 表 1给出了水深的 $L^1$ 误差及收敛阶, 说明格式有二阶精度.

表 1 水深的 $L^1$ 误差及收敛阶

网格数	$L^1$ 误差	收敛阶
50	$6.26*10^{-4}$	-
100	$1.57*10^{-4}$	2.00
200	$3.90*10^{-5}$	2.01
400	$9.73*10^{-6}$	2.00

新窗口打开| 下载CSV

3.1 维持移动稳态算例

本例验证格式能否维持移动稳态解. 取底部函数为

$B\left( x \right) = \left\{ \begin{array}{ll} 0.2 - 0.05{\left( {x - 10} \right)^2}, &8 \leqslant x \leqslant 12, \\ 0, &\mbox{其他}. \\ \end{array} \right.{\nonumber}$

初值条件分为以下三种情形^[6]:

(a) 超临界流:

$q\left( {x, 0} \right) \equiv 24, E\left( {x, 0} \right) \equiv \frac{{{{24}^2}}}{{2 \times {2^2}}} + 9.812 \times 2.{\nonumber}$

(b) 次临界流:

$q\left( {x, 0} \right) \equiv 4.42, E\left( {x, 0} \right) \equiv \frac{{{{4.42}^2}}}{{2 \times {2^2}}} + 9.812 \times 2.{\nonumber}$

$q\left( {x, 0} \right) \equiv 1.53, E\left( {x, 0} \right) \equiv \frac{3}{2}{\left( {9.812 \times 1.53} \right)^{\frac{2}{3}}} + 9.812 \times 0.2.{\nonumber}$

计算区域为 $[0, 25]$ , 输出时间 $t = 20$ , 网格数200. 通过计算发现, 以上三种情形在机器精度下均可维持稳态解. 这也就从数值上验证了格式是具有“和谐性”的.

注3.1 本例的初值条件给的是 $\left\{ {{q_i}\left( 0 \right)} \right\}$ 和 $\left\{ {{E_i}\left( 0 \right)} \right\}$ , 需通过(2.7)式求解 $\left\{ {{h_i}\left( 0 \right)} \right\}$ . 对于情形(a), 在整个计算区域均有 $Fr > 1$ ; 对于情形(b)有 $Fr < 1$ ; 而对于情形(c), 有

$\left\{ \begin{array}{ll} F{r_i}\left( 0 \right) < 1, {x_i} < 10, \hfill \\ F{r_i}\left( 0 \right) > 1, {x_i} > 10. \hfill \\ \end{array} \right.{\nonumber}$

3.2 捕捉移动稳态小扰动算例

本例测试格式捕捉小扰动的能力. 在例3.1的条件下, 当 $x \in [5.75, 6.25]$ 时, 给处于稳态时的水深 $\left\{ {{h_i}\left( 0 \right)} \right\}$ 加一个小扰动^[6] $\eta = 0.05$ , 网格数200, 情形(a) 的输出时间 $t = 1$ , 情形(b)和情形(c)的输出时间 $t = 1.5$ , 图 1–3给出了新格式与旧格式通过计算得到的水深 $\left\{ {{h_i}\left( t \right)} \right\}$ 和处于稳态时的水深 $\left\{ {{h_i}\left( 0 \right)} \right\}$ 之间的差值. 从图中可以看出, 旧格式在凸起地形处会产生较为明显的虚假震荡, 震荡幅度与截断误差成正比. 而新格式实现了几乎无震荡.

图 1

新窗口打开| 下载原图ZIP| 生成PPT

图 1 情形(a): 旧格式(左图)与新格式(右图)在 $t = 1$ 时的水深和处于移动稳态时的水深之间的差值 $\Delta h$

图 2

新窗口打开| 下载原图ZIP| 生成PPT

图 2 情形(b): 旧格式(左图)与新格式(右图)在 $t = 1.5$ 时的水深和处于移动稳态时的水深之间的差值 $\Delta h$

图 3

新窗口打开| 下载原图ZIP| 生成PPT

图 3 情形(c): 旧格式(左图)与新格式(右图)在 $t = 1.5$ 时的水深和处于移动稳态时的水深之间的差值 $\Delta h$

注3.2 本例忽略了(2.4)和(2.6)式中计算 $\widehat E$ 的过程, 而直接求解(2.5)和(2.7)式来替代通过重构得到的水深.

3.3 收敛到移动稳态算例

本例测试格式能否收敛到例3.1中的移动稳态. 初边值条件如下^[6]:

(a) 超临界流:

$h\left( {x, 0} \right) = 2 - B\left( x \right), q\left( {x, 0} \right) \equiv 0, h\left( {0, t} \right) = 2, q\left( {0, t} \right) = 24.$

(b) 次临界流:

$h\left( {x, 0} \right) = 2 - B\left( x \right), q\left( {x, 0} \right) \equiv 0, q\left( {0, t} \right) = 4.42, h\left( {25, t} \right) = 2.$

$h\left( {x, 0} \right) = 0.66 - B\left( x \right), q\left( {x, 0} \right) \equiv 0, q\left( {0, t} \right) = 1.53, h\left( {25, t} \right) = 0.66.$

输出时间 $t = 200$ , 网格数200, 图 4–6给出了新格式与旧格式的水位, 流量以及能量的比较结果. 从图中可以看出, 相较于旧格式, 新格式能够更精确地捕捉流量以及能量, 且有较好的稳健性.

图 4

新窗口打开| 下载原图ZIP| 生成PPT

图 4 情形(a): 水位(左图), 流量(中图), 能量(右图)

图 5

新窗口打开| 下载原图ZIP| 生成PPT

图 5 情形(b): 水位(左图), 流量(中图), 能量(右图)

图 6

新窗口打开| 下载原图ZIP| 生成PPT

图 6 情形(c): 水位(左图), 流量(中图), 能量(右图)

注3.3 情形(c)的下游边值条件( $h\left( {25, t} \right) = 0.66$ )仅在下游为次临界流时施加.

3.4 溃坝算例

本例用格式模拟实际生活中经常遇到的溃坝问题. 第一个算例模拟的是在含矩形凸起的地形上溃坝水流的流动^[9], 取底部函数为

$B\left( x \right) = \left\{ \begin{array}{rl} 8, &\left| {x - 750} \right| \leqslant 1500/8, \\ 0, &\mbox{其他}, \end{array} \right.{\nonumber}$

初值条件为

$h\left( {x, 0} \right) = \left\{ \begin{array}{rl} 20 - B\left( x \right), &x \leqslant 750, \\ 15 - B\left( x \right), &\mbox{其他}, \end{array} \right.\quad q\left( {x, 0} \right) = 0.{\nonumber}$

计算区域为 $[0, 1500]$ , 输出时间 $t = 15$ 以及 $t = 60$ , 网格数400. 图 7–8给出了水位图和速度图. $t = 15$ 时, 从图中可观察到构成齐次浅水方程黎曼问题解的标准稀疏波和激波, 之后波会到达间断底部的边缘, 一部分波会继续向区域外移动, 一部分波被反射成另一种波, 向相反的方向移动, 剩余部分形成静态的接触间断. 当 $t = 60$ 时, 可观察到数值解有6个波, 其中(1)和(3)是稀疏波, (5)和(6)是激波, (2)和(4)是接触间断.

图 7

新窗口打开| 下载原图ZIP| 生成PPT

图 7 $t = 15$ , 水位(左图), 速度(右图)

图 8

新窗口打开| 下载原图ZIP| 生成PPT

图 8 $t = 60$ , 水位(左图), 速度(右图)

第二个算例模拟的是含有干湿界面的溃坝水流的流动^[6], 将例3.3中情形(a)的初值条件改为

$h\left( {x, 0} \right) = \left\{ \begin{array}{ll} 2, x < 5, \\ 0, \mbox{其他}, \\ \end{array} \right.\quad q\left( {x, 0} \right) = \left\{ \begin{array}{ll} 24, x < 5, \\ 0, \mbox{其他}. \\ \end{array} \right.$

边值条件不变, 网格数200. 图 9给出了输出时间 $t = 0.1, 0.15, 0.3, 0.4, 2, 4$ 的数值结果. 从图中可以看出, 水流会平稳地流过凸起地形, 并在 $t = 4$ 时达到与例3.3中情形(a)相同的稳态. 本例也验证了格式具有保正性.

图 9

新窗口打开| 下载原图ZIP| 生成PPT

图 9 水位(左图), 流量(中图), 能量(右图)

4 总结

本文针对一维Saint-Venant方程组, 提出了一种新的非交错中心有限体积格式. 通过对守恒变量以及能量值的分片线性重构, 使得格式具有二阶精度, 且当流体处于移动稳态时, 能够在向前及向后投影时维持流量与能量的不变性, 再求解非线性方程得到基于稳态时的水深; 通过一种新的源项离散方法, 使格式能维持移动稳态解. 格式很好地避免了经典中心格式在求解移动稳态问题时出现较大误差这一问题. 最后通过数值算例, 验证了该格式的和谐性及稳健性.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Saint-Venant

Théorie du mouvement non-permanent des eaux avec application aux crues des rivires et à lintroduction des marées dans leur lit

Comptes Rendus Hebdomadaires Des Séances De Lacadémie Des Sciences, 1871, 73 (99): 148- 154

[本文引用: 1]

[2]

Nessyahu

, Tadmor

Non-oscillatory central differencing for hyperbolic conservation laws

J Comput Phys, 1990, 87 (2): 408- 463

DOI:10.1016/0021-9991(90)90260-8 [本文引用: 3]

[3]

Touma

Central unstaggered finite volume schemes for hyperbolic systems: applications to unsteady shallow water equations

Appl Math Comput, 2009, 213 (1): 47- 59

[本文引用: 1]

[4]

Touma

Well-balanced central schemes for systems of shallow water equations with wet and dry states

Appl Math Model, 2016, 40 (4): 2929- 2945

DOI:10.1016/j.apm.2015.09.073 [本文引用: 3]

[5]

Dong

, Li

D F

An effect non-staggered central scheme based on new hydrostatic reconstruction

Appl Math Comput, 2019, 372

DOI:10.1016/j.amc.2019.124992 [本文引用: 2]

[6]

Cheng

Y Z

, Kurganov

Moving-water equilibria preserving central-upwind schemes for the shallow water equations

Commun Math Sci, 2016, 14 (6): 1643- 1663

DOI:10.4310/CMS.2016.v14.n6.a9 [本文引用: 6]

[7]

Cheng

Y Z

, Chertock

, Herty

, et al.

A new approach for designing moving-water equilibria preserving schemes for the shallow water equations

J Sci Comput, 2019, 80 (1): 538- 554

DOI:10.1007/s10915-019-00947-w [本文引用: 1]

[8]

Liu

, Chen

, Jin

, et al.

Moving-water equilibria preserving partial relaxation scheme for the Saint-Venant system

SIAM J Sci Comput, 2020, 42 (4): A2206- A2229

DOI:10.1137/19M1258098 [本文引用: 1]

[9]

Noelle

, Xing

Y L

, Shu

C W

High-order well-balanced finite volume WENO schemes for shallow water equation with moving water

J Comput Phys, 2007, 226 (1): 29- 58

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

[10]

Xing

Y L

, Shu

C W

, Noelle

On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations

J Sci Comput, 2011, 48 (1): 339- 349

[11]

Xing

Y L

Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium

J Comput Phys, 2014, 257, 536- 553

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

[12]

Kurganov

, Petrova

A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system

Commun Math Sci, 2007, 5 (1): 133- 160

DOI:10.4310/CMS.2007.v5.n1.a6 [本文引用: 2]

[13]

Chen

G X

, Noelle

A new hydrostatic reconstruction scheme based on subcell reconstructions

SIAM J Numer Anal, 2017, 55 (2): 758- 784

DOI:10.1137/15M1053074

[14]

Jiang

G S

, Levy

, Lin

C T

, et al.

High-resolution nonoscillatory central schemes with nonstaggered grids for hyperbolic conservation laws

SIAM J Numer Anal, 1998, 35 (6): 2147- 2168

DOI:10.1137/S0036142997317560

[15]

Touma

Unstaggered central schemes with constrained transport treatment for ideal and shallow water magnetohydrodynamics

Appl Numer Math, 2010, 60 (7): 752- 766

DOI:10.1016/j.apnum.2010.02.006

[16]

Touma

, Khankan

Well-balanced unstaggered central schemes for one and two-dimensional shallow water equation systems

Appl Math Comput, 2012, 218 (10): 5948- 5960

[17]

Bollermann

, Chen

G X

, Kurganov

, et al.

A well-balanced reconstruction for wet/dry fronts for the shallow water equations

J Sci Comput, 2013, 56 (2): 267- 290

DOI:10.1007/s10915-012-9677-5 [本文引用: 1]

Théorie du mouvement non-permanent des eaux avec application aux crues des rivires et à lintroduction des marées dans leur lit

1871

... Saint-Venant方程组由Saint-Venant于1871年提出^[1], 至今仍在浅水流、涌浪传播、水库和明渠流等问题的数值求解中发挥着重要作用. 在一维情形下, 其向量形式可表示为 ...

Non-oscillatory central differencing for hyperbolic conservation laws

1990

... Nessyahu和Tadmor于1990年提出了用于求解双曲守恒律问题的二阶交错中心格式^[2], 相较于迎风格式, 其优势在于避免了求解黎曼问题, 但于相邻时间层, 数值解会在交错单元与非交错单元之间来回更替, 导致了复杂的边界条件设置. 为此, Touma提出了稳健的二阶非交错中心格式^[3], 并将其应用于Saint-Venant方程组, 提出了保正的能维持静稳态解的格式^[4], 但在干湿界面处, “保正性”与“和谐性”无法同时得到满足. Dong和Li^[5]通过重构守恒变量及水位值, 且在干湿界面处构造关于水位的映射, 使格式在干湿界面处经过“保正性”处理后仍可维持静稳态解. 然而在用中心格式求解移动稳态问题时, 在底部地形变化处(

${B_x} \ne 0$

) 流量与能量会产生较为明显的虚假振荡. ...

... 其中

${\left( {{{\overline {{U}} }_i}} \right)^\prime }$

表示对

${{U}}$

在点

$x = {x_i}$

的导数的逼近, 用minmod函数^[2]来计算 ...

... 其中,

$\beta$

为CFL数, 通常取

$\beta \leqslant 1$

. 具体证明可参考文献[2]. ...

Central unstaggered finite volume schemes for hyperbolic systems: applications to unsteady shallow water equations

2009

${B_x} \ne 0$

) 流量与能量会产生较为明显的虚假振荡. ...

Well-balanced central schemes for systems of shallow water equations with wet and dry states

2016

${B_x} \ne 0$

) 流量与能量会产生较为明显的虚假振荡. ...

... 其中

${\left( {{e_i}} \right)^\prime }$

也可通过minmod函数来计算, 再借助sensor函数^[4]使

${{S}}_i^n$

达到源项于

${t^n}$

时刻在非交错单元的二阶近似: ...

... 计算区域为

$[0, 1]$

, 输出时间

$t = 10$

. 该算例描述的是光滑底部函数下的静水问题^[4], 即初始值就是精确解. 表 1给出了水深的

$L^1$

误差及收敛阶, 说明格式有二阶精度. ...

An effect non-staggered central scheme based on new hydrostatic reconstruction

2019

${B_x} \ne 0$

) 流量与能量会产生较为明显的虚假振荡. ...

... 本节主要展示所提出的维持移动稳态的非交错中心格式(新格式)在一些经典算例上的数值结果, 并将其与文献[5]中的格式(旧格式)进行比较, 以展现新格式在求解移动稳态问题时的优势. 取重力加速度

$g = 9.812$

, minmod参数

$\theta = 1.3$

, CFL数

$\beta = 0.5$

. ...

Moving-water equilibria preserving central-upwind schemes for the shallow water equations

2016

... 由于能量与水深之间的非线性关系, 设计能维持移动稳态解的数值格式往往难度更大. 为此, Kurganov等人提出了中心-迎风格式^{[6, 7]}, Liu等人提出了部分松弛格式^[8], Shu等人提出了加权本质非振荡(WENO) 格式及间断有限元(DG)格式^[9-11]. 这些格式很好地避免了中心格式存在的问题. 更多结果参见文献[12-17]. ...

... 注2.1 为了避免在(临近)干湿界面处水深趋近于0而导致速度变得异常的大, 格式采取文献[6]中的去奇异性策略: ...

... 初值条件分为以下三种情形^[6]: ...

... 本例测试格式捕捉小扰动的能力. 在例3.1的条件下, 当

$x \in [5.75, 6.25]$

时, 给处于稳态时的水深

$\left\{ {{h_i}\left( 0 \right)} \right\}$

加一个小扰动^[6]

$\eta = 0.05$

, 网格数200, 情形(a) 的输出时间

$t = 1$

, 情形(b)和情形(c)的输出时间

$t = 1.5$

, 图 1–3给出了新格式与旧格式通过计算得到的水深

$\left\{ {{h_i}\left( t \right)} \right\}$

和处于稳态时的水深

$\left\{ {{h_i}\left( 0 \right)} \right\}$

之间的差值. 从图中可以看出, 旧格式在凸起地形处会产生较为明显的虚假震荡, 震荡幅度与截断误差成正比. 而新格式实现了几乎无震荡. ...

... 本例测试格式能否收敛到例3.1中的移动稳态. 初边值条件如下^[6]: ...

... 第二个算例模拟的是含有干湿界面的溃坝水流的流动^[6], 将例3.3中情形(a)的初值条件改为 ...

A new approach for designing moving-water equilibria preserving schemes for the shallow water equations

2019

Moving-water equilibria preserving partial relaxation scheme for the Saint-Venant system

2020

High-order well-balanced finite volume WENO schemes for shallow water equation with moving water

2007

... 本例用格式模拟实际生活中经常遇到的溃坝问题. 第一个算例模拟的是在含矩形凸起的地形上溃坝水流的流动^[9], 取底部函数为 ...

On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations

2011

Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium

2014

A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system

2007

... 将底部地形函数

$B\left( x \right)$

用分片线性函数

$\widetilde B\left( x \right)$

近似^[12]: ...

A new hydrostatic reconstruction scheme based on subcell reconstructions

2017

High-resolution nonoscillatory central schemes with nonstaggered grids for hyperbolic conservation laws

1998

Unstaggered central schemes with constrained transport treatment for ideal and shallow water magnetohydrodynamics

2010

Well-balanced unstaggered central schemes for one and two-dimensional shallow water equation systems

2012

A well-balanced reconstruction for wet/dry fronts for the shallow water equations

2013

〈

〉