1. 项目概述：从“水往低处流”到数学方程
我们常说“水往低处流”，这背后其实是一整套复杂的物理规律在起作用。当水流过一片宽阔的浅滩，或者洪水漫过平原时，描述其运动的核心工具，就是“浅水方程”。这个项目，我们要做的，就是亲手把“水往低处流”这个朴素的直觉，一步步推导成严谨的二维浅水方程数学模型。这不仅仅是理论推导，更是理解洪水演进、河口潮汐、甚至大气环流等众多自然和工程现象的基础。如果你对流体力学感兴趣，或者从事水利、环境、气象相关的工作，掌握这个方程的来龙去脉，就等于握住了打开一扇大门的钥匙。它看起来复杂，但拆解开来，每一步都有清晰的物理图景支撑。接下来，我会带你从最基本的物理定律出发，看看我们如何用数学语言，精确地描述一片浅水区域的流动。
2. 核心思路拆解：守恒律是基石
推导任何流体力学方程，核心思路都离不开“守恒律”。对于一片水，什么东西是守恒的？质量不会凭空产生或消失，动量（可以粗略理解为“运动的量”）也不会。所以，我们的推导将紧紧围绕“质量守恒”和“动量守恒”这两大物理学基石展开。
2.1 物理图景与核心假设
在动手写方程之前，我们必须先明确我们研究的对象是什么样子，以及做了哪些简化，这直接决定了方程的适用范围。
首先，什么是“浅水”？这里的“浅”不是绝对深度，而是相对于水平尺度而言。想象一个巨大的湖泊或者一片被洪水淹没的广阔平原，它的水平方向可能有几公里甚至上百公里，而水深可能只有几米。这种水平尺度远大于垂直尺度的流动，就是“浅水”流动的典型场景。基于此，我们引入浅水理论最核心的假设：静水压力近似。这意味着，在垂直方向上，我们认为水体的压力分布和静止时一样，只由该点上方水柱的重量决定。这个假设极大地简化了问题，因为它允许我们将三维的流动，在垂直方向上进行“积分”或“平均”，最终降维成二维（仅水平方向）的方程。另一个关键假设是，我们考虑的水体是不可压缩的，密度为常数。对于大多数地表水问题，这个假设是合理的。
所以，我们的物理图景是：一片广阔且相对较浅的水域，水深随着位置和时间变化，我们关注的是水面高度和水平方向流速如何演变。
2.2 控制体方法与积分形式
推导守恒方程，工程师和物理学家最爱用的工具就是“控制体”。我们不是在跟踪某一滴水的命运（拉格朗日视角），而是固定空间中的一个小区域（欧拉视角），看看有多少质量、动量流进流出。我们选取一个底面积为ΔxΔy，高度从底部地形b(x,y)到自由水面η(x,y)的柱状水体作为控制体。这里，h(x,y,t) = η(x,y,t) - b(x,y) 就是我们常说的水深。
对于这个控制体，质量守恒可以表述为：控制体内质量的增加率 = 流入控制体的质量流量 - 流出控制体的质量流量。动量守恒则表述为：控制体内动量的增加率 = 流入的动量流量 - 流出的动量流量 + 作用在控制体上的合外力。合外力通常包括压力、重力和底部摩擦力。
我们将首先针对这个固定的控制体，写出质量和动量守恒的积分形式方程。这是最直观、物理意义最清晰的步骤。
3. 质量守恒方程（连续性方程）推导
我们从质量守恒开始，这是相对简单的一步，但至关重要。
3.1 建立控制体与质量计算
考虑一个固定在空间中的微小控制体，其底面在水平面（x-y平面）上的投影是一个小矩形，边长分别为Δx和Δy。控制体的下边界是固定的河床或地形b(x, y)，上边界是自由水面η(x, y, t)，因此控制体的瞬时高度（水深）为 h(x, y, t) = η(x, y, t) - b(x, y)。假设水的密度ρ为常数。
那么，在时刻t，这个控制体内所含的水体总质量M为： M = ρ * 体积 = ρ * [h(x, y, t) * Δx * Δy]
注意，这里我们隐含了另一个浅水近似：水平速度在垂直方向上均匀分布。也就是说，我们用一个深度平均速度u(x,y,t)和v(x,y,t)来代表整个水柱在x和y方向上的运动。这是将三维问题二维化的关键。
3.2 质量通量与净流入量
质量通过控制体的侧面流入流出。我们以x方向为例。在控制体的左侧面（x处），单位时间流入的质量（质量流量）为：流入质量 = ρ * [u(x, y, t) * h(x, y, t) * Δy]。这里，u是x方向的深度平均速度，u*h可以理解为单宽流量（单位宽度上的体积流量）。
同理，在控制体的右侧面（x+Δx处），单位时间流出的质量为：流出质量 = ρ * [u(x+Δx, y, t) * h(x+Δx, y, t) * Δy]。
因此，在x方向上，净流入控制体的质量流量为： 净流入_x = ρ * [u(x,y,t)h(x,y,t) - u(x+Δx,y,t)h(x+Δx,y,t)] * Δy = -ρ * [ ∂(u h)/∂x * Δx ] * Δy （应用了泰勒展开近似）
同理，在y方向上，净流入控制体的质量流量为： 净流入_y = -ρ * [ ∂(v h)/∂y * Δy ] * Δx
注意：这里符号很重要。我们定义流入为正，流出为负。右侧面速度u为正时表示流出，所以表达式是“左侧流入减去右侧流出”，利用泰勒展开后自然得到负的偏导数项。
3.3 积分形式到微分形式
根据质量守恒定律：控制体内质量的增加率 = 从所有表面净流入的质量流量。 即： ∂M/∂t = 净流入_x + 净流入_y
将各项代入： ∂[ρ h Δx Δy]/∂t = -ρ Δy [∂(u h)/∂x] Δx - ρ Δx [∂(v h)/∂y] Δy
两边同时除以控制体的底面积(ρ Δx Δy)（注意ρ为常数可约去），得到： ∂h/∂t + ∂(u h)/∂x + ∂(v h)/∂y = 0
这就是二维浅水方程组的连续性方程（质量守恒方程）。它的物理意义非常直观：某一点水深的局部变化率（∂h/∂t），是由该点两个方向上的单宽流量散度（∂(uh)/∂x + ∂(vh)/∂y）决定的。如果流出的水多于流入的，该点水深就会下降，反之则增加。
4. 动量守恒方程推导
动量守恒的推导比质量守恒复杂，因为它涉及矢量，并且外力项需要仔细处理。我们分别推导x和y方向的动量方程。
4.1 x方向动量变化率与对流项
同样考虑之前的控制体。控制体内x方向的总动量为： P_x = ρ * (体积 * u) = ρ * (h Δx Δy) * u
其随时间的变化率为：∂P_x/∂t = ρ ∂(h u)/∂t * Δx Δy。这里出现了新的因变量h*u，它代表单位面积上的x方向动量。
动量也会随着水流流入流出控制体，这称为动量对流。在左侧面（x处），流入的x方向动量流量为：[ρ * (u h Δy)] * u(x,y,t) = ρ u (u h) Δy。注意，这里是“质量流量乘以速度”。同理，右侧面流出的动量流量为：ρ u (u h) |{x+Δx} Δy。 因此，x方向净流入的对流动量为：-ρ ∂[u (u h)]/∂x * Δx Δy。 y方向（前后面）也会带来x方向的动量对流。前面（y处）流入：ρ u (v h) Δx；后面（y+Δy处）流出：ρ u (v h) |{y+Δy} Δx。净流入为：-ρ ∂[u (v h)]/∂y * Δx Δy。
所以，对流项带来的总净流入动量为：-ρ [ ∂(u u h)/∂x + ∂(u v h)/∂y ] Δx Δy。
4.2 压力项推导（静水压力近似的应用）
这是浅水方程推导中最巧妙也最关键的一步。作用在控制体上的外力在x方向的分量，主要来自压力。控制体受到的x方向压力，是作用在左右两个侧面的压力差。
在静水压力近似下，水中任意一点(x,y,z)的压力p等于该点上方水柱的重量：p(z) = ρ g [η(x,y,t) - z]。其中g是重力加速度，z是该点的垂直坐标。这意味着压力只随深度线性变化，同一铅垂线上，压力分布和静水时一模一样。
现在计算左侧面受到的压力合力。左侧面是一个垂直平面，其上各点的水深不同，压力也不同。我们需要对这个侧面进行积分。在左侧面（x处），对于从底部b到水面η的任意高度z，该微元受到的压强为ρ g (η - z)，微元面积为Δy * dz。这个压力的方向是沿x轴正方向（指向控制体内）。因此，左侧面受到的总压力为： F_left = ∫_{z=b}^{η} [ρ g (η - z)] dz * Δy 计算这个定积分：∫_{b}^{η} (η - z) dz = [ηz - z^2/2]_{b}^{η} = (η^2 - η^2/2) - (ηb - b^2/2) = (1/2)η^2 - ηb + (1/2)b^2 = (1/2)(η - b)^2 = (1/2) h^2。 所以，F_left = (1/2) ρ g h^2 * Δy。
同理，右侧面（x+Δx处）受到的压力，方向是沿x轴负方向（指向控制体外），其大小为：F_right = (1/2) ρ g h^2 |_{x+Δx} * Δy。
因此，作用在控制体上x方向的净压力为： F_pressure_x = F_left - F_right = -Δy * [ (1/2) ρ g h^2 |{x+Δx} - (1/2) ρ g h^2 |{x} ] ≈ -Δy * ∂[ (1/2) ρ g h^2 ]/∂x * Δx = -ρ g Δx Δy * h (∂h/∂x)。
最后一步用到了链式法则：∂(h^2/2)/∂x = h (∂h/∂x)。这个结果非常优美，它表明净压力驱动项与水深h和水面坡度∂h/∂x成正比。水面越高、坡度越陡，产生的驱动力就越大。这正是“水往低处流”的数学体现。
4.3 重力与底床摩擦力
重力是体积力，垂直向下，在水平x方向没有分量。所以重力不直接贡献于x方向的动量方程。
底部摩擦力是另一个重要的外力项。水流在流经河床或地表时，会受到一个与流动方向相反的阻力，即底床剪应力τ_b。通常采用经验公式来描述，最常用的是曼宁公式或谢才公式。例如，采用曼宁公式，x方向的底床摩擦力可以表示为： F_friction_x = - ρ g n^2 u √(u^2+v^2) / h^{1/3} * (Δx Δy) 其中n是曼宁粗糙系数。这个力是作用在整个控制体底面积上的，方向与流速u相反（故为负号）。这是一个经验项，也是方程中非常重要的耗散项，它决定了流动的最终平衡状态。
4.4 整合得到x方向动量方程
现在，我们将所有项代入x方向的动量守恒定律：控制体内动量的增加率 = 净流入的对流动量 + 合外力。 即： ρ ∂(h u)/∂t ΔxΔy = -ρ [ ∂(u u h)/∂x + ∂(u v h)/∂y ] ΔxΔy + (-ρ g h ∂h/∂x ΔxΔy) + (-ρ g n^2 u √(u^2+v^2) / h^{1/3} ΔxΔy)
两边同时除以ρ Δx Δy，整理后得到： ∂(h u)/∂t + ∂(u u h)/∂x + ∂(u v h)/∂y = -g h ∂h/∂x - g n^2 u √(u^2+v^2) / h^{1/3}
利用连续性方程∂h/∂t = -[∂(uh)/∂x + ∂(vh)/∂y]，可以对上述方程左边进行展开和化简，得到一个更常见的形式。具体地： 左边 = h ∂u/∂t + u ∂h/∂t + ∂(u^2 h)/∂x + ∂(u v h)/∂y 将∂h/∂t用连续性方程替换，并展开导数项，经过一系列运算（这是推导中需要耐心完成的代数步骤），可以消去一些项，最终得到： ∂u/∂t + u ∂u/∂x + v ∂u/∂y = -g ∂η/∂x - g n^2 u √(u^2+v^2) / h^{4/3}
这个形式更为简洁，左边是速度u的物质导数（代表跟随一个水粒子其速度的变化率），右边是驱动力（水面坡度）和阻力（底床摩擦）。注意到这里用了∂η/∂x，因为h = η - b，且假设底床b不随时间变化，所以 -g h ∂h/∂x = -g ∂(h^2/2)/∂x，而当底床坡度平缓时，可以近似为 -g ∂η/∂x。这更直观地表明，流动的驱动力来自于自由水面的坡度。
实操心得：在推导动量方程时，最容易出错的地方是对流项的展开和合并。一个实用的技巧是，始终将因变量视为
h和q_x=u*h、q_y=v*h（即单宽流量）。这样，动量方程可以写成关于q_x和q_y的守恒形式，数值计算时更稳定。上面的推导最终化简成的速度形式（物质导数形式）物理意义清晰，但守恒形式更适合用于构建数值格式。
5. 方程组的完整形式与物理意义
将两个方向的动量方程与连续性方程写在一起，就构成了完整的二维浅水方程组。通常写作以下形式：
连续性方程： ∂h/∂t + ∂(q_x)/∂x + ∂(q_y)/∂y = 0 其中，q_x = u h,q_y = v h。
x方向动量方程（守恒形式）： ∂(q_x)/∂t + ∂(u q_x + g h^2/2)/∂x + ∂(v q_x)/∂y = g h S_{0x} - g h S_{fx}y方向动量方程（守恒形式）： ∂(q_y)/∂t + ∂(u q_y)/∂x + ∂(v q_y + g h^2/2)/∂y = g h S_{0y} - g h S_{fy}
这里我们引入了两个源项：
S_0 = (-∂b/∂x, -∂b/∂y)：底床坡度源项。因为 -g h ∂h/∂x = -g h ∂(η-b)/∂x = -g h ∂η/∂x + g h ∂b/∂x。其中-g h ∂η/∂x是水面坡度驱动力，而g h ∂b/∂x就是底床坡度产生的力。在缓坡假设下，有时会将二者合并为-g ∂η/∂x。S_f = (S_{fx}, S_{fy})：摩擦坡度源项，即前面推导的摩擦阻力项，例如曼宁公式：S_{fx} = n^2 u √(u^2+v^2) / h^{4/3}。
这个方程组是一组非线性双曲型偏微分方程。它的物理意义非常明确：
- 连续性方程：保障了水体的“来龙去脉”清晰，质量不灭。
- 动量方程：左边三项合起来代表了动量的局地变化和对流输运；右边第一项
g h S_0是驱动项（重力的分量，由水面和底床坡度产生），是流动的“发动机”；右边第二项-g h S_f是阻力项（底床摩擦），是流动的“刹车”。
它们共同描述了浅水流动中，水深和流速如何在水面坡度（压力梯度）的驱动下，克服底床摩擦，并通过对流过程相互作用、演化的全部动力学。
6. 数值求解的挑战与核心环节
推导出方程只是第一步，要想用它来模拟真实的洪水或潮汐，必须通过数值方法进行求解。这个过程充满了挑战，也是将理论应用于实践的关键。
6.1 方程的双曲性与特征线
浅水方程是双曲型的，这意味着信息以有限的速度（即波速）传播。对于浅水波，这个波速是√(g h)。这个特性导致了两个重要现象：激波（如潮涌、水跃）和稀疏波。在数值上，这要求我们采用能捕捉激波、保持物理量单调性的格式，比如Godunov类型的格式（Roe, HLLC等）。如果使用不适合的格式（如中心差分），在激波附近会产生非物理的数值振荡，导致解完全失真。
6.2 源项的处理：底坡与摩擦
方程右边的源项处理不当，会导致数值解在平衡状态下无法保持，产生虚假的流动。例如，在静止水体（湖面）情况下，水面是水平的（∂η/∂x=0），底床可以有坡度。此时，驱动项-g h ∂η/∂x和底坡源项g h ∂b/∂x应该精确平衡，合外力为零。这称为“静水平衡”。许多数值格式需要特殊处理（如Well-Balanced Scheme）才能保持这种平衡，否则计算中会出现即使水面初始是平的，也会产生虚假流动的荒谬结果。
底床摩擦项S_f是高度非线性的（与速度的平方成正比，与水深的高次方成反比）。在显式时间推进格式中，它会对时间步长施加非常严格的稳定性限制。通常采用隐式或半隐式方法处理摩擦项，以允许使用更大的时间步长。
6.3 干湿边界处理
实际地形中，水域边界是随时间变化的（比如洪水淹没范围扩大或缩小）。计算域内会存在h=0或h极小的“干单元”。在这些区域，方程会出现奇异性（例如摩擦项分母为0），计算会崩溃。因此，必须设计稳健的“干湿处理”算法。常见的策略包括：
- 设定一个极小水深阈值（如10^-6 m），当单元水深低于该阈值时，将其视为“干”，不计算动量方程，或固定其流速为零。
- 在通量计算中考虑干湿边界，确保质量不会从干单元“渗漏”到更干的单元，动量通量也要做相应限制。
- 处理干湿边界时，还要注意保证质量守恒和动量守恒，这是一个非常棘手但必须解决的问题。
7. 常见问题与排查技巧实录
在实际推导和后续的数值实现中，会遇到各种问题。以下是一些典型问题及解决思路。
7.1 推导过程中的符号混乱
问题：在推导压力项或对流项时，正负号容易搞混。排查技巧：
- 牢记物理图景：对于压力项，水面坡度向下游（x正方向）为正时，压力合力应指向x正方向（推动水流）。所以如果∂η/∂x为负（水面沿x升高），力应该是正的。检查你的方程是否满足：
-g ∂η/∂x，当∂η/∂x为负时，该项为正。 - 检查量纲：每一项的量纲必须一致。动量方程左边是[L/T^2]（速度变化率），右边
g ∂η/∂x的量纲是[L/T^2] * [1] = [L/T^2]，摩擦项g n^2 u^2 / h^{4/3}的量纲也是[L/T^2]。通过量纲分析可以快速发现明显的系数错误。 - 验证特例：用简单的特例验证方程。例如，对于一维均匀定常流（∂/∂t=0, ∂/∂x=0, v=0），动量方程应简化为
S_f = S_0，即摩擦坡度等于底床坡度，这正是曼宁公式描述的情形。如果你的方程不能退化到这个特例，那肯定推导有误。
7.2 数值模拟中的不稳定与崩溃
问题：程序运行时，水深或速度出现NaN（非数字）或异常大的值，计算迅速崩溃。排查技巧：
- 检查初始条件：初始水深场
h必须处处大于0。即使是很小的正数（如0.001m）也比0好。初始速度场最好从静止开始。 - 检查时间步长：双曲方程有严格的CFL稳定性条件：
Δt ≤ CFL * Δx / (|u| + √(g h))，其中CFL数通常小于1（如0.5）。确保你的Δt满足所有网格单元中最严格的条件。计算波速√(g h)时，注意h不能为0。 - 输出中间状态：在崩溃前的一个或几个时间步，将关键变量（h, u, v，通量，源项）输出到文件或屏幕。查看是哪个网格单元先出现异常值，以及异常出现前这些变量的状态。这能帮你定位问题是在通量计算、源项计算还是干湿处理环节。
- 摩擦项处理：如果使用了显式处理摩擦项，尝试将其改为半隐式处理。显式摩擦项要求的稳定时间步长可能比对流项要求的CFL条件还要小得多。
7.3 结果不物理：虚假流动与质量不守恒
问题：模拟一个静止的湖，理论上应该没有任何流动，但结果却出现了速度。排查技巧：
- 静水平衡测试：这是检验格式是否“和谐”的试金石。设置一个水平水面（η=常数），但底床有起伏（b变化）。运行模拟，理论上速度应始终保持为零。如果出现了流动，说明你的格式没有很好地平衡压力梯度项和底坡源项。你需要检查离散格式是否满足“C-性质”或采用专门的Well-Balanced格式。
- 检查边界条件：边界条件设置错误是引入虚假流动的常见原因。对于封闭边界（如岸壁），法向速度应为零。确保你的边界条件实现正确，没有在边界处引入非零通量。
- 质量守恒检查：计算整个域内总水量的变化率。理论上，对于封闭系统（无源汇、无开边界），总水量应严格守恒（仅受机器舍入误差影响）。在每一步或每隔若干步，计算
∑(h * Δx * Δy)，看其变化是否在可接受的误差范围内。如果质量不守恒，问题可能出在通量计算或干湿边界处理上，导致有“漏”水。
7.4 干湿边界处的异常
问题：在水陆交界处，出现“水往高处流”或干区被异常淹没/侵蚀。排查技巧：
- 阈值敏感性测试：调整干湿判断的水深阈值（如从1e-6调到1e-5），观察结果是否发生剧烈变化。如果变化很大，说明算法在干湿边界处不够稳健。一个更稳健的方法是采用“薄层水”处理，即使单元被标记为“干”，也保留一个极小的水深值用于计算波速，但将其速度设为零且不参与通量计算。
- 通量限制：在干湿边界相邻的单元之间计算通量时，需要进行限制。例如，当相邻单元一个很湿、一个很干时，不能直接使用基于两个单元状态计算的通量公式（如Roe通量），因为这可能导致质量从干单元“吸”向湿单元。应采用HLL或HLLC等格式，并妥善处理干湿界面处的波速估计。
- 地形数据精度：检查你的底床高程数据
b。如果地形数据分辨率不够或存在异常噪声，在干湿边界附近会产生虚假的陡坡，驱动不真实的流动。对地形数据进行适当的平滑预处理有时是必要的。
推导二维浅水方程是一个将物理直觉、数学工具和工程实践紧密结合的过程。从最基本的守恒定律出发，通过合理的简化假设，最终得到一套可以描述复杂浅水流动的方程，这个过程的每一步都充满了工程思维的魅力。而将其付诸数值计算，更是对理论理解的深度考验。我个人的体会是，亲手推导一遍胜过读十遍现成的公式，因为在推导中遇到的每一个疑问和障碍，都会迫使你去深入思考其物理本质。而在数值实现阶段，从简单的、理想化的算例（如静水平衡、溃坝波）开始，逐步增加复杂性，是调试代码、理解格式特性的不二法门。最后，永远不要完全相信“黑箱”模型的结果，用这些基本的物理原理和排查技巧去审视你的模拟结果，是成为一个合格的流体模拟工程师的关键。