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

汉堡的 BAW 实验室在 1990 年 6 月需要重现 Jade–Weser 河口 13 天的潮位历史。网格用了 15 万个非结构三角形,时间步长要 5 分钟,而且仿真必须比实时跑得更快,才能在排程里塞下下一个算例。显式求解器里,自由表面波速 $c=\sqrt{gh}$ 决定 CFL 上限。28 m 深的区域 $c \approx 16.6$ m/s,最小网格边长接近 $\sqrt{10}$ m,显式步长上限就是 0.2 秒 — 13 天等于 560 万步。本文跟随 Casulli 与 Walters(2000) 提出的 非结构正交网格半隐式算法 突破这堵墙,并且实地跑一个 1D 简化版。结论先说:只把压力项做隐式处理,得到的系统恰好是 SPD,适合 PCG。

论文元信息#

作者: 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 m 时 $c \approx 10$ m/s,深 100 m 时 $c \approx 31$ m/s。流速 $u$ 通常不到 1 m/s,所以物理时间尺度是 $L/u$ ,显式步长却受制于 $\Delta x/c$ 。两者比值 $c/u \approx 10\sim30$ ,显式代码必须把时间切得比流动需要的细 10 到 30 倍。

非结构网格让情况更糟。河口网格里大小不一的单元混在一起,最小那一格决定全局时间步长。Big Lost River 网格在 1.2 km × 1 km 区域内有 12,622 个三角形,最小边接近 5 m,显式上限大约是 0.5 秒。

办法是「慢项显式,快项隐式」。只把自由表面压力梯度做隐式处理,波速就从稳定性条件中消失。

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

普通有限差分套到非结构网格上有两道坎。单元间距随单元形状变化,梯度近似变得不对称。曲线坐标变换又会引入新项,损害精度和稳定性。

Casulli 与 Walters 走的是绕路。正交非结构网格(orthogonal unstructured grid) — 相邻两单元中心连线与它们共享的边相互垂直。Voronoi 多边形与其对偶 Delaunay 三角剖分(在三角形全为锐角时)严格满足这个条件。

正交条件的回报很清楚。在每条边上只定义法向速度 $u_j$ ,自由表面压力梯度 $\partial \eta / \partial n$ 就化为简单差分 $(\eta_{i_2} - \eta_{i_1})/\delta_j$ ,无需曲率修正,离散保持一致。

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

Casulli 的 θ 方法用隐式因子 $\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$ 是单元 $i$ 的面积, $\lambda_j$ 是边 $j$ 的长度, $s_{i,l}$ 是外法向符号。

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

直接联解上述两个方程会得到 $N_z N_s + N_p$ 个未知量的庞大耦合系统。关键技巧是 从动量方程中把 $\mathbf{U}_j^{n+1}$ 解为 $\eta^{n+1}$ 的函数,再代入自由表面方程。 $\mathbf{A}$ 是 SPD,故其逆也是 SPD。

代入后得到:

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(预处理共轭梯度法)正合适,大网格也能在数十次迭代内收敛。

自由表面解出来后,每条边上的动量方程就化为独立的小三对角系统。最后由连续性方程递归算出垂直速度。

Python 实现 1D 半隐式#

线性化的 1D 浅水方程 $\partial_t \eta + H\partial_x u = 0$ 、 $\partial_t u + g\partial_x \eta = 0$ ,可以用同样的思路求解。

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}$ 立刻发散。同样的 $\Delta t$ ,Casulli 跑到底,相位略有滞后,振幅有些衰减(θ 越接近 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 滑块拉到大于 1。红色显式曲线一过 1 就立刻发散到画面之外。青色 Casulli 曲线即使在 CFL=6 也保持形状 — 不过把 θ 调到 1 时振幅衰减得更快(完全隐式的算术平均效果)。θ 接近 0.5 时振幅保存最好,但那是稳定边界,生产代码通常取 0.55–0.6。

批判性考察#

这个算法并非万能。

第一,时间精度。θ=0.5 是二阶,但正好在稳定边界,重力波波前会略有变形。θ ≥ 0.6 时退到一阶。如果定量精度重要,就要把步长缩小,或换更高阶的 IMEX(Implicit-Explicit Runge-Kutta)。

第二,Eulerian–Lagrangian 对流。对流留在显式部分,但单纯的迎风格式带来强数值粘性。所以要配合反向追踪 Lagrange 轨迹再插值的 EL 方法。这种方法在每步跨越多个单元的大 $\Delta t$ 下仍稳定,但精度对插值阶数和轨迹积分误差敏感。

第三,静水压假设。模型采用 $\partial p/\partial z = -\rho g$ 。垂直加速度大的流动(破波、激流瀑布附近)需要非静水压修正。

要点#

仅把自由表面压力项做隐式处理,波速 $\sqrt{gh}$ 就从稳定性条件中消失。时间步长只受流速 $u$ 限制。
关于 $\eta^{n+1}$ 的离散波动方程是 SPD 且对角占优,适合 PCG。动量方程在每条边上分解为独立的小三对角系统。
Voronoi–Delaunay 正交性让非结构网格上的压力梯度离散保持干净 — 这是 Casulli 模型族 30 年来不断扩展的基础。