基于改进的三阶SPWENO重构方法求解浅水波方程

Solving the Shallow Water Equation Based on the Improved Third-Order SPWENO Reconstruction Method

摘要: 计算了离散情况下跨越间断时的具体熵增量，给出相应的耗散项，从而修正了求解浅水波方程的熵守恒格式。结合Roe的迎风耗散项，构造了一阶精度的熵稳定格式。为了提高格式的精度，利用改进的三阶SPWENO重构方法对熵变量进行了重构，得到了一类高精度高分辨率的熵稳定格式，并通过数值算例验证了新格式的有效性。

Abstract: In discrete situation, the specific entropy production when crossing the discontinuities is calculated, and the corresponding dissipation term is given, thus correcting the entropy conservation scheme for shallow water equation. Adding the Roe-type upwind dissipation term, the first-order entropy stable scheme is constructed. In order to improve the accuracy of the scheme, the entropy variable is reconstructed by the improved third-order SPWENO reconstruction to obtain a kind of high-precision high-resolution entropy stable scheme. Finally, the validity of the new scheme is verified by numerical experiments.

1. 引言

浅水波方程是主要用于描述浅水环境下各类流体活动的一类方程，对于研究潮汐、灌溉和浅水海洋等波动问题有着非常重要的作用。因此，浅水波方程的数值求解方法引起了计算流体力学学者们的重视。

浅水波方程属于非线性双曲守恒律方程，人们对这类方程研究发现：不管给定的初始条件是否光滑，该方程的数值解在某个时刻也会产生间断。这使得数值方法的稳定性结果很难获得。针对这一问题，Lax于1954年首次提出弱解的概念 [1]，允许间断解的存在。然而弱解并不满足唯一性，需要借助额外的限制条件在众多弱解中寻找真正具有物理意义的相关解。从物理角度，双曲守恒律方程可以看成粘性趋于零的极限形式；数学上可以证明：满足“粘性消失”的Cauchy问题的解在任何时刻唯一且具有物理意义，这种解称为熵解。在1973年Lax又证明了熵稳定条件与“粘性消失”机制是等价的 [2]。满足熵稳定条件的弱解就是我们要找的熵解，而不满足熵稳定条件的数值解被称为非物理解，在数值模拟结果中往往表现出如“膨胀激波”、“红斑”等非物理现象 [3]。

自von Neumann和Richtmyer起，人们通过引入额外的数值粘性项来满足熵稳定条件，从而避免非物理现象的产生 [4]。1987年Tadmor引入了熵势和熵变量的概念，定义了一类二阶的熵守恒格式 [5]，记为EC (Entropy Conservation)格式，其数值通量保持总熵不变。该格式能够捕捉光滑区域的解，但是由于总熵并没有发生耗散，在捕捉间断解时产生振荡现象。2006年，Roe在熵守恒格式的基础上通过添加Roe粘性项来满足熵稳定条件，得出一类熵稳定格式 [6]，记为ES (Entropy Stable)格式，该格式在捕捉间断时效果显著。此外，Fjordholm，Mishra及Tadmor在熵守恒通量的基础上引入了一个新的Roe型扩散算子，构造了一种求解浅水波方程的二阶精度能量稳定格式 [7]。Tadmor和Zhong构造出了求解二维浅水波方程的能量稳定格式 [8]。Roe和Ismail在2009年通过进一步分析熵稳定格式的熵耗散的大小，提出了对数值粘性更精确量化的一类熵稳定格式 [9]，并应用该格式对Euler方程进行了数值模拟。该格式进一步控制了激波处的熵增，有效消除了膨胀激波和间断处的伪振荡。

近年来，高精度熵稳定格式的构造方法成为计算流体力学研究者致力突破的一个重要课题。2012年，Lefloch和Rohde提出了通过二阶熵守恒通量的线性组合构造任意偶数阶熵守恒通量的方法 [10]。Fjordholm等人则在文献 [11] 中提出在高阶熵守恒通量的基础上，通过适当的重构获得高阶的数值耗散算子，两者结合得到高阶熵稳定通量。他还强调，重构过程必须满足符号性质，从而保证最后的数值格式是熵稳定的。2015年，Fjordholm和Ray提出了一种满足符号性质的WENO (Weighted Essentially Non-Oscillatory)重构方法，记为SPWENO (Sign Preserving WENO)重构 [12]，并将应用该重构得到的新格式对一维Burgers方程进行了数值模拟。国内的刘友琼、郑素佩等人也通过利用不同的重构方法对浅水波方程的高精度解进行了深入研究 [13] [14]。

2. 数值格式

考虑一维的无粘浅水波方程：

$\frac{\partial }{\partial t}u+\frac{\partial }{\partial x}f\left(u\right)=0,x\in R,t>0,$ (1)

其中：守恒型变量 $u={\left[h,uh\right]}^{\text{T}}$，通量 $f\left(u\right)={\left[hu,\frac{1}{2}g{h}^2+h{u}^2\right]}^{\text{T}}$， $h=h\left(x,t\right)$ 表示水的深度， $u=u\left(x,t\right)$ 是

水的平均速度，g是重力加速度，计算时一般取为1。

本文采用有限体积方法，对空间方向上的求解区域进行均匀网格剖分，网格步长取 $\Delta x$，得到一组网

格点 $x_i=i\Delta x,i=0,1,\cdots ,N$ 和一系列互不相交的单元 $I_i=\left[x_{i-1/2},x_{i+1/2}\right]$， $x_{i\pm 1/2}=x_i\pm \frac{\Delta x}{2}$，则方程(1)的守恒型

半离散格式为：

$\frac{\text{d}}{\text{d}t}{u}_i\left(t\right)=-\frac{1}{\Delta x}\left({\hat{f}}_{i+1/2}-{\hat{f}}_{i-1/2}\right),$ (2)

其中 $\hat{f}$ 是与 $f\left(u\right)$ 相容的数值通量，满足 $\hat{f}\left(u,u,\cdots ,u\right)=f\left(u\right)$。

2.1. 熵守恒格式

2004年，Tadmor等人提出了适用于Navier-Stokes方程、浅水波方程等各种守恒系统的分段线性路径的熵守恒格式，然而该格式在计算每个通量时都包含复数表达式，在计算时容易出现分母为0的警告，而且其数值稳定性还存在一定的争议。于是，2008 年Fjordholm等人根据浅水波方程自身的特点，构造了简单的显式熵守恒EEC (Explicit Energy Conserving)格式，避免了这些缺点，其数值通量如下：

$f_{i+1/2}^{\text{EEC}}=\left[\begin{array}{c}{\bar{h}}_{i+1/2}{\bar{u}}_{i+1/2}\\ {\bar{h}}_{i+1/2}{\left({\bar{u}}_{i+1/2}\right)}^2+\frac{1}{2}g\overline{h^2}{}_{i+1/2}\end{array}\right],$ (3)

这里， ${\overline{(\cdot )}}_{i+1/2}=\frac{{(\cdot )}_i+{(\cdot )}_{i+1}}{2},{\overline{{(\cdot )}^2}}_{i+1/2}=\frac{{(\cdot )}_i^2+{(\cdot )}_{i+1}^2}{2}$。熵守恒格式得到的数值解保持离散总熵不变，熵增为

0。该格式在光滑区域表现良好，但是由于缺少熵的耗散机制，在间断区域表现出强烈的振荡，求出的数值解表现出较强的色散效应。

任意偶数阶熵守恒数值通量可以通过二阶熵守恒通量的线性组合来构造，即

$f_{i+1/2}^{\text{EC},2p}={\displaystyle \underset{r=1}{\overset{p}{\sum }}{\alpha }_r^p}{\displaystyle \underset{s=0}{\overset{r-1}{\sum }}{f}^{\text{EC}}}\left(u_{i-s},\cdots ,u_{i-s+r}\right),$ (4)

式中 ${\alpha }_1^p,{\alpha }_2^p,\cdots ,{\alpha }_p^p$ 满足如下p个线性方程：

${\displaystyle \underset{r=1}{\overset{p}{\sum }}r{\alpha }_r^p}=1,{\displaystyle \underset{i=1}{\overset{p}{\sum }}{i}^{2s-1}{\alpha }_r^p}=0\left(s=2,3,\cdots ,p\right),$ (5)

例如， $p=2$ 时，四阶熵守恒通量为

$f^{\text{EC},4}=\frac{4}{3}{f}^{\text{EC}}\left(u_i,u_{i+1}\right)-\frac{1}{6}\left[f^{\text{EC}}\left(u_{i-1},u_{i+1}\right)+f^{\text{EC}}\left(u_i,u_{i+2}\right)\right];$ (6)

$p=3$ 时，六阶熵守恒通量为

$\begin{array}{l}{f}^{\text{EC},6}=\frac{3}{2}{f}^{\text{EC}}\left(u_i,u_{i+1}\right)-\frac{3}{10}\left[f^{\text{EC}}\left(u_{i-1},u_{i+1}\right)+f^{\text{EC}}\left(u_i,u_{i+2}\right)\right]\\ +\frac{1}{30}\left[f^{\text{EC}}\left(u_{i-2},u_{i+1}\right)+f^{\text{EC}}\left(u_{i-1},u_{i+2}\right)+f^{\text{EC}}\left(u_i,u_{i+3}\right)\right].\end{array}$ (7)

2.2. 熵稳定格式

实际物理问题中，当间断发生时熵是增加的，上述通量无法再保证熵守恒。所以，在对应的数学模型中，我们所采用的数值方法必须有一定的熵耗散，具体耗散的大小取决于间断发生时熵增(熵产)的多少。

对于浅水波方程，我们引入熵函数 $U\left(u\right)$ 和熵通量函数 $F\left(u\right)$，且满足 $F^{\prime }{\left(u\right)}^{\text{T}}=U^{\prime }{\left(u\right)}^{\text{T}}{f}^{\prime }\left(u\right)$，并考虑计算区域 $\Omega$ 及其边界 $\partial \Omega$ 上的熵增

$\stackrel{\dot{}}{U}={\displaystyle {\iint }_{\Omega }\left({\partial }_{t}U+{\partial }_{x}F\right)\text{d}x\text{d}t}={\displaystyle {\int }_{\partial \Omega }\left(U\text{d}x-F\text{d}t\right)}.$ (8)

目前对这个积分直接求解是非常困难的，接下来我们讨论离散熵增。

考虑离散区间 $\left[x_i,x_{i+1}\right]$，浅水波方程的半离散格式为

${\partial }_t{u}_i=\frac{f_i-f^{*}}{h_i},{\partial }_t{u}_{i+1}=\frac{f^{*}-f_{i+1}}{h_{i+1}},$ (9)

$f^{*}=f^{*}\left(u_i,u_{i+1}\right)$ 是单元交界面处的数值通量，单元交界面处的熵增 $\stackrel{\dot{}}{U}$ 定义为熵通量F产生的熵增量与离散总熵之和

$\stackrel{\dot{}}{U}=-{\left[v^{\text{T}}f\right]}_{i+1/2}+{\left[v^{\text{T}}\right]}_{i+1/2}{f}^{*}+{\left[F\right]}_{i+1/2},$ (10)

这里 ${[\cdot ]}_{i+1/2}={(\cdot )}_{i+1}-{(\cdot )}_i$。当间断发生时，我们假设 $f^{*}=\frac{1}{2}\left(f_i+f_{i+1}\right)$，则熵增

$\stackrel{\dot{}}{U}=F_{i+1}-F_i-\frac{1}{2}\left(v_{i+1}^{\text{T}}+v_i^{\text{T}}\right)\left(f_{i+1}-f_i\right)={\displaystyle {\int }_{u_i}^{u_{i+1}}\left(v^{\text{T}}-\overline{v^{\text{T}}}\right)}\frac{\text{d}f}{\text{d}u}\text{d}u,$ (11)

其中 $\overline{v^{\text{T}}}=\frac{1}{2}\left(v_i^{\text{T}}+v_{i+1}^{\text{T}}\right)$。考虑到原始变量 $u$ 与熵变量 $v$ 的一一对应关系，我们对上式进行换元，并考虑最

简单的直线积分路径 $v=\left(1-\xi \right)v_i+\xi v_{i+1},0\le \xi \le 1$ [15]，则有

$\begin{array}{c}\stackrel{\dot{}}{U}={\displaystyle {\int }_{v_i}^{v_{i+1}}\left(v^{\text{T}}-\overline{v^{\text{T}}}\right)}\frac{\text{d}f}{\text{d}u}\frac{\text{d}u}{\text{d}v}\text{d}v\\ =\frac{1}{2}{\left[v\right]}_{i+1/2}^{\text{T}}{\displaystyle {\int }_0^1\left(2\xi -1\right)}\frac{\text{d}f}{\text{d}u}\frac{\text{d}u}{\text{d}v}\text{d}\xi {\left[v\right]}_{i+1/2}\\ \approx \frac{1}{2}{\left[v\right]}_{i+1/2}^{\text{T}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\Lambda R^{-1}R{R}^{\text{T}}\text{d}\xi {\left[v\right]}_{i+1/2}}\\ =\frac{1}{2}{\left[v\right]}_{i+1/2}^{\text{T}}{\displaystyle {\int }_0^1\left(2\xi -1\right)R\Lambda R^{\text{T}}\text{d}\xi {\left[v\right]}_{i+1/2}}.\end{array}$ (12)

熵增总是正的，即 $\stackrel{\dot{}}{U}>0$，故我们将熵增量的表达式进一步修正为：

$\stackrel{\dot{}}{U}=\frac{1}{2}{\left[v\right]}_{i+1/2}^{\text{T}}R\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \text{d}\xi }\right|R^{\text{T}}{\left[v\right]}_{i+1/2},$ (13)

所以，界面 $x_{i+1/2}$ 处的熵耗散项可以表示为

$D_{i+1/2}^{\text{EP}}=\frac{1}{2}{R}_{i+1/2}\left|{\displaystyle {\int }_0^1\left(2\xi -1\right)\Lambda \text{d}\xi }\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}.$ (14)

为了避免积分求解困难，我们采用三点Gauss-Legendre求积公式对上式中的积分近似求解，得到

$D_{i+1/2}^{\text{EP}}=\frac{1}{2}{R}_{i+1/2}\left|{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\omega }_k\left(2S_k-1\right)\Lambda \left(S_k{u}_{i+1}+\left(1-S_k\right)u_i\right)}\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2},$ (15)

其中 $S_k$ 是高斯点， ${\omega }_k$ 是 $\left[0,1\right]$ 区间上高斯求积公式的相关权重，

$S_{1,3}=\frac{1}{2}\pm \frac{\sqrt{15}}{10},S_2=\frac{1}{2},{\omega }_{1,3}=\frac{5}{9},{\omega }_2=\frac{8}{9}.$ (16)

当解发生间断时，整个系统的熵不再是守恒的，所以我们在原熵守恒数值通量(3)的基础上减去熵增引起的耗散 $D_{i+1/2}^{\text{EP}}$，得到修正的有间断时的熵守恒数值通量

$f_{i+1/2}^{\text{EC}}=f_{i+1/2}^{\text{EEC}}-D_{i+1/2}^{\text{EP}}.$ (17)

为了满足熵稳定条件，再引入经典的Roe-型耗散项

$D_{i+1/2}^{\text{ES}}=\frac{1}{2}{R}_{i+1/2}\left|{\Lambda }_{i+1/2}\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2},$ (18)

这样，就得到了一类新的求解浅水波方程的熵稳定数值通量

$\begin{array}{c}{f}_{i+1/2}^{\text{ES}}=f_{i+1/2}^{\text{EEC}}-D_{i+1/2}^{\text{EP}}-D_{i+1/2}^{\text{ES}}\\ =f_{i+1/2}^{\text{EEC}}-\frac{1}{2}{R}_{i+1/2}\left|\left|{\Lambda }_{i+1/2}\right|+{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\omega }_k\left(2S_k-1\right)\Lambda \left(S_k{u}_{i+1}+\left(1-S_k\right)u_i\right)}\right|R_{i+1/2}^{\text{T}}{\left[v\right]}_{i+1/2}.\end{array}$ (19)

3. 高阶熵稳定格式

熵稳定格式有效地消除了膨胀激波和伪震荡，但只有一阶精度。高阶熵稳定通量的构造可以先通过对熵变量进行满足符号性质的重构获得高阶的数值耗散项，再将其与高阶的熵守恒通量相结合得到。所以，为了提高格式的精度，这一节中，我们利用改进的三阶SPWENO重构对熵变量v进行重构，从而构造求解浅水波方程的高精度高分辨率的数值格式。

3.1. 重构过程

设 ${\tilde{v}}_i\left(x\right),{\tilde{v}}_{i+1}\left(x\right)$ 分别为单元 $I_i$ 与 $I_{i+1}$ 上熵变量的多项式重构，并记

${\tilde{v}}_{i+1/2}^{-}={\tilde{v}}_i\left(x_{i+1/2}\right),{\tilde{v}}_{i+1/2}^{+}={\tilde{v}}_{i+1}\left(x_{i+1/2}\right),{\left[\tilde{v}\right]}_{i+1/2}={\tilde{v}}_{i+1/2}^{+}+{\tilde{v}}_{i+1/2}^{-}$

根据三阶WENO重构方法，单元交界面 $x_{i+1/2}$ 左右两侧的重构值可分别表示为：

${\tilde{v}}_{i+1/2}^{+}={\tilde{\omega }}_{0,i+1/2}\left(-\frac{v_{i+2}}{2}+\frac{3v_{i+1}}{2}\right)+{\tilde{\omega }}_{1,i+1/2}\left(\frac{v_i}{2}+\frac{v_{i+1}}{2}\right),$ (20)

${\tilde{v}}_{i+1/2}^{-}={\omega }_{0,i+1/2}\left(\frac{v_i}{2}+\frac{v_{i+1}}{2}\right)+{\omega }_{1,i+1/2}\left(-\frac{v_{i-1}}{2}+\frac{3v_i}{2}\right),$ (21)

这里 ${\omega }_0,{\omega }_1,{\tilde{\omega }}_0,{\tilde{\omega }}_1\in \left[0,1\right]$ 为权重系数(忽略下标 $i+1/2$ )，为了满足三阶精度和符号性质的要求，即 $sign\left({\left[\tilde{v}\right]}_{i+1/2}\right)=sign\left({\left[v\right]}_{i+1/2}\right)$，权重系数应该满足如下条件：

${\omega }_0=\frac{3}{4}+2C_1,{\omega }_1=1-{\omega }_0,{\tilde{\omega }}_0=\frac{1}{4}-2C_2,{\tilde{\omega }}_1=1-{\tilde{\omega }}_0,$ (22)

其中 $C_1,C_2\in \left[-\frac{3}{8},\frac{1}{8}\right]$。Fjordholm在文献 [12] 给出了 $C_1,C_2$ 的一种选取方式，但其表达式中存在多个人为

参数，编程复杂，而且数值实验结果表明这种方法并不能完全消除激波附近的伪振荡。本文将 $C_1,C_2$ 的表达式改进为

$C_1\left({\theta }_i^{+},{\theta }_{i+1}^{-}\right)=\{\begin{array}{l}\frac{1}{8}{\theta }_i^{+}>1,\\ 0{\theta }_i^{+}<1,{\theta }_{i+1}^{-}=1或者{\theta }_i^{+}=1,{\theta }_{i+1}^{-}\le 1,\\ -\frac{3}{8}{\theta }_i^{+}<1,{\theta }_{i+1}^{-}\ne 1或者{\theta }_i^{+}=1,{\theta }_{i+1}^{-}>1.\end{array}$ (23)

$C_2\left({\theta }_i^{+},{\theta }_{i+1}^{-}\right)=C_1\left({\theta }_{i+1}^{-},{\theta }_i^{+}\right)，$ (24)

其中 ${\theta }_i^{-}=\frac{{\left[v\right]}_{i+1/2}}{{\left[v\right]}_{i-1/2}}$， ${\theta }_i^{+}=\frac{1}{{\theta }_i^{-}}=\frac{{\left[v\right]}_{i-1/2}}{{\left[v\right]}_{i+1/2}}$。

3.2. 高阶熵稳定格式

引理1 (Fjordholm, [11] )假设 $D_{i+1/2}\ge 0$， $\tilde{f}$ 代表熵守恒通量，如果重构过程满足符号性质，那么，数值通量

$f_{i+1/2}={\tilde{f}}_{i+1/2}-\frac{1}{2}{D}_{i+1/2}{\left[\tilde{v}\right]}_{i+1/2}，$ (25)

是熵稳定的。

为了提高熵稳定格式(19)的精度，我们采用四阶的熵守恒通量(6)，并用上一节中提出的三阶SPWENO方法重构耗散项中的熵变量，从而得到一类高阶的熵稳定通量，记为 $f_{i+1/2}^{\text{ES-SPWNEO}3}$

$f_{i+1/2}^{\text{ES-SPWNEO}3}=f_{i+1/2}^{\text{EC},4}-\frac{1}{2}{R}_{i+1/2}\left|\left|{\Lambda }_{i+1/2}\right|+{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\omega }_k\left(2S_k-1\right)\Lambda \left(S_k{u}_{i+1}+\left(1-S_k\right)u_i\right)}\right|R_{i+1/2}^{\text{T}}{\left[\tilde{v}\right]}_{i+1/2}.$ (26)

4. 数值算例

本文时间方向上采用强稳定的龙格–库塔方法：

$\{\begin{array}{l}{u}_i^{*}=u_i^k+\Delta tL\left(u_i^k\right),\hfill \\ u_i^{**}=\frac{3}{4}{u}_i^k+\frac{1}{4}{u}_i^{*}+\frac{\Delta t}{4}L\left(u_i^{*}\right),\hfill \\ u_i^{k+1}=\frac{1}{3}{u}_i^k+\frac{2}{3}{u}_i^{**}+\frac{2\Delta t}{3}L\left(u_i^{**}\right).\hfill \end{array}$

符号约定：

ES：熵稳定格式，数值通量为(19)；ES-SPWENO3：加入三阶SPWENO型重构的高阶熵稳定格式，数值通量为(26)；Exact：精确解。

算例1高度接近0的溃坝问题，初始条件为

$h\left(x,0\right)=1,u\left(x,0\right)=\{\begin{array}{l}-4x<0\\ 4x>0\end{array}$

采用Neumann边界条件，计算区域为 $\left[-1,1\right]$，取空间网格数为 $N=200$，时间上计算到 $t=0.1$，该问题的初始速度使得水流向两边扩展，导致中间区域水深很小，需要注意的是，在解决此类问题时，水深可以接近0，但不会产生负值。在图1所示的数值结果中，ES格式和ES-SPWENO3格式都可以完成计算，但明显可以看出，与ES格式相比，ES-SPWENO3格式对间断的捕捉效果更佳，计算结果更加接近精确解。

Figure 1. Numerical results for Dam Break problem of Near-Zero Height

图1. 高度接近0的溃坝问题的数值结果

算例2大型溃坝问题，初始条件为

$h\left(x,0\right)=\{\begin{array}{l}2x<0\\ 1.5x>0\end{array},u\left(x,0\right)=0,$

采用Neumann边界条件，计算区域为 $\left[-1,1\right]$，取空间网格数为 $N=200$，时间上计算到 $t=0.4$，从图2给出的数值结果来看，ES格式在稀疏波和激波附近出现了严重的抹平现象，相比之下，ES-SPWENO3

Figure 2. Numerical results for the Dam Break problem 图2. 大型溃坝问题的数值结果

格式能够更准确地捕捉解的结构，体现了其高精度、高分辨率的特点。从总熵随时间的变化图可以看出，与ES格式相比，ES-SPWENO3格式的耗散较少，避免了在间断处、稀疏波的波头和波尾处出现过度抹平的现象。