SWE 的原理解析

二、SWE 的运动机制

SWE 的运动机制可以理解成：水面坡度改变流速，流动的汇聚和发散改变水深。

假设你在平静水面中央制造一个隆起：

所以，SWE 需要同时保存“水有多少”和“水如何运动”。

这也解释了为什么：仅根据当前高度图，通常不能唯一决定下一时刻的水面。 两片形状完全相同的水面，可能一个正在向外流动，另一个正在向内汇聚，它们下一刻的变化会不同。

1、质量守恒

对于密度恒定的水，质量守恒等价于体积守恒。

不可压缩流体的三维连续性方程（散度形式）：

$$\frac{\partial u}{\partial x}+\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}=0$$

考虑水平面积为 A 的一个网格单元，内部水的体积为：

$$V=Ah$$

若流入、流出的体积流量分别为 Qin、Qout，那么：

$$\frac{dV}{dt}=Q_{in}-Q_{out}$$

因此：

$$\frac{dh}{dt}=\frac{Q_{in}-Q_{out}}{A}$$

例如，面积为1m^2的单元，每秒流入 0.03m^3、流出 0.01m^3，水深就每秒增加 0.02m。

把这个关系写成连续空间中的方程：

$$\frac{\partial h}{\partial t}+\frac{\partial (hu)}{\partial x}+\frac{\partial (hv)}{\partial z}=0$$

其中：

• hu：沿 x 方向的单位宽度体积流量；
• hv：沿 z 方向的单位宽度体积流量；
• 后两项之和：一个位置附近净流出的程度。

净流出为正，水深下降；净流入为正，水深上升。这就是没有外部水源时的连续性方程。

展开可以得到：

$$\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial x}+v\frac{\partial h}{\partial z}=-h\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial z}\right)$$

右侧括号内是水平速度的散度：

$$\nabla \cdot \mathbf{u}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial z}$$

跟随水流观察时：

• 速度场向外发散，水柱被摊开，水深减小；
• 速度场向内汇聚，水被堆积，水深增大。

水本身没有被压缩，改变的是水柱的形状。

2、水平动量方程

另外两条核心方程是水平动量方程，它们决定流速怎样变化。

先考虑固定底床，并忽略摩擦、风力和地球自转：

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial z}=-g\frac{\partial \eta }{\partial x}$$

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial z}=-g\frac{\partial \eta }{\partial x}$$

这两条式子结构相同，只是分别描述两个水平方向。

$$\underset{\text{固定位置的流速变化}}{\frac{\partial u}{\partial t}}+\underset{\text{流动带来的对流加速度}}{u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial z}}=\underset{\text{水面坡度引起的加速度}}{-g\frac{\partial \eta }{\partial x}}$$

“对流加速度”可以这样理解：水流到了另一个速度不同的位置，跟随这团水观察，它的速度就会发生变化。

右边的负号表示，重力所产生的压力作用会使水沿水面降低的方向加速。水面坡度决定的是加速度，而不是直接决定速度。 因此，水面刚变平时，水仍然可能继续流动。

压力项的来源也很直观。静水压力近似为：

$$p=\rho g(\eta -Y)$$

这里使用相对于大气压的压力。在相同竖直高度 Y 上，水面越高，上方压着的水越多，压力越大：

$$-\frac{1}{\rho }\frac{\partial p}{\partial x}=-g\frac{\partial \eta }{\partial x}$$

这就得到了动量方程中的驱动力。

2.1 动量

动量（Momentum）：是流体带有的“运动的力量”。

$$p=mv$$

2.2 动量通量

通量（Flux）：​单位时间内，垂直穿过单位面积的某种物理量​。

动量通量（Momentum Flux）就是指：单位时间内，通过单位面积的动量。

$$动量通量=\rho u\otimes u$$

其中 ρ 是密度，u 是速度矢量，⊗ 是张量积。

如果把它展开，在二维（x, z 平面）的浅水方程中，动量通量张量大概长这样：

$$\left[\begin{array}{cc}\rho u^2& \rho uw\\ \rho wu& \rho w^2\end{array}\right]$$

• 对角线上的 ρu² 和 ρw² 叫做​正应力​（相当于水流动带来的动压）。
• 非对角线上的 ρuw 叫做​剪应力​（不同方向水流交错摩擦、拉扯带来的力）。

2.3 三维 N-S 方程

对于不可压缩、无粘性的流体，三维 N-S 方程在 xx 方向的动量方程为：

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+f\mathrm{x}$$

其中 u, v, w 分别是 x, y, z 方向的速度， y 是竖直方向，p 是压强, fx 是外力。

​2.4 x ​方向动量守恒方程。

2.4.1 从三维N-S方程出发(x方向)

利用连续性方程（质量守恒的散度形式），我们可以将左边的对流项（非线性项）改写为​守恒形式​（Flux form，把所有东西都塞进偏导数里面的形式）：

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$$

代入 N-S 方程：

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}+u(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z})=-\frac{1}{\rho }\frac{\partial p}{\partial x}+f\mathrm{x}$$

$$\frac{\partial u}{\partial t}+2u\frac{\partial u}{\partial x}+(v\frac{\partial u}{\partial y}+u\frac{\partial v}{\partial y})+(w\frac{\partial u}{\partial z}+u\frac{\partial w}{\partial z})=-\frac{1}{\rho }\frac{\partial p}{\partial x}+f\mathrm{x}$$

整理得：

$$\frac{\partial u}{\partial t}+\frac{\partial u^2}{\partial x}+\frac{\partial uv}{\partial y}+\frac{\partial uw}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+f\mathrm{x}$$

2.4.2 沿水深 y 方向积分

由于浅水方程得核心假设是：水平尺度远大于垂直尺度，那么我们假设垂直方向的加速度可以忽略不记，那么有：

$$v\frac{\partial u}{\partial y}\approx u\frac{\partial v}{\partial y}\approx 0$$

将上述方程从水底 b(x, z)积分到水面 η(x,z,t)：

$${\int }_b^{\eta }\frac{\partial u}{\partial t}dy=\frac{\partial }{\partial t}{\int }_b^{\eta }udy=\frac{\partial hu}{\partial t}$$

$${\int }_b^{\eta }\frac{\partial (u^2)}{\partial x}dy=\frac{\partial }{\partial x}{\int }_b^{\eta }{u}^{2}dy=\frac{\partial (h{u}^2)}{\partial x}$$

$${\int }_b^{\eta }\frac{\partial (uw)}{\partial z}dy=\frac{\partial }{\partial z}{\int }_b^{\eta }uwdy=\frac{\partial (huw)}{\partial z}$$

将 w 改为 v，表示 z 方向上的水流速度。

2.4.3 处理压强梯度项

根据​静水压强假设​，任意深度 yy 处的压强由上方水柱的重量产生：

$$p(y)=\rho g(\eta -y)$$

(其中 ​η ​是水面高度，ρ ​是密度，g ​是重力加速度)

将其代入积分（利用莱布尼茨积分法则，并考虑到水面 ηη 和水底 bb 随 xx 的变化）：

$$-\frac{1}{\rho }{\int }_b^{\eta }\frac{\partial p}{\partial x}dy=-g{\int }_b^{\eta }\frac{\partial (\eta -y)}{\partial x}dy=-gh\frac{\partial \eta }{\partial x}$$

2.4.4 处理底坡/外力项

方程右边还有一个重力在底坡上的分量 fx，重力沿坡面的分力为：

$$f\mathrm{x}=-g\frac{\partial b}{\partial x}$$

积分后（考虑到 hh 的厚度）：

$${\int }_b^{\eta }f\mathrm{x}dy=-gh\frac{\partial b}{\partial x}$$

2.4.5 合并与化简

最终得到：

$$\frac{\partial (hu)}{\partial t}+\frac{\partial }{\partial x}(h{u}^2+\frac{1}{2}g{h}^2)+\frac{\partial (huv)}{\partial z}=-gh\frac{\partial b}{\partial x}$$

如上我们就得到了包含固定底床起伏的二维非线性 SWE 方程组：

$$①\frac{\partial h}{\partial t}+\frac{\partial (hu)}{\partial x}+\frac{\partial (hv)}{\partial z}=0$$

$$②\frac{\partial (hu)}{\partial t}+\frac{\partial }{\partial x}(h{u}^2+\frac{1}{2}g{h}^2)+\frac{\partial (huv)}{\partial z}=-gh\frac{\partial b}{\partial x}$$

$$③\frac{\partial (hv)}{\partial t}+\frac{\partial (huv)}{\partial x}+\frac{\partial }{\partial z}(h{u}^2+\frac{1}{2}g{h}^2)=-gh\frac{\partial b}{\partial z}$$

3、初始化项

有地形时，一个关键性质是：水深不同，也可以完全静止。
静水平衡满足：

$$u=v=0,\eta =常数$$

由于：

$$h=\eta -b$$

底床起伏会使水深变化，但水平水面不会因此产生水平加速度。

在守恒方程中，这种静止依靠两部分相互抵消：

$$\frac{\partial }{\partial x}\left(\frac{1}{2}g{h}^2\right)=-gh\frac{\partial b}{\partial x}$$

如果程序分别近似压力通量和底床源项，却没有让它们正确平衡，就可能出现“没有施加扰动，水也自己流动”的假象。

能够在离散计算中保持这类平衡的数值方法称为 ​well-balanced，良平衡方法​。静水重构是实现这一性质的一类经典方法。

因此，做 Unity 地形水体时，静水初始化通常应当是：

$$h(x,z,0)=\max (\eta \_0-b(x,z),0)$$

4、SWE 波速

4.1 连续性方程

考虑平底、静止背景水体，平均水深为常数 H。令：

$$h=H+\zeta$$

为了推导方便，我们现在只看一维（x 方向），忽略 z 方向的流动（或者假设流动在 z 方向均匀），方程变成：

$$\frac{\partial h}{\partial t}+\frac{\partial (hu)}{\partial x}=0$$

利用乘积求导法则展开第二项：

$$\frac{\partial h}{\partial t}+h\frac{\partial u}{\partial x}+u\frac{\partial h}{\partial x}=0$$

引入“微小扰动（小振幅波）”假设，假设水面本来平静（静止水深为 H），你往水里丢了一颗小石子。水面产生了一个非常微小的波澜（高度为 ζ）。

此时：总水深 h(x,t)=H+ζ(x,t)，并且水流速度 u(x,t) 也非常小。

把 h=H+ζ 代入刚才的一维方程：

$$\frac{\partial (H+\zeta )}{\partial t}+(H+\zeta )\frac{\partial u}{\partial x}+u\frac{\partial (H+\zeta )}{\partial x}=0$$

由于H 是不随时间变化的常数，则有：

$$\frac{\partial H}{\partial t}=\frac{\partial H}{\partial x}=0$$

又因为 ζ 和 u 都是微小量，可以忽略，也有：

$$\zeta \approx 0,u\approx 0$$

所以原式等于：

$$\frac{\partial \zeta }{\partial t}+H\frac{\partial u}{\partial x}=0$$

4.2 动量方程

在只考虑一维（x方向）且忽略底部摩擦的情况下，浅水方程的标准动量方程是：

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=-g\frac{\partial \eta }{\partial x}$$

• ​左边​：速度随时间的变化率 + ​对流加速度​（你之前问过的，流体跑到了流速不同的地方产生的加速度）。
• ​右边​：由水面高度差（压力梯度）产生的加速度（水往低处流）。

现在，假设水面非常平静，轻轻吹了一口气，产生了一个极小的涟漪。

• 水流速度 u 非常非常小（趋近于 0）。
• 水面高度的变化量 ζ 也非常非常小。

则有：

$$u\frac{\partial u}{\partial x}\approx 0$$

那么原式就变为：

$$\frac{\partial u}{\partial t}+g\frac{\partial \eta }{\partial x}=0$$

4.3 计算波速

对连续性方程两边求 t 的偏导：

$$\frac{{\partial }^2\zeta }{\partial t^2}+H\frac{{\partial }^{2}u}{\partial x\partial t}=0$$

对动量方程两边求 x 的偏导：

$$\frac{{\partial }^{2}u}{\partial x\partial t}+g\frac{\partial \eta }{{\partial }^{2}x}=0$$

联立得到：

$$\frac{{\partial }^2\zeta }{\partial t^2}=Hg\frac{{\partial }^2\eta }{\partial x^2}$$

又因为：

$$\eta =\zeta +H$$

则有：

$$\frac{{\partial }^2\eta }{\partial t^2}=Hg\frac{{\partial }^2\eta }{\partial x^2}$$

波动方程标准形式为：

$$\frac{{\partial }^2\eta }{\partial t^2}=c^2\frac{{\partial }^2\eta }{\partial x^2}$$

所以得到物理波速公式：

$$c=\sqrt{gH}$$

• g 是重力加速度（约 9.81 m/s29.81m/s2）
• HH 是​静水深​（或者当地实际水深 hh）

​这个公式非常直观​：

• 水越深，波传得越快（太平洋里的海啸速度接近飞机）。
• 水越浅，波传得越慢（冲上海滩的波浪会减速）。

三、数值离散方法

从方程到程序，需要选择数值离散方法。SWE 本身并不是某一种固定算法。

常见的层次关系是：

3.1 物理模型

3.1.1 选择物理模型的目的

真实水体包含大量现象：

• 水平流动、竖直运动；
• 压力变化；
• 黏性、摩擦和湍流；
• 波浪传播、破碎、飞溅；
• 不同波长的波以不同速度传播。

如果完整保留这些现象，求解成本会非常高。因此，需要选择一个适合研究目标的简化模型。

SWE 的主要简化是：沿水深方向进行平均，用水深和深度平均水平速度描述水体，并采用静水压力近似。它适用于水深相对水平运动尺度较小的情况。

选择物理模型，本质上是在决定：哪些现象需要计算，哪些现象可以近似或忽略。

3.1.2 一维非线性 SWE

对于平底、没有摩擦和外力的一维水流，非线性 SWE 为：

$$\frac{\partial h}{\partial t}+\frac{\partial (hu)}{\partial x}=0$$

$$\frac{\partial (hu)}{\partial x}+\frac{\partial }{\partial x}(h{u}^2+\frac{1}{2}g{h}^2)=0$$

第一条描述质量守恒，第二条描述水平动量守恒。

定义状态向量：

$$U=\left[\begin{array}{c}h\\ m\end{array}\right],m=hu$$

其中，m 是省略常数密度因子的动量变量。

再定义物理通量：

$$F(U)=\left[\begin{array}{c}m\\ \frac{m^2}{h}+\frac{1}{2}g{h}^2\end{array}\right]$$

于是，两条方程可以统一写成：

$$\frac{\partial U}{\partial t}+\frac{\partial F(U)}{\partial x}=0$$

这里的“非线性”，意味着状态量之间存在乘积、平方等关系。例如：

$$hu，h{u}^2，h^2$$

水深和速度会相互影响，波的传播速度也依赖当前状态：

$$u-\sqrt{gh}，u+\sqrt{gh}$$

因此，波浪的不同部分可能具有不同传播速度，逐渐形成陡峭波面，甚至出现水跃等间断结构。

3.1.3 一维线性 SWEH

假设水体最初静止，平均水深为常数 H，扰动为 ζ：

$$h=H+\zeta$$

线性化时忽略扰动量之间的高阶乘积，化简连续性方程得到：

$$\frac{\partial \zeta }{\partial t}+H\frac{\partial u}{\partial x}=0$$

动量方程相应简化为：

$$\frac{\partial u}{\partial t}+g\frac{\partial \zeta }{\partial x}=0$$

它们描述围绕静止背景的小扰动传播，波速为：

$$c=\sqrt{gH}$$

3.1.4 两种物理模型的区别

把网格加密，不能补回物理模型已经忽略的现象。

例如，线性 SWE 即使用很细的网格，也不会自动变成非线性 SWE；基础 SWE 即使求解得非常精确，也不会自动产生真实水滴和翻卷浪头。

此外，物理问题还需要配套的初始条件、边界条件和外部作用。例如：最初的水深、墙壁位置、入口流量，以及物体怎样向水体施加扰动。

3.2 空间离散

3.2.1 空间离散的含义

连续模型中的：

$$h(x,t)，u(x,t)$$

在任意位置都具有值。计算机无法保存无限多个位置，因此需要用有限个数据近似整个空间。

例如，将一条水槽划分为 N 个区域：

$$\Delta x=\frac{L}{N}$$

但“划分网格”只是开始，还需要明确：

• 每个位置保存什么；
• 一个区域内部的状态怎样表示；
• 空间导数或区域积分怎样近似。

有限差分、有限体积和有限元，在这些问题上采用不同思路。下面通过一维网格来介绍这些方法。

3.2.2 有限差分法

有限差分法通常保存若干位置上的函数值：

$$h_i(t)\approx h(x_i,t)$$

例如，用中心差分近似水面坡度：

∂η∂x∣xi≈ηi+1−ηi−12Δx
\\left. \\frac{\\partial \\eta}{\\partial x} \\right |_{x_i} \\approx \\frac{\\eta_{i+1}-\\eta_{i-1}}{2\\Delta x}
∂x∂η​​xi​​≈2Δxηi+1​−ηi−1​​

如果采用线性 SWE：

$$\frac{\partial u}{\partial t}=-g\frac{\partial \zeta }{\partial x}$$

空间离散后就得到：

$$\frac{d{u}_i}{dt}\approx -g\frac{{\zeta }_{i+1}-{\zeta }_{i-1}}{2\Delta x}$$

原本连续的空间导数，变成了相邻数据之间的运算。注意，此时只完成了空间离散，还没有决定如何推进时间。

有限差分法的直观特点是：

• 对规则网格很容易实现；
• 相邻采样和纹理操作比较自然；
• 可以构造高阶导数近似；
• 稳定性和守恒性取决于具体格式。

空间离散（算空间导数）和时间离散（算时间推进）必须满足稳定性条件（CFL 条件），且不同组合的稳定性完全不同。 不要以为随便用什么方法（比如中心差分 + 显示时间推进）都能算出正确结果。

3.2.3 有限体积法

有限体积法把每个网格视为一个有大小的控制体。

1、有限体积法区域内存量变化计算

一维情况下，单元 i 的左右边界为：

$$x_{i-\frac{1}{2}}，x_{i+\frac{1}{2}}$$

这里：

• 整数下标 i：单元；
• 半整数下标 i+1/2：两个单元之间的界面。

保存的状态是单元平均值：

$$\overline{U_i}=\frac{1}{\Delta x}{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}U(x,t)dx$$

公式两边同时对时间 t 求导数：

Ui‾=ddt(1Δx∫xi−12xi+12U(x,t) dx)
\\overline{U_i} = \\frac{d}{dt} \\left ( \\frac{1}{\\Delta x}\\int_{x_{i-\\frac{1}{2}}}^{x_{i+\\frac{1}{2}}} U(x, t)\\,dx \\right )
Ui​​=dtd​​Δx1​∫xi−21​​xi+21​​​U(x,t)dx​

积分区间固定，可以由莱布尼茨积分法则把时间导数直接移到积分括号里面：

$$\overline{U_i}=\frac{1}{\Delta x}{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\frac{\partial U(x,t)}{\partial t}dx$$

由一维浅水方程核心部分：

$$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}=0$$

$$其中U=[h,hu{]}^T，F=[hu,h{u}^2+\frac{1}{2}g{h}^2{]}^T$$

将其带入刚才的积分中：

$$\frac{d\overline{U_i}}{dt}=\frac{1}{\Delta x}{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left(-\frac{\partial F(x,t)}{\partial x}\right)dx$$

最终得到：

$$\frac{d\overline{U_i}}{dt}=-\frac{F(x_{i+\frac{1}{2}})-F(x_{i-\frac{1}{2}})}{\Delta x}$$

实际计算时，再用数值通量近似界面通量：

$$\frac{d\overline{U_i}}{dt}=-\frac{F(x_{i+\frac{1}{2}})-F(x_{i-\frac{1}{2}})}{\Delta x}$$

上式表达的是：单元内部存量的变化。由左右两侧的净流入决定。

2、有限体积元的内部守恒变量的处理

① 常数

有限体积法并不意味着每个单元内部永远只能是一个常数。最简单的做法是分片常数，即整个单元都使用同一个状态。这通常对应基础的一阶空间近似。即：

$$U(x)\approx {\overline{U}}_i$$

② 线性重构

使用分片线性重构：

$$U_i(x)\approx {\overline{U}}_i+{\sigma }_i(x-x_i)$$

其中 σi 是单元内的斜率。

这样可以更好地表示平滑波形。但在陡峭波面附近，直接使用较大的斜率可能产生过冲，例如计算出不合理的负水深。因此，通常需要限制器；对于 SWE，还可能需要额外的正水深约束与静水平衡处理。

3.2.4 有限元法

有限元法的基本形式是：

$$U_h(x,t)=\sum _j{a}_j(t){\varphi }_j(x)$$

其中：

• ϕj(x)ϕj​(x)：预先选定的基函数；
• aj(t)aj​(t)：需要求解的系数；
• UhUh​：对真实状态的近似。

例如，最简单的线性单元，可以通过端点值在单元内部进行线性插值。

随后，将微分方程乘以测试函数并进行积分，使方程在这种积分意义下成立。这称为​弱形式​。

有限元也有不同类型。连续 Galerkin 有限元通常要求相邻单元的解满足相应连续性；间断 Galerkin（DG）允许单元之间出现状态跳跃，因此同样需要数值界面通量来连接相邻单元。

3.2.5 三种离散方法的区别

3.3 界面通量

3.3.1 物理通量和数值通量的关系

物理模型给出了：

$$F(U)=\left[\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\end{array}\right]$$

但在离散网格中，同一个界面有两个来源不同的状态：

$$U_L，U_R$$

它们分别来自：

• 左单元的右边界；
• 右单元的左边界。

两侧水深、速度可能不同，因此：

$$F(U_L)\ne F(U_R)$$

而一个共享界面需要有统一的交换量，因此要根据左右两侧状态计算：

$$\overline{F}=\mathcal{F}(U_L,U_R)$$

这里：

• F：由物理模型规定的物理通量；
• \\mathcal F：选择的数值通量算法；
• \\overline F：实际用于网格更新的界面通量。

3.3.2 局部 Riemann 问题

可以把两个单元之间的界面暂时看成一道隔板：

• 左侧具有状态 UL；
• 右侧具有状态 UR；
• 在某一时刻移除隔板。

随后会产生波，并在界面附近形成新的流动状态。

这种左右初始状态各自为常数的问题，称为 ​Riemann 问题​。

精确求解它可以得到界面处的状态和通量。但在大量网格、每个时间步都要计算的情况下，通常会采用近似求解器。

3.3.3 Rusanov：平均物理通量，加上数值耗散

Rusanov 通量也称局部 Lax-Friderichs 通量：

$$F_{Rus}=\frac{F_L+F_R}{2}-\frac{a}{2}(U_R-U_L)$$

其中：

$$F_L=F(U_L)，F_R=F(U_R)$$

对全湿 SWE，一种常见的局部速度估计为：

$$a=max\left(|u_L|+\sqrt{g{h}_L}，|u_R|+\sqrt{g{h}_R}\right)$$

公式分为两部分：

$$\frac{F_L+F_R}{2}$$

取两侧物理通量的平均值；

$$-\frac{a}{2}(U_R-U_L)$$

根据状态差异加入数值耗散。

这种耗散有助于抑制振荡，但也可能把陡峭变化和细小波纹抹平。因此，Rusanov 实现简单，但计算出的波浪细节也与网格大小、重构阶数有关。

例如，设：

$$\begin{gathered} h_L=2，h_R=1 \\ {u}_L=u_R=0，g=9.81 \end{gathered}$$

则：

$$U_L=\left[\begin{array}{c}2\\ 0\end{array}\right]，U_R=\left[\begin{array}{c}1\\ 0\end{array}\right]$$

物理通量为：

$$F_L=\left[\begin{array}{c}0\\ 19.62\end{array}\right]，F_R=\left[\begin{array}{c}0\\ 4.905\end{array}\right]$$

最大速度估计：

$$a=\sqrt{9.81\times 2}\approx 4.429$$

最终：

$$F_{Rus}\approx \left[\begin{array}{c}2.215\\ 12.263\end{array}\right]$$

第一项是单位宽度体积流量，正号表示向右输水。

尽管两侧初始速度都是零，水深差引起的压力不平衡仍会产生流动。 数值通量是在近似界面问题随即发生的演化。

3.3.4 HLL：分别估计向左和向右的波速

HLL 使用两个波速：

$$S_L，S_R$$

分别表示最左侧和最右侧传播信号的速度估计。

对全湿 SWE，一种基础估计是：

$$\begin{gathered} S_L=min\left(u_L-\sqrt{g{h}_L}，u_R-\sqrt{g{h}_R}\right) \\ {S}_L=max\left(u_L+\sqrt{g{h}_L}，u_R+\sqrt{g{h}_R}\right) \end{gathered}$$

具体求解器也可能采用更稳健的速度估计，干床情况尤其需要专门处理。

HLL 通量为：

$$F_{HLL}=\{\begin{array}{ll}{F}_L,& S_L\ge 0\\ \frac{S_R{F}_L-S_L{F}_R+S_L{S}_R(U_R-U_L)}{S_R-S_L},& S_L<0<S_R\\ F_R,& S_R\le 0\end{array}$$

它体现了传播方向：

• 如果所有信号都向右传播，界面由左侧状态决定；
• 如果所有信号都向左传播，界面由右侧状态决定；
• 如果信号向两侧传播，就使用中间表达式。

HLL 将两个外侧波之间的复杂结构近似为一个平均状态。

Rusanov 与 HLL 还有一个直接关系。令：

$$S_L=-a，S_R=a$$

HLL 就会化简成 Rusanov 通量。

所以，​HLL 并不保证在每一种状态下都比 Rusanov 更精细​；采用对称速度估计时，两者可以完全相同。

3.3.5 Roe：分解不同的波，再分别处理传播

Roe 的思路是：在界面附近构造一个具有特定守恒性质的局部线性系统，用它近似原来的非线性问题。

对一维 SWE，常用 Roe 平均为：

$$\begin{gathered} u=\frac{u_L\sqrt{h_L}+u_R\sqrt{h_R}}{\sqrt{h_L}+\sqrt{h_R}} \\ h=\frac{h_L+h_R}{2}，c=\sqrt{gh} \end{gathered}$$

两个特征速度为：

$${\lambda }_1=u-c，{\lambda }_2=u+c$$

把状态差分解成两个波：

$$U_R-U_L={\alpha }_1{\mathbf{r}}_1+{\alpha }_2{\mathbf{r}}_2$$

再计算：

$$F_{Roe}=\frac{F_L+F_R}{2}-\frac{1}{2}\sum _{p=1}^2|{\lambda }_p|{\alpha }_p{\mathbf{r}}_p$$

可以把它理解为：每一种波按照自己的传播速度提供相应的数值耗散。

Roe 能较细致地处理波结构，但标准 Roe 方法还存在需要处理的问题，例如跨声速稀疏波的熵修正，以及某些状态下出现非物理负水深。

3.3.6 三种通量应怎样比较

​一个单元界面的平均水面是平的，也不意味着动量通量为零​。静水中仍有压力通量；静止状态依靠两侧压力作用相互抵消来保持。

3.4 时间推进

3.4.1 空间离散后，得到的是变化率

完成空间离散和通量计算后，可以写出：

$$\frac{d{\overline{\mathbf{U}}}_i}{dt}={\mathcal{L}}_i(\overline{\mathbf{U}})$$

其中：

$${\mathcal{L}}_i(\overline{\mathbf{U}})=-\frac{{\mathbf{F}}_{i+\frac{1}{2}}-{\mathbf{F}}_{i-\frac{1}{2}}}{\Delta x}+{\overline{\mathbf{S}}}_i$$

\\mathcal L 是离散空间算子，它负责根据当前状态计算：

• 各界面的流动；
• 单元内水深的变化率；
• 单元内动量的变化率；
• 必要的源项作用。

此时，整个问题成为一组关于时间的常微分方程。

时间推进方法决定：如何利用这些变化率，把当前状态推进到下一时刻。

这种先离散空间、再求解时间变化的方法通常称为线方法。

3.4.2 显式 Euler：使用当前变化率前进一步

显示 Euler 为：

$$U^{(}n+1)=U^n+\Delta t\mathcal{L}(U^n)$$

含义是：假设在这一个短时间段内，状态按照当前变化率变化。它只计算一次空间算子，形式最简单，时间精度为一阶。

使用 3.3.3 中的例子，设两个单元宽度均为：

$$\Delta x=1m$$

它们之间的体积通量约为：

$$F_h=2.215m^2/s$$

若外侧边界没有水流通过，取：

$$\Delta t=0.05s$$

则左单元水深减少：

$$\Delta h_L=-\frac{\Delta t}{\Delta x}{F}_h\approx -0.1107m$$

右单元增加相同水深：

$$\Delta h_R\approx +0.1107m$$

于是：

$$\begin{gathered} h_L^{n+1}\approx 1.8893 \\ {h}_R^{n+1}\approx 1.1107 \end{gathered}$$

这里展示了质量分量；实际求解器还必须同步更新动量分量。

4.2.3 Runge-Kutta：在一个时间步内多次评估变化率

Runge–Kutta 是一类方法的统称。它通过构造中间状态，多次评估变化率，再组合得到下一时刻的结果。

常见的二阶 SSP Runge–Kutta 为：

$$\begin{gathered} U^{(1)}=U^n+\Delta t\mathcal{L}(U^n) \\ {U}^{(n+1)}=\frac{1}{2}{U}^n+\frac{1}{2}\left[U^{(1)}+\Delta t\mathcal{L}(U^{(1)})\right] \end{gathered}$$

第一阶段得到预测状态；第二阶段重新计算预测状态下的变化率，再进行组合。第二阶段必须重新计算通量。

4.2.4 SSPRK3 与“强稳定保持”

三阶 SSP Runge-Kutta 的常见形式为：

$$\begin{gathered} U^{(1)}=U^n+\Delta t\mathcal{L}(U^n) \\ {U}^{(2)}=\frac{3}{4}{U}^n+\frac{1}{4}\left[U^{(1)}+\Delta t\mathcal{L}(U^{(1)})\right] \\ {U}^{(n+1)}=\frac{1}{3}{U}^n+\frac{2}{3}\left[U^{(2)}+\Delta t\mathcal{L}(U^{(2)})\right] \end{gathered}$$

它需要三次空间算子评估。

SSP 的含义是：如果空间离散与前向 Euler 在一定步长限制下具有某种合适的稳定性质，SSP 时间方法可以在相应限制下保留该性质。它不能自动修复本来就不稳定或会产生负水深的空间离散。

3.2.5 三种离散方式的区别

时间阶数更高，并不自动意味着更适合你的实现。还需要考虑稳定性、空间误差、显存访问和总计算成本。

3.2.6 时间精度与空间精度的耦合

假设空间离散只有一阶，而时间推进使用四阶 RK，整体误差仍可能主要来自空间离散：

$$误差\sim C_x\Delta x+C_t\Delta t^4$$

这是平滑解条件下的示意关系。

因此：

• 一阶空间方法加高阶时间方法，不能自动恢复被空间耗散抹掉的波纹；
• 高阶空间方法加过于粗糙的时间推进，也可能损失精度；
• 存在间断时，局部精度还会受到限制器等因素影响。

一个常见的升级方向是：先建立一阶有限体积与显式 Euler 基线，再增加受限线性重构和 SSPRK2。

3.2.7 时间步长的稳定性要求

对常见一维显式格式，CFL 条件具有如下形式：

$$\Delta t\le C_{CFL}\frac{\Delta x}{ma{x}_i(|u_i|+\sqrt{g{h}_i})}$$

这里的 Ccfl 上限依赖具体格式。对于显式 Euler 时间推进，一维情况下 ​Cmax=1；但在多维或高阶格式下，通常取 0.5∼0.9 ​以确保安全。

二维非分裂格式常见的保守步长形式为：

$$\Delta t\le \frac{C_{CFL}}{\frac{a_x}{\Delta x}+\frac{a_z}{\Delta z}}$$

其中 ax, az 是两个方向的最大传播速度估计。

CFL 约束将前面各层联系起来：物理模型决定传播速度，空间网格决定信息需要跨越的距离，时间方法影响允许的稳定步长。显式 SWE 求解器通常据此选择时间步。

对于二维规则网格，如果其他条件不变，将两个方向的网格间距都减半：

• 单元总数约增加到 4 倍；
• 显式时间步通常也需要减半；
• 模拟同样物理时长的基础计算量可能接近 8 倍。