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

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

其中

(1.2) $ \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} $

$ 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) $已知, 若求解区域并非整个实数域, 还需添加适当的边值条件.

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

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

其中$ 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节, 对格式做出总结.

针对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}}] $上积分并利用高斯定理得

(2.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} $

在$ {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)式利用中点积分公式求通量的时间积分, 得到其离散形式

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

(2.3) $ \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} $

$ {\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} $

(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}}} $如下

(2.4) $ \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} $

其中, $ {\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}}} $:

(2.5) $ \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} $

由此得到$ \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} $:

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

如果$ {\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} $:

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

(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) $, 就选取较小(较大)的解, 使得流体能保持状态不变.

(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 } $定义为

(2.8) $ \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} $

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

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

其中

$ \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} $

其中

(2.10) $ \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)式第一项用的是中心差分及中点积分公式, 且对于光滑解, 第二项有

$ {\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 $于离散水平满足移动稳态条件约束, 即

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

则(Ⅰ)由(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)式

(2.12) $ \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.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), $

其中

(2.13) $ \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} $

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

(2.14) $ \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.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}} $.

本节主要展示所提出的维持移动稳态的非交错中心格式(新格式)在一些经典算例上的数值结果, 并将其与文献[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 $误差及收敛阶, 说明格式有二阶精度.

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

$ 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} $

(c) 跨临界流:

$ 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.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.2  本例忽略了(2.4)和(2.6)式中计算$ \widehat E $的过程, 而直接求解(2.5)和(2.7)式来替代通过重构得到的水深.

本例测试格式能否收敛到例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. $

(c) 跨临界流:

$ 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给出了新格式与旧格式的水位, 流量以及能量的比较结果. 从图中可以看出, 相较于旧格式, 新格式能够更精确地捕捉流量以及能量, 且有较好的稳健性.

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

本例用格式模拟实际生活中经常遇到的溃坝问题. 第一个算例模拟的是在含矩形凸起的地形上溃坝水流的流动^[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)是接触间断.

第二个算例模拟的是含有干湿界面的溃坝水流的流动^[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)相同的稳态. 本例也验证了格式具有保正性.

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