线性浅水离散格式的能量守恒条件

约 2515 字大约 8 分钟

2026-05-28

线性浅水方程的能量守恒不能只看空间离散，也不能只看时间格式。一个数值格式要真正保持能量，至少要同时满足两层结构：

空间离散中，离散梯度和离散散度互为负伴随；
时间推进中，更新格式能保持反对称线性系统的二次能量。

这篇文章统一说明连续系统、半离散矩阵、Crank--Nicolson、forward Euler，以及 overlap 权重里的常见陷阱。

连续系统里的能量守恒

考虑闭曲面上的线性旋转浅水系统

\partial_t \bm{u}+f\bm{u}^\perp+g\nabla h=0,

\partial_t h+H\nabla\cdot\bm{u}=0.

这里 $\bm{u}$ 是速度， $h$ 是自由面高度扰动， $H$ 是平均流体深度。物理能量定义为

E(t)=\frac12\int_\Omega \left(H|\bm{u}|^2+gh^2\right)\,dx.

对时间求导并代入方程：

\frac{dE}{dt} = H(\partial_t \bm{u},\bm{u}) +g(\partial_t h,h).

再利用

\partial_t \bm{u}=-f\bm{u}^\perp-g\nabla h, \qquad \partial_t h=-H\nabla\cdot\bm{u},

得到

\frac{dE}{dt} = -H(f\bm{u}^\perp,\bm{u}) -gH(\nabla h,\bm{u}) -gH(h,\nabla\cdot\bm{u}).

其中

(\bm{u}^\perp,\bm{u})=0,

而在没有边界通量时，

(\nabla h,\bm{u})=-(h,\nabla\cdot\bm{u}).

两项抵消，所以

\frac{dE}{dt}=0.

因此连续能量守恒的核心有两点：

Coriolis 项对应反对称结构；
梯度和散度在 $L^2$ 内积下互为负伴随。

半离散系统的空间条件

空间离散后，方程通常写成

M_u\frac{d\bm{u}}{dt}=-R\bm{u}-gGh,

M_h\frac{dh}{dt}=-HD\bm{u}.

其中 $M_u$ 和 $M_h$ 是速度、高度空间的质量矩阵， $R$ 是 Coriolis 项对应的反对称块， $G$ 是高度到速度方程的离散梯度耦合， $D$ 是速度到高度方程的离散散度耦合。

离散能量定义为

E_h=\frac12 H\bm{u}^T M_u \bm{u}+\frac12 g h^T M_h h.

对时间求导并代入半离散方程：

\frac{dE_h}{dt} = -H\bm{u}^T R\bm{u} -gH\bm{u}^T Gh -gH h^T D \bm{u}.

由于反对称矩阵满足

\bm{u}^T R\bm{u}=0,

并且

h^T D\bm{u}=\bm{u}^T D^T h,

所以

\frac{dE_h}{dt} = -gH\bm{u}^T(G+D^T)h.

要让任意 $\bm{u},h$ 下都没有能量产生或耗散，必须有

R^T=-R, \qquad G+D^T=0,

也就是

D=-G^T.

这就是空间离散层面的能量守恒条件。它对应连续层面的分部积分公式

(\nabla h,\bm{u})=-(h,\nabla\cdot \bm{u}).

为什么 DG 高度空间下梯度形式和散度弱形式不等价

在 $H(\mathrm{div})$ 速度空间和 DG 高度空间中，使用散度弱形式

(\partial_t \bm{u},\bm{v})+(f\bm{u}^\perp,\bm{v})-g(h,\nabla\cdot \bm{v})=0

通常更自然，因为它直接把离散负伴随关系写进耦合结构。

连续层面上，

(\nabla h,\bm{v})=-(h,\nabla\cdot \bm{v})

在没有外边界通量时成立，所以

\partial_t \bm{u}+g\nabla h=0

既可以写成

(\partial_t \bm{u},\bm{v})+g(\nabla h,\bm{v})=0,

也可以写成

(\partial_t \bm{u},\bm{v})-g(h,\nabla\cdot \bm{v})=0.

但这个等价隐含了一个前提： $h$ 是全局足够光滑的连续函数， $\nabla h$ 是普通意义下的全局梯度。对 DG 高度空间，这个前提不成立。

在常见的兼容有限元离散里，

\bm{u}_h\in H(\mathrm{div}), \qquad h_h\in DG.

$H(\mathrm{div})$ 速度的法向通量在单元面上连续，但 DG 高度允许在单元面上跳跃。因此 $h_h$ 的分布意义梯度不是单纯的逐单元梯度，而是

\nabla h_h = \nabla_h h_h + \text{internal face jump terms}.

如果只计算

\sum_K \int_K \bm{v}_h\cdot\nabla_h h_h\,dx,

就把内部面上的跳跃贡献丢掉了。

逐单元分部积分给出

\int_K \bm{v}_h\cdot\nabla h_h\,dx = -\int_K h_h\nabla\cdot \bm{v}_h\,dx + \int_{\partial K} h_h\,\bm{v}_h\cdot n_K\,ds.

把所有单元求和后，闭合曲面只能消去外边界项，内部面项仍然保留。设内部面 $F$ 两侧高度为 $h^+$ 和 $h^-$ ，由于 $\bm{v}_h\cdot n$ 连续，内部面贡献可写成

\sum_F \int_F (\bm{v}_h\cdot n)[h_h]\,ds.

DG 高度一般满足 $[h_h]\neq 0$ ，因此这个项不会自动消失。于是

\sum_K \int_K \bm{v}_h\cdot\nabla_h h_h\,dx \ne -\sum_K \int_K h_h\nabla\cdot \bm{v}_h\,dx.

这说明在 $H(\mathrm{div})$ - $DG$ 配对里，逐单元梯度和散度弱形式不是同一个离散算子。

从矩阵角度看，若定义散度耦合

C_{ij} = \int_\Omega (\nabla\cdot v_i)\,h_j\,dx,

那么使用散度弱形式时，速度方程里的高度块和高度方程里的速度块天然满足负转置关系

G=-C, \qquad D=C^T, \qquad G=-D^T.

这正是能量守恒所需的反伴随结构。

但如果把速度方程的高度耦合改成逐单元梯度

G_{ij} = \int_\Omega v_i\cdot\nabla_h h_j\,dx,

而高度方程仍然使用

C_{ij} = \int_\Omega (\nabla\cdot v_i)\,h_j\,dx,

则由于缺少内部面跳跃项，一般有

G\ne -C.

于是负伴随结构被破坏，CN 也就不再能自动保持离散二次能量。

更重要的是，这不只是“少了一点高阶修正”。对 DG 高度来说，原方程中的梯度项本来就应当包含界面跳跃贡献。若只保留逐单元梯度，相当于改写了离散波传播算子本身，因此误差可能很快放大到 $O(1)$ 。常见现象是质量约束仍然看起来稳定，但能量明显漂移，整体误差快速增长。

CN 格式为什么能守恒

空间条件只说明半离散连续时间系统守恒。时间离散以后，还要看一步更新是否继承这个结构。

设半离散系统写成

M_u \frac{d\bm{u}}{dt}+R\bm{u}+gGh=0,

M_h \frac{dh}{dt}+HD\bm{u}=0.

Crank--Nicolson 格式为

M_u(\bm{u}^{n+1}-\bm{u}^n)+\Delta t\left(R\bm{u}^{n+1/2}+gGh^{n+1/2}\right)=0,

M_h(h^{n+1}-h^n)+\Delta t\,HD\bm{u}^{n+1/2}=0,

其中

\bm{u}^{n+1/2}=\frac{\bm{u}^{n+1}+\bm{u}^n}{2}, \qquad h^{n+1/2}=\frac{h^{n+1}+h^n}{2}.

第一式左乘 $H(\bm{u}^{n+1/2})^T$ ，第二式左乘 $g(h^{n+1/2})^T$ 。利用质量矩阵对称性，两式相加得到

E^{n+1}-E^n = -\Delta t\,H(\bm{u}^{n+1/2})^T R\bm{u}^{n+1/2} -\Delta t\,gH(\bm{u}^{n+1/2})^T(G+D^T)h^{n+1/2}.

若空间离散满足

R^T=-R, \qquad D=-G^T,

则

E^{n+1}=E^n.

所以 CN 的作用是把能量差精确写成时间中点处的耦合项；空间反伴随结构再让这个耦合项精确抵消。

为什么 forward Euler 不守恒

即使空间离散满足反伴随关系，forward Euler 也通常不保能量。

把半离散系统抽象成

y_t=Ay.

若 $A^T=-A$ ，连续系统保持二次能量。但 forward Euler 给出

y^{n+1}=(I+\Delta t A)y^n.

于是

\|y^{n+1}\|^2 = (y^n)^T(I+\Delta t A)^T(I+\Delta t A)y^n.

利用 $A^T=-A$ ，

(I+\Delta t A)^T(I+\Delta t A) = (I-\Delta t A)(I+\Delta t A) = I-\Delta t^2A^2.

对反对称矩阵， $-A^2$ 通常是半正定的，因此

\|y^{n+1}\|^2 = \|y^n\|^2+\Delta t^2\|Ay^n\|^2.

除非 $Ay^n=0$ ，否则每一步都会增加能量。空间反伴随只能保证半离散系统守恒，不能让 forward Euler 继承守恒律。

其他保能量时间格式

CN 不是唯一选择。对线性反对称系统，隐式中点格式和 CN 等价，因此同样严格保持二次能量。CN 也可以看成反对称算子的 Cayley 变换；这个更新在能量内积下是正交的。

Leapfrog 或 Störmer--Verlet 常用于波动方程和 Hamiltonian 系统。它们通常不逐步精确保原始能量，但能保持 modified energy，或者让原始能量长期有界振荡。对非线性 Hamiltonian 系统，Average Vector Field 和 Discrete Gradient 方法则通过离散链式法则直接保持 Hamiltonian。