[论文评论] 在风暴来临前算完仿真 — Casulli–Walters(2000) 非结构正交网格的浅水方程

1990年6月，汉堡 BAW 得把 Jade–Weser 河口13天的潮位重新对上。三角形15万个。时间步长5分钟。必须比实时更快结束，下一次预报才有计算时间。

显式推进时，自由面波速 $c=\sqrt{gh}$ 抓住 CFL。水深28米， $c \approx 16.6$ m/s。最短边接近 $\sqrt{10}$ 米，允许步长就是0.2秒。13天是1,123,200秒。用0.2秒来切是5,616,000步，和论文里的560万步一致。5分钟则是3,744步。比值正好是1500倍。因为300秒除以0.2秒，换哪一段日期，这个倍数都不变。

撑住这个倍数的是 Casulli–Walters(2000) 的半隐式。只有自由面梯度放在新时刻，对流留在旧时刻。

把 CFL 滑块拉过1，红色显式曲线会炸开，而把 θ 调高只会削掉青色曲线的振幅。

CFL ratio (Δt / Δt_explicit): 2.00θ (implicitness): 0.60

CFL > 1 ⇒ explicit blows up; semi-implicit stays bounded for any CFL when θ ≥ 1/2.

CFL=6 时青色的峰还在。把 θ 推到1，峰很快死掉。振幅还站得住的实务区间是0.55到0.6。

论文元信息

作者: Vincenzo Casulli (Univ. of Trento), Roy A. Walters (USGS)
期刊: International Journal for Numerical Methods in Fluids, 32, 331–348 (2000)
题目: "An unstructured grid, three-dimensional model based on the shallow water equations"
关键词: semi-implicit, shallow water, unstructured orthogonal grid, wetting/drying

显式求解器为何在河口失灵

水深10米时 $c \approx 10$ m/s，100米时 $c \approx 31$ m/s。流速 $u$ 通常不到1 m/s。潮水走过一段河的时间是 $L/u$ ，显式步长却卡在 $\Delta x/c$ 。 $c/u \approx 10\sim30$ 时，水过一次，求解器把时间切得比这更碎。

非结构网格把这个限制全部推到最短边上。Big Lost River 在1.2 km × 1 km 里有12,622个三角形，最短边约5米，显式上限大约0.5秒。5米除以0.5秒是10 m/s。停住时钟的不是深槽里的16 m/s，而是浅处那条短边。

慢的项显式，快的项隐式。只把自由面压力梯度做成隐式， $\sqrt{gh}$ 就退出稳定条件。剩下的限制在流速那边。

正交非结构网格 — Voronoi–Delaunay 的对偶#

有限差分直接放到非结构网格上，有两处会坏。单元间距随形状变，梯度就歪向一侧。改走曲线坐标，变换项会把精度和稳定性一起削掉。

Casulli 与 Walters 限制的是网格。正交非结构网格：相邻中心的连线垂直于共享边。Voronoi 多边形和它的对偶 Delaunay 三角剖分，在全是锐角三角形时满足这个条件。混进一个钝角，外心就跑出三角形。中心连线即使仍垂直于边，也不再穿过单元内部。让外接圆保持为空的插点能把这个锐角维持多久，这个差分就成立多久。

正交换来的是一行差分。边上只放法向速度 $u_j$ ， $\partial \eta / \partial n$ 就是 $(\eta_{i_2} - \eta_{i_1})/\delta_j$ 。没有曲率修正。

半隐式分裂: θ 方法分离两个时间尺度

θ 方法用隐式因子 $\theta \in [1/2, 1]$ 把时间切开。动量里自由面梯度写成 $\theta \eta^{n+1} + (1-\theta)\eta^n$ ，自由面方程里速度写成 $\theta u^{n+1} + (1-\theta)u^n$ 。

离散动量如下。垂直粘性和风应力都在里面。

\mathbf{A}_j^n \mathbf{U}_j^{n+1} = \mathbf{G}_j^n - \theta g \frac{\Delta t}{\delta_j}\!\left[\eta_{i(j,2)}^{n+1} - \eta_{i(j,1)}^{n+1}\right]\Delta\mathbf{Z}_j^n

$\mathbf{U}_j^{n+1}$ 是边 $j$ 上的法向速度， $\mathbf{A}_j^n$ 是垂直粘性、风和底摩擦合成的三对角， $\mathbf{G}_j^n$ 是对流、Coriolis 和水平摩擦。每条边各自独立，于是分成 $N_s$ 个 $N_z \times N_z$ 的三对角。

自由面方程是这样。

P_i \eta_i^{n+1} = P_i \eta_i^n - \theta \Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n+1} - (1-\theta)\Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n}

$P_i$ 是单元面积， $\lambda_j$ 是边长， $s_{i,l}$ 是朝外的符号。

化为波动方程,SPD 的回报#

两式原样联立，未知量有 $N_z N_s + N_p$ 个。从动量里把 $\mathbf{U}_j^{n+1}$ 解成 $\eta^{n+1}$ 的函数，再代进自由面。 $\mathbf{A}$ 正定，逆也正定。

代进去，掉出来的是离散波动方程。

P_i \eta_i^{n+1} - g \Delta t^2 \theta^2 \sum_{l=1}^{S_i} \frac{s_{i,l}\lambda_{j(i,l)}}{\delta_{j(i,l)}}\!\left[(\Delta\mathbf{Z})^\top \mathbf{A}^{-1}\Delta\mathbf{Z}\right]_{j(i,l)}^{n}\!\!\left(\eta^{n+1}_{i[j(i,l),2]} - \eta^{n+1}_{i[j(i,l),1]}\right) = R_i

未知量只剩 $\eta^{n+1}$ 。规模是单元数 $N_p$ 。矩阵对称，对角强，正定。PCG 合适。论文写的是，大网格上迭代也在几十步内停下。

没有垂直结构的一维里， $\mathbf{A}$ 就是1，代码里的 $k = g H \Delta t^2 \theta^2 / \Delta x^2$ 就是邻居系数。 $\Delta t$ 取显式上限的4倍，且 $c^2 = gH$ ，则 $k = 16\theta^2$ 。θ=0.6 时 $k=5.76$ 。内部对角 $1+2k$ 是12.52，两条非对角之和是11.52。对角多出来的宽度永远是1，这1来自储存项 $P_i \eta$ 。步长加大时， $k$ 同时贴在对角和非对角上，盈余并不增加。前面滑块上的 CFL=6、θ=0.6，给出 $k=12.96$ ，对角26.92，非对角之和25.92。宽度仍是1。再隐一些，矩阵也不会更安全。死掉的是振幅。

$\eta$ 解出以后， $\mathbf{U}_j^{n+1}$ 回到每条边上的小三对角。垂直速度沿不可压关系往回补。

Python 实现 1D 半隐式#

同一个分裂。 $N=200$ ， $L=1$ ， $H=1$ ， $g=9.81$ ，θ=0.6。

import numpy as np
from scipy.linalg import solve_banded
 
class CasulliShallowWater1D:
    """1D 线性浅水 — 半隐式 θ 方法。"""
 
    def __init__(self, N=200, L=1.0, H=1.0, g=9.81, theta=0.6):
        self.N, self.L, self.H, self.g, self.theta = N, L, H, g, theta
        self.dx = L / N
        self.eta = np.zeros(N)            # 单元中心的自由表面
        self.u = np.zeros(N + 1)          # 边(face)上的法向速度,两端 = 0
 
    def initial_bump(self, x0=0.3, sigma=0.05, amp=0.05):
        x = (np.arange(self.N) + 0.5) * self.dx
        self.eta = amp * np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
        self.u[:] = 0.0
 
    def step(self, dt):
        N, dx, g, H, th = self.N, self.dx, self.g, self.H, self.theta
        k = g * H * dt * dt * th * th / (dx * dx)   # 波动项系数
        alpha = H * dt / dx
 
        # G_j: 收集到边上的显式部分(两端 = 0)
        G = np.zeros(N + 1)
        G[1:N] = self.u[1:N] - g * (dt / dx) * (1 - th) * (self.eta[1:] - self.eta[:-1])
 
        # 三对角: -k η_{i-1} + (1+2k) η_i - k η_{i+1} = R_i
        ab = np.zeros((3, N))
        ab[0, 1:] = -k                         # 上对角
        ab[2, :-1] = -k                        # 下对角
        ab[1, :] = 1.0 + 2 * k
        ab[1, 0] = 1.0 + k                     # 左壁
        ab[1, -1] = 1.0 + k                    # 右壁
 
        rhs = self.eta.copy()
        rhs -= alpha * th * (G[1:] - G[:-1])
        rhs -= alpha * (1 - th) * (self.u[1:] - self.u[:-1])
 
        new_eta = solve_banded((1, 1), ab, rhs)
 
        new_u = self.u.copy()
        new_u[1:N] = self.u[1:N] - g * (dt / dx) * (
            th * (new_eta[1:] - new_eta[:-1]) + (1 - th) * (self.eta[1:] - self.eta[:-1])
        )
        new_u[0] = new_u[-1] = 0.0
        self.eta, self.u = new_eta, new_u
        return self.eta
 
if __name__ == "__main__":
    sim = CasulliShallowWater1D(N=200, theta=0.6)
    sim.initial_bump()
    c = np.sqrt(sim.g * sim.H)
    dt_explicit_limit = sim.dx / c
    dt = 4.0 * dt_explicit_limit                # 显式上限的 4 倍
    for _ in range(int(0.6 / dt)):
        sim.step(dt)
    print(f"max|η| = {np.max(np.abs(sim.eta)):.4f}  (仍然有界)")

显式在 $\Delta t = 4 \Delta t_\text{ex}$ 上很快发散。Casulli 用同一间隔走完长度0.6的时段。相位滞后。θ 越靠近1，振幅死得越快。θ=0.5 时权重各半，是二阶。超过0.6就变成一阶。代码默认的0.6，已经站在这一阶一侧。

调大时间步长试试

渠道长1、单元200，则 $\Delta x = 0.005$ 。 $c=\sqrt{9.81}\approx 3.13$ m/s，显式上限约0.0016秒，四倍约0.0064秒。0.6除以这个间隔再丢掉小数，剩下93步。输出只有 max|η| 有限这一句，初始隆起0.05衰减了多少，这一句并不说。

这0.05是水深 $H=1$ 的5%。线性化 $\partial_t \eta + H\partial_x u = 0$ 要在 $\eta$ 远小于 $H$ 时才闭合。于是波速只有 $\sqrt{gH}$ 这一个数，三对角的 $k$ 在每个单元上也相同。把隆起加到水深量级， $c$ 就随位置变，这个常数 $k$ 不再能把波动方程闭上。

θ=0.5 配上同一个四倍步长， $k=16\times 0.25=4$ 。对角盈余仍是1。所谓稳定边界不来自矩阵，而来自时间权重 $(θ-1/2)$ 。它为零时数值粘性关掉，网格模态来回走，不衰减。

批判性考察

不是万能算法。

时间阶数先就不一样。θ=0.5 是二阶，可它站在稳定边界上，重力波波前容易弯。过了0.6就是一阶。验潮记录的振幅若是交付物，就缩小步长，或换更高阶的 IMEX。

对流留在显式项里。单纯迎风的数值粘性很黏。所以才配上沿 Lagrange 轨迹倒着走再插值的 EL。 $\Delta t$ 跨过不止一个单元也还能稳住，但结果会随插值阶数和出发点误差晃。

静水压 $\partial p/\partial z = -\rho g$ 也是前提。破波和跌水这类垂向加速度大的地方，需要非静水压修正。

混进一个钝角，五分钟步长就作废

5分钟配 θ=0.6，换来的是步数变成原来的1500分之一，付出的是波高的一阶衰减。验潮振幅若是交付物，就不要把 θ 往1那边推。对角多出来的宽度已经是1，再隐，矩阵也不会更安全。死的是波。

这个差分等于 $\partial \eta / \partial n$ ，仅当中心连线垂直于面。一个钝角把外心推出单元，等式就结束。面上会留下不垂直的分量。这个分量能不能按扩散通量的方式补回去，正是非正交面上的扩散修正在问的事。