浅水波方程是流体动力学中描述流体运动规律的数学模型, 在海洋、清洁能源、环境与水利工程等众多领域具有广泛应用, 该模型属于非线性双曲守恒律方程, 对应解的结构十分复杂, 很难得到解析解, 因此构造合适的数值格式成为一个研究热点. 目前已有多种数值方法可用来求解浅水波方程, 如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], 从而得到一类新的高阶熵稳定格式. 为了证明格式的保平衡性, 先介绍了二阶熵守恒格式以及具有一阶耗散的低阶熵稳定格式的保平衡性, 再考虑对源项进行数值离散, 离散方式与所采用的偶数阶熵守恒通量的组合方式相对应, 进而证明高阶熵稳定格式的保平衡性定理, 并通过理论推导、数值实验验证新格式的保平衡性.

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

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

其中, 守恒变量${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$常数时, 非零通量与源项保持平衡.

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

其中守恒变量${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$常数时, 各个方向的非零通量与源项保持平衡.

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

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

${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^[10] 与${F}$相容的通量${F}_{i+1/2}$是熵守恒数值通量, 如果等式

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

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

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

那么称格式(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].

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

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

其中系数$\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}$处的左右重构值为

定义高阶熵稳定数值通量

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

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

其中$\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]为

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

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

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

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

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

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

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

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

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

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

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

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

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

其中, $\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 )$.

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

其中, $\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)的有限体积格式为

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

根据恒等式

进一步地, 有

又由于对任意的$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.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}$处的重构跳跃值为

通过熵变量${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)具有保平衡性.证毕.

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

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

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

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

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

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

该算例选自文献[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格式有效地改善了抹平现象.

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

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

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

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

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

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

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

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

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

该算例选自文献[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]的数值结果, 这里用更少的网格数达到了不错的效果, 避免了由于网格数过大而造成的计算困难.

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

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