Ripa方程的高阶平衡有限差分A-WENO格式

High-Order Well-Balanced Finite Difference A-WENO Scheme for the Ripa Equations

摘要: 在本文中，我们提出了一种基于自适应加权本质无振荡(Adaptive Weighted Essentially Non-Oscillatory, A-WENO)方法的高阶有限差分格式，用于求解在浅水方程中引入水平温度梯度后形成的Ripa方程。该方程属于带源项的非齐次双曲守恒律，其稳态解在通量梯度与源项之间存在严格平衡。所提出的数值格式在有限差分框架下，通过对源项进行特殊分裂、结合改进的Lax-Friedrichs通量分裂及特征分解，并将由WENO重建得到的离散导数系数“冻结”后应用于源项离散，从而在保持五阶空间精度的同时，实现了严格的稳态保持性质。与传统的高阶WENO格式相比，本文方法在复杂底形和强非线性条件下能够有效抑制数值振荡，并显著提高了计算的稳定性与精度。数值算例验证了该方法在保持平衡、高阶精度及捕捉间断方面的优越性能。

Abstract: In this paper, a high-order finite difference scheme based on the Adaptive Weighted Essentially Non-Oscillatory (A-WENO) method is proposed for solving the Ripa equations, which are derived from the shallow water equations by introducing horizontal temperature gradients. The Ripa system represents a class of nonhomogeneous hyperbolic conservation laws with source terms, whose steady-state solutions exhibit an exact balance between the flux gradients and the source terms. Within the finite difference framework, the proposed method reformulates the source term, incorporates a modified Lax–Friedrichs flux splitting and characteristic decomposition, and applies the frozen discrete derivative coefficients obtained from WENO reconstruction to the source term discretization. This approach achieves strict well-balanced properties while maintaining fifth-order spatial accuracy. Compared with conventional high-order WENO schemes, the proposed A-WENO method effectively suppresses numerical oscillations under complex topographies and strong nonlinear conditions, enhancing both stability and accuracy. Numerical experiments demonstrate that the method exhibits excellent well-balanced behavior, high-order accuracy, and strong capability in resolving discontinuities.

1. 引言

浅水方程作为描述浅水区自由水面流动的数学模型，在模拟河流、渠道等水动力学问题中具有重要的应用价值。该方程在水利工程和海岸工程的数值模拟中得到广泛应用[1] [2]，而Ripa方程在海洋学与环境流体力学中具有重要应用，能够描述密度分层流体的运动。由于地形源项的存在，标准的高阶格式往往在稳态附近产生非物理的寄生振荡，从而污染微小扰动波的传播。开发高阶且能精确保持平衡的数值格式，对于长时间模拟以及捕捉小振幅扰动波至关重要。该模型在[3]-[5]中提出，用于模拟海洋洋流。其一维形式可写作带源项的非齐次双曲守恒律方程：

其中$h,$ $u,$ $\theta$ 分别表示水深、流速和温度，$g$ 为重力加速度常数，$B$ 表示底部地形函数；$h+B$ 对应自由水面高度，而$\frac{1}{2}g{h}^2\theta$ 表示压力项。该模型可保持如下两类稳态解：

$u=0,B=常�,h^2\theta =常�.$

或

$u=0,h+B=常�,\theta =常�.$

在稳态下，通量梯度与源项严格平衡。Bermudez与Vazquez [6]于1994提出“精确守恒性”(well-balanced)以刻画数值格式与稳态解的完全相容；自Greenberg等人[7] [8]开创性工作以来，相关方法统称well-balanced格式。近年来，well-balanced格式的研究活跃，已形成多种构造思路，可参见文献[9]-[14]。因而，发展兼具高阶精度与平衡保持能力的有限差分格式具有重要意义。对Ripa方程而言，若离散层面不能维持平衡，即使解析解为稳态也会出现非物理振荡，且多维中难以靠加密网格抑制误差。为克服这一困难，本文在有限差分框架下提出了一种高阶自适应加权本质无振荡(Adaptive Weighted Essentially Non-Oscillatory, A-WENO)格式，其主要思想包括：

源项分裂：通过代数恒等式：

$-gh\theta B_x=\theta {\left(\frac{1}{2}g{B}^2\right)}_x-g\theta \left(h+B\right)B_x,$

将源项重构为可由同一差分算子离散的形式；

修改的Lax-Friedrichs通量分裂：引入修正变量

$\tilde{U}={\left(h+B,hu,\left(h+B\right)\theta \right)}^{\text{T}}$

在稳态下该变量为常向量，可实现黏性项的逐点抵消；

同一算子离散：采用同一差分算子离散通量与源项，保证离散层面的严格平衡。

该方法在保持五阶空间精度的同时，能精确维持稳态解并抑制数值振荡。数值实验验证了其在保持平衡、高阶精度以及复杂底形下稳定性方面的优越性。

2. 高阶有限差分A-WENO格式

在本节中，我们回顾下有限差分A-WENO格式[15] [16]，并对一些标准符号加以说明。区间$I=\left[a,b\right]$ 被分成$N$ 个子区间并且我们用$I_j=\left[x_{j-1/2},x_{j+1/2}\right],j=1,\cdots ,N$ 来表示单元。为方便起见，我们先假设网格均匀分布，单元格的大小为$\Delta x=x_{j+1}-x_j$ ，并且我们使用$x_{j+1/2}=x_j+\Delta x/2$ 来表示单元$I_i$ 的中心。在不丧失一般性的情况下，采用有限差分方法对空间变量进行离散，Ripa方程可表示为：

$u_t+f{\left(u\right)}_x=0$

半离散格式为

$\frac{\text{d}{U}_j}{\text{d}t}=-\frac{1}{\Delta x}\left({\widehat{f}}_{j+\frac{1}{2}}-{\widehat{f}}_{j-\frac{1}{2}}\right)+S_j$

其中${\widehat{f}}_{j+\frac{1}{2}}$ 和${\widehat{f}}_{j-\frac{1}{2}}$ 为数值通量。A-WENO方法旨在利用高阶插值对数值通量进行重构，以便在解的光滑区域保持高阶精度，并在间断附近避免产生非物理振荡。进一步地，在有限差分WENO框架中，界面数值通量由通量分裂后的点值序列通过WENO非线性重构得到。对分裂通量点值序列$\left\{f_j^{\pm }\right\}$ ，A-WENO在界面$x_{j+1/2}$ 给出${\widehat{f}}_{j+\frac{1}{2}}^{\pm }={ℛ}_{\text{WENO}}\left(f_{j-r}^{\pm },\cdots ,f_{j+r}^{\pm }\right)$，${\widehat{f}}_{j+\frac{1}{2}}={\widehat{f}}_{j+\frac{1}{2}}^{+}+{\widehat{f}}_{j+\frac{1}{2}}^{-}$ 进而采用通量差分形式近似通量导数：${\left(D_{f}f\right)}_j=\frac{1}{\Delta x}\left({\widehat{f}}_{j+1/2}-{\widehat{f}}_{j-1/2}\right)$ 。

本文所称的“源项离散算子$D_f$ ”即与通量导数共享的通量差分算子。构造良好平衡格式的关键在于：对源项中可改写为空间导数的部分，同样采用该$D_f$ 进行离散(见第3.3节)，从而在所考虑的静水平衡解上，使离散层面的通量梯度与源项项在网格点处严格抵消。需要强调的是，本文采用纯有限差分的点值框架：离散未知量$U_j$ 表示网格点$x_j$ 处的点值(point-wise values)，而非单元平均值(cell averages)。下文中所提到的“重构/插值”是指对点值序列$\left\{U_j\right\}$ (或分裂通量点值$\left\{f^{\pm }\left(U_j\right)\right\}$ )进行A-WENO非线性重构，在界面$x_{j+1/2}$ 构造界面数值通量${\widehat{f}}_{j+\frac{1}{2}}^{\text{num}}$(若采用状态型实现，则先得到$U_{j+\frac{1}{2}}^{L/R}$ 再构造通量)。因此，后续静水重构亦仅作用于界面重构值，对界面左右量进行修正后再计算数值通量，并最终通过通量差分形式给出导数近似。

2.1. 五阶WENO插值

对于一维标量守恒律

$u_t+f{\left(u\right)}_x=0$

A-WENO方法通过在单元边界$x_{j+\frac{1}{2}}$ 处进行加权插值来构造数值通量。以左侧插值${\widehat{f}}_{j+\frac{1}{2}}^L$ 为例，定义三个三点子模板：

$S_0=\left\{x_{j-2},x_{j-1},x_j\right\},S_1=\left\{x_{j-1},x_j,x_{j+1}\right\},S_2=\left\{x_j,x_{j+1},x_{j+2}\right\}$

每个子模板构造一个三阶插值多项式，其在边界$x_{j+\frac{1}{2}}$ 的值记为$q_k\left(x_{j+\frac{1}{2}}\right)$ 。

五阶A-WENO重建为它们的加权和：

${\widehat{f}}_{j+\frac{1}{2}}^L={\displaystyle \sum _{k=0}^2{\omega }_k{q}_k\left(x_{j+\frac{1}{2}}\right)}$

其中${\omega }_k$ 为非线性权重，反映各子模板的光滑性。

非线性权重定义为

${\alpha }_k=C_k{\left(\epsilon +I_k\right)}^{-p},{\omega }_k=\frac{{\alpha }_k}{{\displaystyle \sum _{\ell =0}^2{\alpha }_{\ell }}}$

其中线性权重$C_k$ 由最小二乘条件确定，满足${\displaystyle \sum C_k}=1$ 。

其中$I_k$ 为光滑度指标：

$I_k={\displaystyle \sum _{m=1}^2{\displaystyle {\int }_{x_{j-1}}^{x_{j+1}}{\left(\frac{{\text{d}}^m{q}_k}{\text{d}{x}^m}\right)}^2\text{d}x}}$

$\epsilon$ 为防止除零的小常数，$p$ 控制非线性权重对光滑度的敏感性，$m$ 为导数阶数。光滑度指标$I_k$ 越小表示子模板越光滑，其权重越大。其右侧插值${\widehat{f}}_{j+\frac{1}{2}}^R$ 的构造与左侧完全对称。最终的界面数值通量可由双边平均给出。

2.2. 生成数值通量

为了在求解Ripa方程时兼顾高阶精度与保持平衡特性，本节采用改进的Lax-Friedrichs通量分裂与 A-WENO非线性重建相结合的方法。对于Ripa方程而言，特征速度主要由流速$u$ 与重力项$\sqrt{g\theta h}$ 决定，因此可取

$\alpha =\left|u\right|+\sqrt{g\theta h}$

这样处理可保证数值通量的稳定性与无振荡性。然而，在稳态$u=0,h+B=常�,\theta =常�$ 下，若直接采用$U$ 参与通量分裂，由于$h$ 随地形$B$ 变化，$\alpha U$ 项中仍包含空间变化，导致通量残差不为零，破坏了稳态平衡。为避免此问题，本文引入修正变量：

$\tilde{U}={\left(h+B,hu,\left(h+B\right)\theta \right)}^{\text{T}}$

并将通量分裂修改为

$f^{\pm }\left(U\right)=\frac{1}{2}\left[f\left(U\right)\pm \alpha \tilde{U}\right]$

此时，在平衡态下$\tilde{U}$ 为常向量，所有与$\alpha \tilde{U}$ 相关的项在离散层面可逐点抵消，实现真正的“保持平衡(well-balanced)”特性。为了保证高阶精度，A-WENO重建过程采用五阶中心差分模板，通过非线性权重自动调节通量插值。最终，通量导数由统一差分算子$D_f$ 表示为：

${\left(D_{f}f\right)}_j=\frac{{\widehat{f}}_{j+\frac{1}{2}}-{\widehat{f}}_{j-\frac{1}{2}}}{\Delta x}$

由于源项也会使用同一个差分算子进行离散(详见3.3节中的“统一算子离散”)，故可以实现在离散层面对稳态平衡关系的逐点匹配。

2.3. 推广到方程组

对于Ripa方程这类多变量系统，A-WENO格式可在特征空间中进行重建，从而更好地刻画不同波模的传播特性。具体而言，首先在单元界面处利用Roe平均状态构造相应的雅可比矩阵：

并对其进行特征分解，可求得其右特征向量矩阵$R$ 和左特征向量矩阵$L=R^{-1}$ ，通过这一过程，可以将守恒变量及其通量表示为一组相互独立的特征模态。在此基础上，将守恒变量投影到特征空间：

$V=LU$

然后在每个特征分量上独立执行A-WENO插值与通量重建。最后将重建值投影回物理$L=R^{-1}$ 空间：

这种“特征分解–分量重建–物理回投影”的策略，能够在离散层面有效解耦不同物理波动，显著减弱变量间的数值耦合误差，避免直接在物理变量上重建时容易出现的伪波和振荡现象。

3. 良好平衡的A-WENO格式

3.1. 改写方程

为在离散层面保持Ripa方程静水平衡，需对方程作代数分裂，使稳态下通量与源项逐点精确抵消。本文给出等价重写并据此构造良好平衡有限差分A-WENO格式；针对一维非齐次形式直接离散会破坏静水稳态，故首先对源项分裂：

据此可将方程改写为

$U_t+\tilde{f}{\left(U\right)}_x=\tilde{S}\left(U,x\right),$

其中

$\tilde{f}\left(U\right)=\left[\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g\left(h^2+2hB\right)\theta \\ hu\theta \end{array}\right],\tilde{S}\left(U,x\right)=\left[\begin{array}{c}0\\ -g\theta \left(h+B\right)B_x\\ 0\end{array}\right].$

这样在稳态条件下，${\tilde{f}}_x$ 与$\tilde{S}$ 可在离散层面逐点匹配。通过本节的代数变换，Ripa方程被改写为通量项与源项可由同一差分算子离散的等价形式，从而为后续两节的数值格式构造与分析奠定了基础。

为避免有限体积与有限差分表述上的混淆，下面给出Audusse等人的静水重构在本文有限差分通量差分框架中的实现要点(以界面$x_{j+1/2}$ 为例)。首先，我们对点值序列$\left\{h_j,b_j,{\left(hu\right)}_j,{\left(h\theta \right)}_j\right\}$ 施加A-WENO非线性重构，得左右重构值$h_{j+\frac{1}{2}}^{L/R},b_{j+\frac{1}{2}}^{L/R},{\left(hu\right)}_{j+\frac{1}{2}}^{L/R},{\left(h\theta \right)}_{j+\frac{1}{2}}^{L/R}$ ，并组装左右状态$U_{j+\frac{1}{2}}^L,U_{j+\frac{1}{2}}^R$ 。随后我们采用静水修正：取$b_{j+\frac{1}{2}}^{*}=\max \left(b_{j+\frac{1}{2}}^L,b_{j+\frac{1}{2}}^R\right)$ ，$h_{j+\frac{1}{2}}^{L/R,*}=\max \left(0,h_{j+\frac{1}{2}}^{L/R}+b_{j+\frac{1}{2}}^{L/R}-b_{j+\frac{1}{2}}^{*}\right)$ ，并按比例缩放$hu,h\theta$ 以保持界面速度与温度一致性，从而得到修正后的左右状态$U_{j+\frac{1}{2}}^{L,*},U_{j+\frac{1}{2}}^{R,*}$ 。将其代入本文采用的“修正变量 + 改进Lax-Friedrichs分裂”构造界面数值通量${\widehat{F}}_{j+\frac{1}{2}}=F^{\text{num}}\left(U_{j+\frac{1}{2}}^{L,*},U_{j+\frac{1}{2}}^{R,*}\right)$，并仍以通量差分骨架更新：${\dot{U}}_j=-\frac{1}{\Delta x}\left({\widehat{F}}_{j+\frac{1}{2}}-{\widehat{F}}_{j-\frac{1}{2}}\right)+S_j=-{\left(D_f\widehat{F}\right)}_j+S_j$ ，其中$D_f$ 为与通量导数共享的差分算子。与第3.3节的源项一致性离散配合时，可在静水平衡解上实现离散层面的严格抵消，从而构造良好平衡格式。

3.2. 应用静水重构方法构造数值通量

为保持静水平衡，我们引入修正变量：

$\tilde{h}=h+B,\tilde{U}={\left(\tilde{h},hu,\tilde{h}\theta \right)}^{\text{T}}$

在稳态下$\tilde{U}$ 为常向量。通量采用修改的Lax–Friedrichs分裂：

在单元边界$x_{j+\frac{1}{2}}$ 处，用静水重构法[14] [17]修正左右水深：

$h_{j+\frac{1}{2},L}^{*}=\max \left(0,h_j+B_j-\max \left(B_j,B_{j+1}\right)\right),h_{j+\frac{1}{2},R}^{*}=\max \left(0,h_{j+1}+B_{j+1}-\max \left(B_j,B_{j+1}\right)\right)$

对应的界面通量为

${\widehat{f}}_{j+\frac{1}{2}}=\frac{1}{2}\left[f\left(U_{j+\frac{1}{2},L}^{*}\right)+f\left(U_{j+\frac{1}{2},R}^{*}\right)-{\alpha }_{j+\frac{1}{2}}\left({\tilde{U}}_{j+\frac{1}{2},R}^{*}-{\tilde{U}}_{j+\frac{1}{2},L}^{*}\right)\right]$

静水重构结合修正变量与Lax–Friedrichs分裂，使稳态下界面通量逐点抵消、黏性项无残差；并在干湿界面避免负水深，可增强复杂底形鲁棒性。

3.3. 源项离散

本节的核心在于：将源项中可改写为空间导数的平衡量，与对流通量导数采用完全一致的差分算子$D_f$ 进行离散，从而在静水平衡附近实现离散层面的精确抵消。为此，我们利用有限差分WENO在界面重构中产生的局部线性组合结构：界面数值量$q_{j\pm \frac{1}{2}}^{\text{num}}$ 由A-WENO作用于点值序列得到，一旦界面处的非线性权重$\left\{{\omega }_k\right\}$ 确定，$q_{j\pm \frac{1}{2}}^{\text{num}}$ 以及通量差分导数${\left(D_{f}q\right)}_j$ 都可等价写成邻域点值的线性组合。为明确“同一算子离散”的含义，我们将通量导数离散算子记为$D_f$ ，其基本形式为${\left(D_{f}q\right)}_j=\frac{1}{\Delta x}\left(q_{j+\frac{1}{2}}^{\text{num}}-q_{j-\frac{1}{2}}^{\text{num}}\right)$ 。

其中$q_{j\pm \frac{1}{2}}^{\text{num}}$ 由A-WENO从点值序列重构得到。本文所谓“冻结系数”是指：先在通量端由$\left\{f\left(U_j\right)\right\}$ 的重构确定每个界面$x_{j\pm \frac{1}{2}}$ 处的WENO权重，随后在同一界面将这组权重原封不动复用于源项端，对源项中导数型平衡量$G\left(U,b\right)$ 做同阶重构得到${\widehat{G}}_{j\pm \frac{1}{2}}$ ，并用同一算子$D_f$ 计算${\left(D_{f}G\right)}_j=\frac{1}{\Delta x}\left({\widehat{G}}_{j+\frac{1}{2}}-{\widehat{G}}_{j-\frac{1}{2}}\right)$ 。由此，源项离散与通量导数离散在模板与权重层面严格同构，使得静水平衡附近通量梯度与源项项可在离散层面逐点抵消，从机制上抑制因离散不一致引发的非物理寄生振荡。

依据上述分裂，将动量方程源项写为：

接下来，对流量通量和源项使用统一差分算子进行离散。在一维情况下，采用高阶差分方法(如五阶WENO或中心差分)进行通量和源项的离散化。对流量通量的离散表示为：

$\frac{\partial }{\partial x}\left(h{u}^2+\frac{1}{2}g{h}^2\right)=D_x\left(h{u}^2+\frac{1}{2}g{h}^2\right)$

对于流量通量的离散，可以表示为：

$\frac{\partial }{\partial x}\left(h{u}^2+\frac{1}{2}g{h}^2\right)$

对于源项的离散，假设源项形式为：

$\frac{\partial }{\partial x}\left(g{h}^2\right)$

该“同一算子离散”保证通量与源项在离散层面的形式一致。

故源项的离散式表达为：$S_j=\left[{\theta }_j{D}_j{\left(\frac{1}{2}g{B}^2\right)}_j-g{\theta }_j\left(h_j+B_j\right)D_f{\left(B\right)}_j\right]$

这样源项的形式与通量的离散层面一致，确保了源项与通量在数值离散时的平衡性。

3.4. 证明良好平衡性质

我们可设Ripa方程的离散形式为

$\frac{\text{d}{U}_j}{\text{d}t}=-\frac{1}{\Delta x}\left({\widehat{f}}_{j+\frac{1}{2}}-{\widehat{f}}_{j-\frac{1}{2}}\right)+S_j,$

其中

${\widehat{f}}_{j+\frac{1}{2}}=\frac{1}{2}\left[f\left(U_{j+\frac{1}{2},L}\right)+f\left(U_{j+\frac{1}{2},R}\right)-{\alpha }_{j+\frac{1}{2}}\left({\tilde{U}}_{j+\frac{1}{2},R}-{\tilde{U}}_{j+\frac{1}{2},L}\right)\right]$

$S_j=\left[\begin{array}{c}0\\ {\theta }_j{D}_f{\left(\frac{1}{2}g{B}^2\right)}_j-g{\theta }_j\left(h_j+B_j\right)D_f{\left(B\right)}_j\\ 0\end{array}\right]$

静水稳态为

$u=0,h+B=H_0,\theta ={\theta }_0,$

由此有

$\tilde{U}=\left[\begin{array}{c}h+B\\ hu\\ \left(h+B\right)\theta \end{array}\right]=\left[\begin{array}{c}{H}_0\\ 0\\ H_0{\theta }_0\end{array}\right]=常向量,f\left(U\right)=\left[\begin{array}{c}0\\ \frac{1}{2}g{h}^2{\theta }_0\\ 0\end{array}\right]$

在稳态下

${\tilde{U}}_{j+\frac{1}{2},R}={\tilde{U}}_{j+\frac{1}{2},L}\Rightarrow {\alpha }_{j+\frac{1}{2}}\left({\tilde{U}}_{j+\frac{1}{2},R}-{\tilde{U}}_{j+\frac{1}{2},L}\right)=0$

于是

${\widehat{f}}_{j+\frac{1}{2}}=f{\left(U\right)}_{j+\frac{1}{2}},{\widehat{f}}_{j+\frac{1}{2}}-{\widehat{f}}_{j-\frac{1}{2}}=D_f\left(f\left(U\right)\right)$

取动量方程第二分量：

${\left(D_f\left(\frac{1}{2}g{h}^2{\theta }_0\right)\right)}_j-{\theta }_0{\left(D_f\left(\frac{1}{2}g{B}^2\right)\right)}_j+g{\theta }_0\left(h_j+B_j\right)D_f{\left(B\right)}_j$

由

$\frac{1}{2}g{h}^2=\frac{1}{2}g{\left(H_0-B\right)}^2=\frac{1}{2}g\left(H_0^2-2H_{0}B+B^2\right),$

代入上式得

$D_f\left(\frac{1}{2}g{H}_0^2{\theta }_0\right)-g{\theta }_0{H}_0{D}_f\left(B\right)+\frac{1}{2}g{\theta }_0{D}_f\left(B^2\right)-{\theta }_0{D}_f\left(\frac{1}{2}g{B}^2\right)+g{\theta }_0\left(H_0-B\right)D_f\left(B\right)=0,$

整理后所有项完全抵消：

$D_f\left(\frac{1}{2}g{\theta }_0{\left(H_0-B+B\right)}^2\right)=0\Rightarrow \frac{\text{d}{U}_j}{\text{d}t}=0$

这意味着格式是良好平衡的。

3.5. 平衡格式的总结

最后，我们针对Ripa方程的数值求解，给出良好平衡有限差分A-WENO格式的完整过程。A-WENO格式详细的算法如下所示：

1) 我们对源项进行代数分裂，通过合适的分裂形式将源项改写为能与通量使用相同差分算子进行离散的部分以及其他部分，应用高阶A-WENO方法对数值通量进行重建；

2) 我们对改写后的源项中的导数项进行近似处理。

3) 为了推进时间，我们采用龙格–库塔方法[18]对数值通量和源项进行时间推进，最后重复步骤2至步骤5，直到达到预定的终止条件(如最大时间步数或误差范围)。

4. 数值结果

在本节中，我们通过进行广泛的数值实验，展示了五阶(r = 2)有限差分A-WENO方案的性能。所有实验中，CFL数值设定为0.6，采用$\epsilon ={10}^{-6}$ 且重力常数g取值为1。为了优化时间步长，时间步长被修改为$\Delta t=\text{CFL}\cdot \left(\Delta x^{5/3}/\alpha \right)$ ，而非传统的$\Delta t=\text{CFL}\cdot \left(\Delta x/\alpha \right)$ ，从而减小了时间步长，其中$\alpha =\max \left|\lambda \left(u\right)\right|$ ，$\lambda \left(u\right)$ 是由雅可比矩阵的特征值。

4.1. 测试平衡性质

首先，测试了所采用A-WENO方案的平衡性质。我们定义了两种不同的底部地形，分别是光滑底部和不连续底部。光滑底部由以下公式给出：

$B\left(x\right)=5\exp \left[-\frac{2}{5}{\left(x-5\right)}^2\right]$

而不连续底部则定义为：

$B\left(x\right)=\{\begin{array}{ll}4,\hfill & 4\le x\le 8,\hfill \\ 0,\hfill & 其他情�.\hfill \end{array}$

在这些底部地形下，初始条件被设置为和。本算例取θ为常数(即无水平温度梯度)。更一般地，当存在“温度间断 + 变化地形”时，θ的跃迁会通过斜压压力项改变局部平衡，可能产生内部重力波并在坡折处发生物理反射/透射。为了验证该A-WENO方案是否能够维持平衡性质，我们使用不同的精度进行数值模拟，模拟时间为t = 0.5，网格大小为200个单元。计算不同精度下和的误差，见表1和表2。

Table 1. The total error of the steady-state solution on the slippery surface under different precisions

表1. 不同精度下在光滑底面地形上稳态解的和误差

精度	L1误差				L∞误差
	$h$	$hu$	$h\theta$		$h$	$hu$	$h\theta$
单精度	2.17e−07	3.13e−07	3.35e−07		9.54e−07	1.05e−07	6.12e−07
双精度	8.70e−16	6.34e−16	2.61e−16		7.11e−15	2.19e−16	1.07e−16
四倍精度	5.01e−33	4.32e−31	9.15e−32		3.08e−33	1.42e−31	4.57e−31

Table 2. The total error of the steady-state solution on the discontinuous surface under different precisions

表2. 不同精度下在不连续底面地形上稳态解的和误差

精度	L1误差				L∞误差
	$h$	$hu$	$h\theta$		$h$	$hu$	$h\theta$
单精度	7.15e−08	2.37e−07	1.50e−07		9.54e−07	5.01e−07	1.49e−08
双精度	2.66e−16	2.12e−16	9.11e−15		3.55e−15	1.07e−16	5.27e−16
四倍精度	1.06e−33	1.46e−31	2.35e−32		5.16e−33	1.44e−31	5.33e−32

4.2. 测试精度阶数

接下来，我们测试了该A-WENO方案在光滑解条件下的五阶精度。我们选择了一个简单的底部地形：$B\left(x\right)={\sin }^2\left(\pi x\right)$ ，定义了计算区域为[0, 1]，初始条件如：

$h\left(x,0\right)=5+\exp \left(\cos \left(2\pi x\right)\right),\left(hu\right)\left(x,0\right)=\sin \left(\cos \left(2\pi x\right)\right),\theta \left(x,0\right)=\cos \left(2\pi x\right).$

边界条件为周期性边界，计算时间t = 0.1，采用五阶A-WENO方案，网格大小为6400，计算L1误差和精度阶数。本节算例底形为平坦($b\equiv 0$ )，温度间断的演化主要体现对流输运以及其通过压力项对动量波系的调制，不涉及由地形引起的散射或反射。$h,hu$ 和$h\theta$ 的L1误差和精度阶数，见表3。

Table 3. The error and convergence order in Section 4.2

表3. 在4.2节算例中，和的L1误差和精度阶数

网格大小	$h$			$hu$			$h\theta$
	L1误差	阶数		L1误差	阶数		L1误差	阶数
25	1.7566E−02			1.0990E−01			2.0208E−02
50	2.2028E−03	3.00		1.9714E−02	2.48		2.7216E−03	2.89
100	3.3138E−04	2.73		2.8273E−03	2.80		3.5233E−04	2.95
200	2.3271E−05	3.83		2.0103E−04	3.81		2.3724E−05	3.89
400	9.3899E−07	4.63		8.1350E−06	4.63		9.5787E−07	4.63
800	3.1516E−08	4.90		2.7319E−07	4.90		3.1021E−08	4.95
1600	8.5864E−10	5.20		7.4405E−09	5.20		8.2517E−10	5.23

4.3. 平底上的黎曼问题

我们进一步利用平底地形($b\equiv 0$ )下的黎曼问题来检验本文五阶A‑WENO格式在捕捉间断、分辨波结构以及保持无振荡性方面的性能。这里采用水坝破裂(dam‑break)型初值作为测试。

$\left(h,u,\theta \right)\left(x,0\right)=\{\begin{array}{ll}\left(5,0,3\right),\hfill & x<0,\hfill \\ \left(1,0,5\right),\hfill & x\ge 0.\hfill \end{array}$

Ripa方程中的温度变量θ通过压力项引入斜压效应：在一维情形中，θ主要随流场对流输运，其跃迁通常表现为随流移动的接触间断，图1中θ的跃迁位置与水深/压力波系结构保持一致，且间断附近未出现明显高频振荡，说明本文对温度梯度的离散能够在保持高阶精度的同时避免伪波产生。我们计算在时间t = 0.2时的数值结果，并将其与参考解进行了对比，见图1。

(a)

(b) (c)

Figure 1. Numerical solution of the flat-bottomed Riemann problem in Section 4.3 at t = 0.2 with 200 grid cells. (a) Water depth; (b) Temperature; (c) Pressure

图1. 4.3节中平底黎曼问题在t = 0.2时，200个网格单元下的数值解。(a) 水深；(b) 温度；(c) 压力

4.4. 平底上的水坝破裂问题

最后，我们进一步探讨了平底上的水坝破裂问题。在这一测试中，初始条件如下：

$\left(h,u,\theta \right)\left(x,0\right)=\{\begin{array}{ll}\left(2,0,1\right),\hfill & \left|x\right|\le 0.5,\hfill \\ \left(1,0,1.5\right),\hfill & 其他情�\text{.}\hfill \end{array}$

计算区域为[−1, 1]。该格式对温度相关项具有良好的稳定性。本算例在含不连续底形的稳态保持算例中未观察到此类伪波，侧面验证了所提方法对地形–温度耦合的鲁棒性。我们的数值计算在时间t = 0.2时进行，并与参考解进行了对比，见图2。

(a) (b)

(c) (d)

Figure 2. Numerical results for the flat-bottomed dam-break problem in the example of Section 4.4. (a) Water depth; (b) Discharge; (c) Temperature; (d) Pressure

图2. 4.4节算例平底水坝破裂问题的数值结果。(a) 水深；(b) 流量；(c) 温度；(d) 压力

5. 结论与展望

本文针对含源项的Ripa方程，提出一种严格平衡的高阶有限差分A-WENO格式。通过源项代数分裂、修正变量及统一差分离散，使通量梯度与源项在离散层面逐点匹配，严格保持静水稳态。结合自适应权重与改进Lax-Friedrichs分裂，数值结果表明，该格式在光滑与间断解中均具高精度、稳定性与非振荡性，较传统WENO在波结构分辨和稳态保持方面更具优势。同时，提出的高阶良好平衡有限差分A-WENO格式可自然推广至二维Ripa方程。二维情形下，可在x，y两个方向分别构造同阶的通量差分算子$D_f^x,D_f^y$ ，并将地形源项中的${\partial }_{x}b,{\partial }_{y}b$ 按“同一算子离散”原则分别与对应方向的通量导数相匹配。与此同时，可在每个网格面上实施与一维完全一致的静水重构：以$b^{*}=\max \left(b^L,b^R\right)$ 修正水深h，并在源项平衡量的界面重构中复用该面的重构权重/线性组合系数，从而在二维中延续一维格式的严格平衡性质与正性保持。另一方面，二维情况下温度梯度将引入更丰富的斜压效应与内部波结构。针对“温度间断跨越复杂地形”这一典型场景，后续可在数值实验中比较坡折前后的波幅与能量变化(例如定义反射/透射系数，或对自由面与速度扰动进行谱分析)，并结合网格加密与高精度参考解对照，系统地区分真实的物理反射/透射与由离散不平衡引发的数值伪反射，从而进一步验证方法在复杂二维地形与热力耦合问题中的鲁棒性与可靠性。

基金项目

本研究得到了山东省自然科学基金面上项目(No. ZR2021MA072, ZR2023MA012)的资助。

NOTES

^*通讯作者。