多维Ripa系统的高阶良平衡保正性加权紧致非线性格式

《Computers & Mathematics with Applications》：High-order well-balanced and positivity-preserving weighted compact nonlinear schemes for the multidimensional Ripa system

　　

**摘要**

本文提出了一种用于二维Ripa系统的高阶加权紧致非线性格式（WCNS），该格式具有良好的平衡性和保正性。我们采用源项分裂技术和界面重建方法来保持静水平衡态。同时，引入一种"排空"时间步技术来确保水深和温度的正性。该方法能够有效处理干湿前沿问题，这也是本格式的主要贡献之一。为保持等压平衡态，我们应用平衡态保持参数来修改原始的Lax-Friedrichs通量。通过重建平衡变量而非守恒变量，并结合与通量离散相匹配的源项离散化，来保持动水平衡态。在Ripa系统的一维和二维经典问题上进行了若干数值实验。数值结果表明，我们所提出的格式具有良好的平衡性和保正性、高阶精度以及优秀的激波捕捉性能。

**引言**

本文旨在为Ripa系统开发一种良好平衡的、保正的、高阶加权紧致非线性格式。Ripa模型被用于模拟洋流运动[1]、[2]、[3]。它将浅水方程与描述温度波动的项耦合在一起。一维Ripa模型的 governing 方程结构如下：

$$\begin{cases} h_t + (hu)_x = 0, \\ (hu)_t + (hu^2 + \frac{1}{2}g\theta h^2)_x = -gh\theta B_x, \\ (h\theta)_t + (hu\theta)_x = 0, \end{cases}$$

其中 $h(x, t)$ 表示水深，$u(x, t)$ 表示深度平均速度，$g$ 为重力常数，$\theta$ 表示势温场，$B(x)$ 表示底部地形函数。底部地形与时间变量无关。

系统(1.1)允许如下形式的稳态解：

$$\begin{cases} (hu)_x = 0, \\ (hu^2 + \frac{1}{2}g\theta h^2)_x = -gh\theta B_x, \\ (hu\theta)_x = 0. \end{cases}$$

(1.2)中描述的稳态解包含两种不同的静水平衡态：

(1) **静水平衡态**：$u = 0$，$w = \text{Const}$，$\theta = \text{Const}$，

(2) **等压平衡态**：$p := \frac{1}{2}g\theta h^2 = \text{Const}$，$u = 0$，$B = 0$，其中 $w = h + B$ 表示水面高度，$p$ 为静水压力。

对于稳态解，传统数值格式往往无法精确捕捉它们，并倾向于引起虚假振荡。因此，必须在离散层面上构造能够保持通量梯度与源项之间平衡的数值格式。已有多种数值格式被设计用于以离散形式保持稳态解。浅水方程的良好平衡格式示例可参见[4]、[5]、[6]、[7]、[8]、[9]、[10]、[11]。

在浅水方程的数值模拟方面，间断Galerkin (DG) 方法已被系统地发展。特别是，Aizinger和Dawson [12] 建立了二维深度积分浅水方程的DG框架，实现了高分辨率的流动和输运模拟。随后，Dawson和Aizinger [13] 将该方法扩展到三维浅水方程，并引入了一种网格平滑策略来处理动态自由表面。此外，Aizinger和Dawson [14] 进一步提供了局部间断Galerkin方法应用于三维浅水流模型的稳定性分析，这是对该复杂方程组任何数值格式的首次分析。

稳态解(1.2)还包含一个光滑的运动稳态解，其形式为：

$$\begin{cases} q := (hu) = \text{Const}, \\ E := \frac{q^2}{2h^2} + g\theta(h + B) = \text{Const}, \\ \theta = \text{Const}, \end{cases}$$

其中 $q$ 表示流量，$E$ 表示能量。当流速 $u = 0$ 时，模型退化为静水平衡态(1.3)。由于能量与水深之间存在的非线性关系，设计一个能够精确保持动水平衡态(1.5)的良好平衡格式极具挑战性。已有多种格式被开发用于处理动水平衡态(1.5)，包括有限体积WENO格式[15]、[16]、通量全局化方法[17]、平衡变量重建方法[18]、[49]、间断Galerkin方法[19]、[20]、[21]、路径守恒格式[22]以及一些其他格式[23]、[24]、[25]。

二维Ripa系统如下：

$$\begin{cases} h_t + (hu)_x + (hv)_y = 0, \\ (hu)_t + (hu^2 + \frac{1}{2}g\theta h^2)_x + (huv)_y = -gh\theta B_x, \\ (hv)_t + (huv)_x + (hv^2 + \frac{1}{2}g\theta h^2)_y = -gh\theta B_y, \\ (h\theta)_t + (hu\theta)_x + (hv\theta)_y = 0. \end{cases}$$

对于二维情况，我们关注(1.6)的两种稳态解，具体为：

(1) **静水平衡态**：$u = 0$，$v = 0$，$w = \text{Const}$，$\theta = \text{Const}$，

(2) **等压平衡态**：$p := \frac{1}{2}g\theta h^2 = \text{Const}$，$u = 0$，$v = 0$，$B = 0$。

此外，已有大量数值方法被提出用于Ripa系统。这些包括带有界面跟踪方法的中心迎风格式[26]、带有界面梯度方法的中心有限体积方法[27]、路径守恒HLLC格式[28]、松弛格式[29]、时空CESE格式[30]、Riemann求解器[31]、动力学通量向量分裂格式[32]、非交错中心格式[18]、[33]、高阶紧致有限差分格式[34]以及其他方法[19]、[35]。由于声速表示为 $\sqrt{g\theta h}$，这些格式必须确保水深和温度的正性。

水深的非负性在浅水方程中起着重要作用，并已被广泛研究。Azerad等人[36]通过将保证最大波速的人工粘性方法与集中质量矩阵以及CFL条件结合，直接为浅水方程保证非负水深。当包含Manning摩擦时，Guermond等人[37]通过一种通量校正输运技术实现二阶保正性，该技术强制执行从低阶格式导出的局部约束。此外，Hajduk和Kuzmin [38] 采用代数通量校正和整体凸限制，并引入地形校正因子来构造非负的水深中间态。该格式通过凸组合在CFL条件下保证水深非负。相比之下，Bollermann等人[4]通过限制出流通量并采用局部排空时间来确保水深的正性，从而避免了在干湿前沿附近需要全局时间步缩减。

Ripa系统源项的处理尤其具有挑战性，因为源项可能是一个非守恒乘积。为了高效处理浅水方程中的非守恒乘积，Audusse等人[6]提出了静水压重建(HR)格式，该格式引入中间底部地形来平滑单元边界处的非守恒乘积，从而实现良好平衡性质。后来，Dong [39] 为带有间断地形和Manning摩擦的浅水方程提出了一种二阶界面重建格式，该格式更准确地捕捉了间断地形上的流动行为，并稳定地保持稳态解。随后，Dong和Qian [49] 将该格式扩展到了Ripa系统。

高阶良好平衡格式已被广泛研究，以提高计算精度和分辨率。Xing和Shu [40] 为浅水方程开发了一种加权本质无振荡(WENO)格式，该格式保持了良好平衡性质并实现了高阶精度。其主要贡献是将源项分解为两个分量，并分别进行离散化。Han和Li [41] 后来扩展了这一方法，为Ripa系统开发了一种高阶良好平衡有限差分WENO格式。此外，Liu等人[42]为Ripa和污染物输运模型系统提出了一种高阶良好平衡的替代加权本质无振荡格式。Britton和Xing [19] 提出了一种高阶DG方法，能够同时保持Ripa系统的静水与动水平衡态。上述高阶良好平衡格式大多基于WENO或DG框架。最近，Cheng等人[43]为预平衡形式的浅水方程开发了一种良好平衡的WCNS。然而，[43]中考虑的浅水方程主要在湿地区域，并未讨论水深正性保持性质。

在本文中，受[39]、[40]、[43]启发，我们为Ripa系统开发了一种良好平衡的、保正的、高阶WCNS。本文的主要贡献如下：

- 我们将界面重建和源项分裂技术扩展到Ripa系统的WCNS框架中，后者具有更复杂的稳态结构。所提出的格式在保持静水平衡态的同时实现了高阶精度。

- 本文的关键贡献是：在WCNS框架内严格证明了水深和温度的非负性。这是通过采用"排空"时间步技术实现的。所提出的格式能够稳定捕捉干湿界面流动，克服了[43]中WCNS仅适用于深水区域的局限性。

- 引入一种平衡态保持参数来维持等压平衡态(1.4)。此外，我们通过重建平衡变量并采用与通量离散完全兼容的源项离散化来保持动水平衡态(1.5)。

本文的组织结构如下。在第2节中，我们简要回顾了用于双曲守恒律的五阶WCNS-JS格式。在第3节中，我们介绍一种良好平衡的、保正的、高阶加权紧致非线性格式，用于保持一维Ripa系统的静水平衡态。在第3.1节中，我们介绍一种界面重建方法，该方法重建底部地形和水深。在第3.3节中，我们证明该格式可以保持静水平衡态。在第3.4节中，我们介绍一种"排空"时间步技术来确保水深和温度的正性。在第4节中，我们引入平衡态保持参数并证明该格式维持等压平衡态。在第5节中，我们使用仔细的源项离散化来保持动水平衡态(1.5)。在第6节中，我们将所构造的格式扩展到二维情况。在第7节中，我们给出该格式在Ripa系统一维和二维情况下的数值结果。

**第五阶WCNS-JS格式**

首先，我们考虑以下一维标量守恒律：

$$u_t + f(u)_x = 0,$$

其中 $u$ 表示守恒量，$f(u)$ 为相应的通量。该格式使用交错通量点和解点的网格。为描述WCNS，将区间 $[a, b]$ 划分为 $N$ 个单元，通量点位于 $x_{j+1/2} = a + j\Delta x$，其中 $j = 0, 1, \ldots, N$，$\Delta x = (b - a)/N$，解点位于单元中心 $x_j = a + (j - 1/2)\Delta x$，其中 $j = 1, 2, \ldots, N$。域边界设置为通量点。

**静水保持格式**

在本节中，我们考虑Ripa系统的一维情形。为构造良好平衡格式，我们将源项 $-gh\theta B_x$ 分解为两个分量 $\theta(\frac{1}{2}gB^2)_x - g\theta(h+B)B_x$，如[40]中所做。因此，原始方程(1.1)可以改写为向量形式：

$$\frac{\partial \mathbf{U}}{\partial t} = -\frac{\partial \mathbf{F}(\mathbf{U})}{\partial x} + \mathbf{S_1}(\mathbf{U}, B) + \mathbf{S_2}(\mathbf{U}, B),$$

其中

$$\mathbf{U} = \begin{pmatrix} h \\ hu \\ h\theta \end{pmatrix}, \quad \mathbf{F}(\mathbf{U}) = \begin{pmatrix} hu \\ hu^2 + \frac{1}{2}g\theta h^2 \\ hu\theta \end{pmatrix},$$

$$\mathbf{S_1}(\mathbf{U}, B) = \begin{pmatrix} 0 \\ \theta(\frac{1}{2}gB^2)_x \\ 0 \end{pmatrix}, \quad \mathbf{S_2}(\mathbf{U}, B) = \begin{pmatrix} 0 \\ -g\theta(h+B)B_x \\ 0 \end{pmatrix}.$$

进一步地，我们使用完全相同的最高阶线性方法来离散通量梯度和源项。

**等压稳态保持格式**

第3节中提出的数值方法能够精确保持静水稳态(1.3)，但无法维持等压稳态(1.4)。为保持等压稳态，我们引入[49]中的平衡态保持参数，记为

$$\varphi_{j+1/2} := \max(\varphi^-_{j+1/2}, \varphi^+_{j+1/2}),$$

其中

$$\varphi^-_{j+1/2} = \max\left(\frac{\min(u^-_{j+1/2}, 0)}{a^-_{j+1/2}}, \frac{-\min((hu)^*_{-j+1/2}, 0)}{\max(|(hu)^*_{-j+1/2}|, 1)}\right),$$

$$\varphi^+_{j+1/2} = \max\left(\frac{\max(u^+_{j+1/2}, 0)}{a^+_{j+1/2}}, \frac{\max((hu)^*_{+j+1/2}, 0)}{\max(|(hu)^*_{+j+1/2}|, 1)}\right),$$

其中 $a^+_{j+1/2}$, $a^-_{j+1/2}$ 由

$$a^+_{j+1/2} = \max(u^-_{j+1/2} + \sqrt{g\theta^-_{j+1/2}})$$

计算得出。

**动水稳态保持格式**

在本节中，我们提出一种用于动水稳态(1.5)的良好平衡格式。需要注意的是，静水稳态(1.3)的格式与动水稳态(1.5)的格式之间的主要区别如下：

(1) 水深的重建；

(2) 源项的离散化。

**二维扩展**

在本节中，我们按维逐个维度将所提出的数值格式扩展到二维Ripa系统。首先，我们将二维Ripa系统(1.6)改写为向量形式：

$$\frac{\partial \mathbf{U}}{\partial t} = -\frac{\partial \mathbf{F}(\mathbf{U})}{\partial x} - \frac{\partial \mathbf{G}(\mathbf{U})}{\partial y} + \mathbf{S_1}(\mathbf{U}, B) + \mathbf{S_2}(\mathbf{U}, B),$$

其中

$$\mathbf{U} = \begin{pmatrix} h \\ hu \\ hv \\ h\theta \end{pmatrix}, \quad \mathbf{F}(\mathbf{U}) = \begin{pmatrix} hu \\ hu^2 + \frac{1}{2}g\theta h^2 \\ huv \\ hu\theta \end{pmatrix},$$

$$\mathbf{G}(\mathbf{U}) = \begin{pmatrix} hv \\ huv \\ hv^2 + \frac{1}{2}g\theta h^2 \\ hv\theta \end{pmatrix},$$

$$\mathbf{S_1}(\mathbf{U}, B) = \begin{pmatrix} 0 \\ \theta(\frac{1}{2}gB^2)_x \\ \theta(\frac{1}{2}gB^2)_y \\ 0 \end{pmatrix}, \quad \mathbf{S_2}(\mathbf{U}, B) = \begin{pmatrix} 0 \\ -g\theta(h+B)B_x \\ -g\theta(h+B)B_y \\ 0 \end{pmatrix}.$$

我们重建变量 $\mathbf{V} = (w, u, v, \theta, B)^\top$ 而非 $\mathbf{U}$ 来保持系统(1.1)的静水稳态(1.3)。

**数值结果**

在本节中，我们进行经典的一维和二维数值实验，以验证所提格式的优良良好平衡性和保正性。这些实验包括验证该格式保持静水、等压和动水稳态能力的测试算例。

**结论**

本文针对一维和二维Ripa系统提出了一种良好平衡且保正的加权紧致非线性格式。我们证明了所提格式能够保持静水平衡态、等压平衡态和动水平衡态。为保持静水平衡态，我们采用源项分裂方法来离散化源项，并使用界面重建方法来重建水深和底部地形。

**利益冲突声明**

作者声明不存在利益冲突。

**作者信息**

毛欣悦 | 董建 | 刘豪 | 钱旭

国防科技大学，理学院，湖南长沙，410073，中国