浅水波由Saint-Venant方程组^[4]来描述,这个方程组大约在130年前就被提出来了,至今仍然广泛的应用于河流、湖泊或是海洋的水的流动.在一维的情形下, Saint-Venant方程组为

初始条件为

这里$h(x, t)$表示水的深度, $u(x, t)$是速度, $g$是重力加速度常数,函数$B(x)$表示底部地形,假设它是和时间无关的并且可能是不连续的.我们在一个确定的空间区域$X$考虑方程组$(1.1)$,如果$X\neq R$那么Saint-Venant方程组需要加上合适的边界条件.

我们主要考虑方程组的稳态解

在许多的实际应用当中,我们考虑其中一类特殊的稳态解称为"Lake at rest",也称为"静稳态解"

当$h=0$时,这样的解我们称为"Dry lake",求解这类问题的数值格式我们可以参考^[5-12].数值格式很大的一个挑战就是解决这样的稳态解以及它们的一个小的扰动问题.如果数值格式能够很好地维持这样的解,我们称这个格式就具有"和谐性",同时称格式是和谐格式.和谐格式的核心问题是当求解稳态解时,数值通量可以精确得与数值源项达到平衡.由于物理和数学因素的限制,要求格式具有"保正性",即$h$非负(在文献[10-11]里可以了解更多"保正性"格式).为了格式书写方便,我们将浅水波方程组写成如下向量形式:

本论文比较研究了:第2节, KP格式: KP格式^[1]主要是基于半离散的中心迎风格式,在干湿单元进行重构以使得$h\geq0$,但是当格式具有了保正性后,没有了和谐性;第3节, BCKN格式: BCKN格式^[2]解决了KP格式的不足,弥补了KP格式不具有和谐性的缺陷,核心就是建立了$\omega_{j}$和$\overline{h}_{j}$之间的映射;第4节, T格式: T格式^[3]是建立在交错网格上的中心格式,格式在干湿单元具有保正性,但是不足之处依然是和谐性被破坏;第5节,展示格式的算例.

KP格式也称为半离散的中心迎风格式

这里$H_{j+\frac{1}{2}}$表示中心迎风格式的数值通量, $\overline{S}_{j}$表示源项的积分平均的合适的离散.

$\overline{U}_{j}(t)$表示近似的单元积分平均值

离散的源项$(2.2)$的第二组成部分(细节可以看文献[11]和[12])

中心迎风数值通量为

这里用以下通量记号

这里的值$U_{j+\frac{1}{2}}^{\pm}=(w_{j+\frac{1}{2}}^{\pm}, h_{j+\frac{1}{2}}^{\pm}\cdot u_{j+\frac{1}{2}}^{\pm})$表示解在点$x_{j+\frac{1}{2}}$左右两边的值

这个解由分片线性重构给出

斜率限制器选择minmod限制器.

最后$(2.4)$式的局部速度$a^{\pm}_{j+\frac{1}{2}}$由雅克比矩阵$\frac{\partial F}{\partial U}$的特征值得到

为了和谐性,在不影响守恒性的同时,人们往往改用对自由面$w$进行重构,从而间接达到对水深$h$进行重构的目的.然而这种方法会破坏$h$的保正性. KP格式的核心就是为了解决这个问题,其具体实现步骤如下

如果$\omega_{j+\frac{1}{2}}^{-}<B_{j+\frac{1}{2}}$,令$(\omega_{x})_{j}:=\frac{B_{j+\frac{1}{2}}-\overline{\omega}_{j}}{\Delta x/2}$那么

如果$\omega_{j-\frac{1}{2}}^{+}<B_{j-\frac{1}{2}}$,令$(\omega_{x})_{j}:=\frac{\overline{\omega}_{j}-B_{j-\frac{1}{2}}}{\Delta x/2}$那么

这种重构方法,虽然保证了在重构过程中水深非负性,但是在干湿界面附近的和谐性出现问题,这从图 2.1中可以体现.

因为如果满足$\omega_{j+\frac{1}{2}}^{-}<B_{j+\frac{1}{2}}$,那么

进而可以得到$\omega_{j+\frac{1}{2}}^{+}=2\overline{\omega}_{j+1}-B_{j+\frac{3}{2}}$, $\omega_{j-\frac{1}{2}}^{-}=B_{j-\frac{1}{2}}$显然从$(2.7)$式中可以看出$\omega_{j+\frac{1}{2}}^{+}$和$\omega_{j+\frac{1}{2}}^{-}$不一定相等以及$\omega_{j-\frac{1}{2}}^{+}$和$\omega_{j-\frac{1}{2}}^{-}$同样不一定相等.如果满足$\omega_{j-\frac{1}{2}}^{+}<B_{j-\frac{1}{2}}$,同样的方法我们可以得到$\omega_{j+\frac{1}{2}}^{+}=B_{j+\frac{1}{2}}$, $\omega_{j-\frac{1}{2}}^{-}=2\overline{\omega}_{j-1}-B_{j-\frac{3}{2}}$.从$(2.8)$式中同样可以看出$\omega_{j+\frac{1}{2}}^{+}$和$\omega_{j+\frac{1}{2}}^{-}$不一定相等以及$\omega_{j-\frac{1}{2}}^{+}$和$\omega_{j-\frac{1}{2}}^{-}$同样不一定相等.所以KP格式的不足就是重构的过程中破坏了和谐性.

BCKN格式可以看作是KP格式的一个改进.在研究了KP格式之后,我们可以发现KP格式的不足之处是在干湿交错单元没有和谐性, BCKN格式解决了这个问题.接下来,我们研究BCKN格式的重构过程.

当出现干湿界面的情况时,中心迎风格式在重构阶段有可能使得格式产生负的水深, 图 3.1可以帮我们很好的理解这一点.为了叙述简单只考虑以下情形: $w_{j}>B_{j}$,而且在分片线性重构中的斜率$(w_{x})_{j}$和$(u_{x})_{j}$已经使用了非线性的限制器.

假设在某个几乎干的单元

而且在单元$j+1$, $w$的重构的值满足

也就是单元$j+1$是全淹没的.这意味着单元$j$位于干的边界的附近(海滨).在单元$j$上设计一个和谐的重构.

计算单元$j$的自由面(用$w_{j}$表示),在这个单元假设水面是平静的.

如果单元$j$是全淹没单元,函数$w_{j}(x)$的定义

否则是连续的分片线性函数

这样定义以后,在单元$j$的积分平均值为

其中$x^{\star}_{w}$表示干湿分界点.进而可以导出自由面$w_{j}$在干/湿单元的计算公式

至此可以得到$\omega_{j}$和$\overline{h}_{j}$之间的映射,这是BCKN格式的核心之处.

关于BCKN格式,在半淹没单元$j$修正$h$的重构以确保格式具有和谐性.首先设置$w^{-}_{j+\frac{1}{2}}=w^{+}_{j+\frac{1}{2}}$ (这样的设置可以保证在提高一边水的高度时另一边的和谐性不会遭到破坏), $w_{j}(x)$在单元$j$的重构由$\overline{h}_{j}$的守恒性决定.如果单元$j$的水足够多(图 3.2左边所示),那么有唯一解$h^{+}_{j-\frac{1}{2}}\geq0$满足

可得$w^{+}_{j-\frac{1}{2}}=h^{+}_{j-\frac{1}{2}}+B_{j-\frac{1}{2}}$,这样在单元$j$的重构就完成了.

根据守恒性要求,如果$(3.5)$式求得的$h^{+}_{j-\frac{1}{2}}$的值是负的,那么我们就在单元$j$对$h(x, t)$进行分段线性重构.$\!\!$"湿"和"干"的分界点由$x^{\ast}_{j}$来表示

这里

综合两种情形我们可以得到重构的值

从重构的整个过程不难理解, BCKN格式的核心思想就是建立了$\omega_{j}$与$\overline{h}_{j}$之间的映射,使得格式具有和谐性.

T格式虽然也是中心格式,但是与KP格式和BCKN格式都有着很大的区别. KP格式和BCKN格式使用的网格都是非交错网格,而T格式的网格则是交错网格,因为方程的解同样是分片线性的一个近似解,格式同样具有二阶精度. T格式的主要思想大致可以分为三个部分: "向前投影", "预测步", "向后投影".将方程$(1.2)$在矩形$R_{j}^{n}=[x_{j}, x_{j+1}]\times[t^{n}, t^{n+1}]$上进行积分并利用高斯定理可得

其中

通量满足

因为对方程的解进行了线性化处理,即$U(x, t^{n})\approx \iota_{j}(x, t^{n})=U_{j}^{n}+(x-x_{j})U'_{j}$,其中的$U'_{j}$表示有限的数值导数是对$U$在点$x_{j}$的导数的逼近,限制器选择的是minmod函数.

通过求解积分

从而可以得到"向前投影"

同理可以得到"向后投影"

通过对上式积分等式的处理我们可以得到交错网格上的中心格式

其中

源项的离散为

求解这个步骤还需要"预测步"

其中$S_{j}^{n}$表示底部函数的离散,这里就要用到sensor函数.那么为了解决"向前投影"可能产生负的水深这个问题.

T格式的修正过程为

为了叙述的简洁,省略时间记号并且我们只描述两种情形: (1) $h_{j}\geq 0, h_{j+1}=0$; (2) $h_{j}=0, h_{j+1}\geq 0$,首先定义一个足够小的数$\epsilon=10^{-9}$.

第一种情形

如果$h_{j}>\epsilon$,那么令$h'_{j}=-\frac{h_{j}}{\Delta x/4}$,同时为了保持格式具有和谐性底部函数的斜率为$z'_{j}=\frac{z_{j+\frac{1}{2}}-z_{j}}{\Delta x/2}.$

如果$h_{j}<\epsilon$,那么令$h_{j}=0, h'_{j}=0$.

第二种情形

如果$h_{j+1}>\epsilon$,那么令$h'_{j+1}=\frac{h_{j+1}}{\Delta x/4}$,同时为了保持格式具有和谐性底部函数的斜率为$z'_{j+1}=\frac{z_{j+1}-z_{j+\frac{1}{2}}}{\Delta x/2}.$

如果$h_{j+1}<\epsilon$,那么令$h_{j+1}=0, h'_{j+1}=0$.

这样的修正的过程可以使得格式具有保正性,但是却破坏了和谐性.考虑如图 4.1中所示有下面定理.

定理4.1  考虑计算区域出现干湿界面:取$\varepsilon=10^{-9}$,如果当$h_{j}>\varepsilon, h_{j+1}<\varepsilon$时, $h'_{j}=-\frac{h_{j}}{\Delta x/4}, h'_{j+1}=0, B'_{j}= \frac{B_{j+\frac{1}{2}}-B_{j}}{\Delta x/2}$,那么$(hu)_{j}^{n+1/2}=(hu)_{j}^{n}$不一定成立.

证  不妨假设$B_{j+1}>B_{j}$,由所给条件易得$\overline{h}_{j+\frac{1}{2}}=0$以及$B_{j+\frac{1}{2}}=2\overline{w}_{j+\frac{1}{2}}-\frac{1}{2}(B_{j+1}+B_{j})$.

结论得证.

T格式在交错单元解的更新公式中关于源项的离散不能保证守恒性,如图 4.2中所示.水深的积分平均值

事实上水深的积分平均值为

基于以上两种原因易知格式不满足和谐性要求.

我们在这一部分主要研究: KP格式, BCKN格式, T格式的部分应用.这里我们取KP格式和BCKN格式的斜率限制器为$\theta =0.36$以及$CFL$条件: $CFL=0.5$; T格式的斜率限制器为$\theta=1.5$以及$CFL$条件: $CFL=0.46$,所有格式的重力加速度$g=9.8$.

我们首先考虑静水情形,即在有水的部分$W:=h+B=const.$取底部函数为

水深的初始值以及水的初始速度为

我们取三种格式的统一输出时间为: $t=0.1$,计算的区域长度为$L=1$,单元个数$cells=800$. 图 5.1中展示的是水平面$W:=h+B$的数值解. 图 5.2-5.4展示的是三种格式计算的$h*u$的数值解.我们从图 5.2-5.4中可以明显的看到, BCKN格式具有和谐性,可以很好的逼近稳态解$h*u=0$. KP格式与T格式都没有和谐性,但是KP格式要比T格式计算的更精确一些.

其次我们考虑当水面有一定扰动的情形.同样取底部函数为

水深的初始值以及水的初始速度为

图 5.5中展示的是三种中心格式在输出时间均为$t=0.1$,区间长度$[0, 1]$,单元个数为$cells=800$的时候$h+B$的计算结果. 图 5.6计算的是$h*u$的数值解.从两幅图中可以看到,当水面有一定扰动的时候, T格式计算的误差较大,这主要是因为T格式不满足守恒性, BCKN格式和KP格式计算的精度较高一些.

在这一部分我们主要考虑以下黎曼问题.底部函数为$B(x)=0$,水深的初始值为

水的初始速度$u(x, 0)=0$从所给的初始条件我们可以知道,在解的结构中同时会有$\mbox{"稀疏波"}$和"激波"的产生.对于高精度的数值格式往往会在激波附近产生振荡,产生振荡的原因往往是格式不具有TVD性质,因此需要合适的限制器来解决振荡问题.我们用三种格式来解决此情况的黎曼问题时,选择单元个数$cells=800$,输出时间为$t=0.1$,区域长度$[0, 1]$.这里需要指出的是KP格式和BCKN格式我们选择的斜率限制器为$\theta=0.36$, T格式的斜率限制器$\theta=1.5$. 图 5.7和图 5.8是计算的结果. 图 5.7的展示的是稀疏波和激波的数值解与精确解的比较,从图中可以看出T格式的拟合程度要比KP格式和BCKN格式都好,即T格式的数值解要更精确一些,注意到这这个算例当中底部是平的,因此T格式在这个情形中是满足守恒性的,而且计算精度要更高. 图 5.8展示的是$h*u$的计算结果,从图中看以看出向左传播的激波与向两边传播的稀疏波.

其次我们考虑以下黎曼问题.底部函数为$B(x)=0$,水深的初始值为

水的初始速度$u(x, 0)=0.$.我们从所给的初始条件容易判断出,方程的解只有"稀疏波".这里我们取三种格式的计算单元个数均为$cells=400$,输出时间$t=0.1$,计算的区域长度为$[0, 1]$. 图 5.9展示的是三种格式的数值解与精确解的比较, 图 5.9计算的是水面$h+B$的数值解与精确解的比较, 图 5.10是$h*u$的数值解与精确解的比较.从图中可以看出T格式的数值解要比KP格式和BCKN格式的数值解更精确一些, BCKN格式的数值解要比KP格式的数值解更精确一些.

论文总结研究了受关注的浅水波方程的离散格式: KP格式, BCKN格式, T格式. KP格式和BCKN格式中都需要局部速度的值,因此又称为中心-迎风格式. T格式可以看作是NT格式在带有源项的浅水波方程上的应用.我们在论文中主要是分析研究三种格式在求解带有源项的浅水波方程上时的优势以及不足之处.浅水波方程对于其数值格式有较高的要求,在实际应用中我们更关心在稳态解附近的行为,所以要求格式具有和谐性,特别是当计算区域出现干湿界面的时候需要保持水深恒为非负,同时又要求数值格式具有较高的精度.通过一些算例的展示我们可以看到BCKN格式可以满足所有的要求.这里需要特别指出的是我们在应用BCKN格式和KP格式时斜率限制器的选择为$\theta=0.36$显然这里还有问题需要解决.我们今后的工作将关注于此.