SPH粒子流体模拟：从物理原理到溃坝仿真Python实战

简介：SPH流体仿真项目资源包，面向流体力学模拟、计算机图形学与物理仿真方向的开发者。项目基于光滑粒子动力学（SPH）方法，用C++与OpenGL实现流体运动过程，通过大量粒子求解动力学方程并实时渲染，可帮助学习者理解粒子法流体模拟的完整流程。压缩包共75个文件，包含25个头文件、13个C++源文件，另有若干CUDA源文件与头文件用于GPU加速计算，以及项目工程文件、静态库和动态链接库，能直接支撑编译与调试，整体大小仅5.64MB。目前已有1057人学习下载，适合具备一定C++和图形学基础的读者进行源码研读与实践。资源中提供了完整的工程配置与可执行程序，便于快速运行观察效果，也可在此基础上修改粒子数量、边界条件或外力场，进一步掌握SPH核函数估算、邻居搜索、压强求解等关键环节。

1. 粒子代替网格的流体模拟：SPH 能算什么、不能算什么

第一次接触 sph流体仿真，大多数人都是冲着“水花四溅、飞沫乱飞”来的；真正跑起来之后才会意识到，这个算法的价值不在画面，而在它把流体方程换成了一种完全不同的离散方式：没有网格、没有拓扑重建，每个采样点自带质量和速度，流体的大变形、飞溅和破碎都变成粒子位置的重新排列。SPH 的基本思路是用一组“光滑核”把邻居粒子的物理量加权平均，把连续场方程变成一组可迭代的求和。它特别适合自由表面剧烈变化的问题，比如溃坝、泼溅、液滴撞击液面；缺点也写在明面上——粒子数上去后算力陡增，参数没调好会出现粒子炸飞或扎堆。

这套方法适合三类人：需要自己掌握流体内核的仿真工程师、做物理动画的技术美术、以及研究粒子方法的算法人员。如果你只是想要一个“一键出效果”的流体工具，SPH 不是首选；如果你要的是把物理过程完整握在手里，它值得投入。下文按我自己的落地路径讲：先立住物理，再给最小可运行实现，然后谈参数和排坑，最后说怎么验证结果。

2. 从 Navier-Stokes 到粒子求和：SPH 的物理内核与核函数

2.1 拉格朗日视角：粒子真正承载的是速度场而不是“水分子”

先弄清一个常被误解的点：SPH 粒子不是水分子，它是携带物质属性的体积单元。每个粒子在初始化时继承了质量、位置、速度、密度，后续所有物理量都跟着粒子走。这和网格法的欧拉视角正好相反——欧拉法盯着固定的空间点，流体流过时该点的速度、压力随时间变化；SPH 盯着采样点本身，点移动到哪里，属性就带到哪里。

这个视角带来的第一个好处是质量守恒变得极其自然。每个粒子质量在初始化后保持不变，总质量就是所有粒子质量之和，不存在网格法中“单元间通量计算误差”导致的质量流失。第二个好处是自由表面不需要额外追踪：自由表面就是“有粒子和没粒子”的边界，溃坝里的波浪翻卷、液滴脱离水面，本质上都是粒子集合的重新分布，没有网格畸变问题。

代价在于：粒子的邻居是动态变化的，任何物理量的计算都要靠自己周围的邻居来插值。这就引出了 SPH 的核心工具——核函数。你可以把核函数理解为一种“带半径的加权平均模板”：距离越近贡献越大，超过一定半径直接为零。所有连续方程到这里都被改写成“对邻居求和”的形式，这正是 SPH 名字里“光滑粒子”的由来。

2.2 核函数与光滑长度 h：把连续积分写成粒子求和

SPH 的数学起点是一个恒等式：任意连续场 A(r) 可以用核函数 W 表示成积分，再用粒子离散化成求和：

A(r) ≈ Σ_j (m_j / ρ_j) A_j W(|r - r_j|, h)

其中 W 就是核函数，h 是光滑长度，决定 W 的衰减半径。工程中最常用的是 Wendland C2 核，因为它比传统三次样条更不容易产生粒子成团现象，是近年多数 SPH 实现的首选。它的二维形式如下：

W(q) = α (1 - q/2)^4 (2q + 1), 0 ≤ q < 2 q = r / h α = 7 / (4π h²) # 2D 归一化系数

核函数必须满足三个性质：归一化，权重在支持域内积分等于 1；对称性，W 只取决于距离；紧支撑，距离超过 2h 权重直接为零。紧支撑决定了近邻搜索的半径是 2h，也就是每个粒子只需要跟半径 2h 内的粒子通信，这直接决定了性能模型。

光滑长度 h 是一个全局敏感参数。h 太小，支持域内邻居粒子太少，权重分布抖动大，密度场和压力场都会出现高频噪声；h 太大，所有权重都趋向均匀，流体细节被抹平，看起来像一团“缓慢膨胀的果冻”。常见做法是取初始粒子间距 dp 的 1.2 到 2.0 倍，我一般从 1.5 起步，具体影响在第 4 章展开。

2.3 密度、压力与粘性：最简 SPH 需要哪三项

一套能稳定跑起来的最小 SPH，至少包含三项计算。第一项是密度，直接对邻居求和：

ρ_i = Σ_j m_j W(|r_i - r_j|, h)

注意求和要包含粒子自身，因为 r=0 时 W 取最大值，这是密度公式里不可省略的自贡献项。第二项是压力，弱可压缩 SPH 使用状态方程而不是求解压力泊松方程：

p_i = (c0² ρ0 / γ) [ (ρ_i / ρ0)^γ - 1 ]

这是 Tait 方程，γ 通常取 7，ρ0 是参考密度，c0 是人工声速，量级取预估最大流速的 10 倍。c0 越大，流体越“硬”，可压缩性越小，但时间步长会被迫缩小，计算成本上升。第三项是动量方程，也就是加速度怎么算：

dv_i/dt = -Σ_j m_j (p_i/ρ_i² + p_j/ρ_j²) ∇W_ij + Σ_j m_j Π_ij ∇W_ij + g

第一项是压力梯度力，压力高的粒子把压力低的推开，这是流体流动和波浪产生的主要驱动；第二项是人工粘性项，抑制粒子之间的高频相对运动，防止数值震荡；g 是重力。人工粘性 Π_ij 只在粒子相互靠近时激活，作用类似物理粘性，但它属于数值手段，调整幅度要克制，调大了会让流动显得“稠”。

我不建议一上来就上投影法。虽然不可压缩 SPH 在精度上更优，但每一步要解压力方程组，实现复杂度上升一个量级。先用弱可压缩形式跑通物理流程，确认参数和边界处理没问题，再考虑换更高级的压力求解器，工程节奏更合理。

3. 最小可运行 SPH：用 Python 写一个 2D 溃坝模拟

选择溃坝作为首个场景，原因很直接：初始条件简单，一个大水柱放在矩形区域一侧，撤去支撑后水流自由坍塌，物理现象丰富，又是公开文献里的标准基准测试。下面这套基于 NumPy 的实现是我的教学惯用版本，性能不是目标，但结构完整，能直接跑出肉眼可见的水流效果。

3.1 参数与粒子初始化：先把“水柱子”摆进计算域

先用一组常量把全局参数固定下来，再写初始化函数。粒子参数集中定义，后续调参会很顺手：

import numpy as np # ---------- 模型参数 ---------- DP = 0.02 # 初始粒子间距 0.02 m H = 1.5 * DP # 光滑长度，h/dp=1.5 是常用起点 RHO0 = 1000.0 # 参考密度 kg/m^3 C0 = 20.0 # 人工声速，按最大流速 2m/s 的 10 倍取 GAMMA = 7.0 # 水的 Tait 状态方程指数 ALPHA = 0.1 # 人工粘性系数 G = np.array([0.0, -9.8]) # 重力，2D 场景只有 y 分量 DT = 2e-4 # 时间步长，先保守取小 BOUNDARY_MARGIN = 3 * DP # 边界斥力作用距离 K_BOUNDARY = 500.0 # 边界斥力强度 def init_dambreak(water_w, water_h, domain_w, domain_h): # 在计算域左侧生成一块矩形水体 xs = np.arange(0.5 * DP, water_w, DP) ys = np.arange(0.5 * DP, water_h, DP) gx, gy = np.meshgrid(xs, ys) pos = np.stack([gx.ravel(), gy.ravel()], axis=1).astype(float) # 加微小随机抖动，打破规则排列 pos += np.random.uniform(-0.05 * DP, 0.05 * DP, pos.shape) vel = np.zeros_like(pos) mass = np.full(len(pos), RHO0 * DP * DP) return pos, vel, mass

参数里值得重点说明的是人工声速 C0。很多人误以为它应该取真实水中的声速 1500 m/s，那样时间步长会被压到不可用的程度。SPH 中 C0 是数值参数，取预估最大流速的 10 倍就能把密度变化控制在 1% 以内，同时保证时间步长在合理范围。溃坝实验的最大流速常在两三米每秒，取 20 是个稳妥起点。

质量赋值 m = ρ0 · dp² 对应的是 2D 场景：每个粒子代表一个面积元，而不是体积元。如果你改成 3D，质量要换成 ρ0 · dp³。最后一行的 np.full 确保所有粒子质量一致，这是 SPH 里“每个粒子质量相等”的常见简化——非均匀质量的 SPH 也能做，但没必要在起步阶段给自己加复杂度。

3.2 核函数与近邻搜索：用 KDTree 拿到邻居点对

接下来定义核函数和梯度。我用 Wendland C2，它在粒子分布不均匀时的表现优于三次样条：

def wendland(q): # q = r / h，返回核函数数值 alpha = 7.0 / (4.0 * np.pi * H * H) out = np.zeros_like(q) mask = q < 2.0 qm = q[mask] out[mask] = alpha * (1.0 - qm / 2.0) ** 4 * (2.0 * qm + 1.0) return out def grad_wendland(dx, q): # dx 为从粒子 i 指向粒子 j 的位移向量，返回梯度向量 alpha = 7.0 / (4.0 * np.pi * H * H) g = np.zeros_like(dx) mask = q < 2.0 qm = q[mask, None] r = np.linalg.norm(dx[mask], axis=1, keepdims=True) + 1e-12 # dW/dr 的解析式 dq = -5.0 * alpha * qm * (1.0 - qm / 2.0) ** 3 / H g[mask] = dq * dx[mask] / r return g

梯度公式里的符号要格外小心。dx 从 i 指向 j，dW/dr 为负值（距离越远权重越小），两者相乘后梯度方向实际指向 j。这个方向约定直接关系到压力梯度的正负，后面主循环里所有力的方向都以它为准。

近邻搜索用 SciPy 的 cKDTree 实现。它一次性返回所有距离小于 2h 的点对，省去自己写空间哈希的麻烦：

from scipy.spatial import cKDTree def build_neighbors(pos): tree = cKDTree(pos) # output_type='array' 直接返回 Nx2 的点对数组 return tree.query_pairs(2.0 * H, output_type='array')

query_pairs 返回的点对不含“自己和自己的配对”，这是符合预期的，因为密度求和里粒子的自贡献需要单独补上。用 KDTree 的教学版本在几千到几万粒子时够用；到了十万级，我会换成均匀网格配合空间哈希，原因在第 5 章的近邻搜索瓶颈里展开。

3.3 密度与加速度计算：主循环的每次迭代

密度求和要特别注意自贡献项。每个粒子对自身密度的贡献是 m_i · W(0)，W(0) 在二维 Wendland 下等于 7/(4πh²)，必须手动加上，否则密度整体偏低，压力也随之失真：

def compute_density(pos, mass, pairs): n = len(pos) rho = np.full(n, 0.0) dx = pos[pairs[:, 0]] - pos[pairs[:, 1]] r = np.linalg.norm(dx, axis=1) q = r / H w = wendland(q) # 对称累加：每对邻居的贡献同时加到双方 np.add.at(rho, pairs[:, 0], mass[pairs[:, 1]] * w) np.add.at(rho, pairs[:, 1], mass[pairs[:, 0]] * w) # 自贡献项，最容易漏 rho += wendland(np.zeros(1))[0] * mass return rho

密度求出后，压力由状态方程直接计算，然后组装加速度。压力梯度、人工粘性和重力依次加上：

def compute_acceleration(pos, vel, mass, pairs, domain_w, domain_h): n = len(pos) rho = compute_density(pos, mass, pairs) p = C0 * C0 * RHO0 / GAMMA * ((rho / RHO0) ** GAMMA - 1.0) acc = np.tile(G, (n, 1)) # 重力项 dx = pos[pairs[:, 0]] - pos[pairs[:, 1]] r = np.linalg.norm(dx, axis=1, keepdims=True) q = r / H grad = grad_wendland(dx, q) # 梯度指向 j # 压力梯度力 coef = mass[pairs[:, 1]] * ( p[pairs[:, 0]] / rho[pairs[:, 0]] ** 2 + p[pairs[:, 1]] / rho[pairs[:, 1]] ** 2 ) np.add.at(acc, pairs[:, 0], -coef[:, None] * grad) np.add.at(acc, pairs[:, 1], +coef[:, None] * grad) # 人工粘性，只在粒子逼近时激活 dv = vel[pairs[:, 0]] - vel[pairs[:, 1]] prod = (dv * dx[:, 0]).sum(axis=1) if False else (dv * dx).sum(axis=1) r2 = r.squeeze() ** 2 + 0.01 * H * H mu = H * prod / r2 rho_bar = 0.5 * (rho[pairs[:, 0]] + rho[pairs[:, 1]]) Pi = np.where(prod < 0.0, -ALPHA * C0 * mu / rho_bar, 0.0) np.add.at(acc, pairs[:, 0], mass[pairs[:, 1]] * Pi[:, None] * grad) np.add.at(acc, pairs[:, 1], -mass[pairs[:, 1]] * Pi[:, None] * grad) acc += compute_boundary_force(pos, domain_w, domain_h) return acc

人工粘性里的 prod 是 (v_i − v_j)·(r_i − r_j)，当它小于 0 表示粒子正在靠近，此时叠加一个与逼近速度成正比的阻力。mu 表达式里 0.01h² 是为了避免分母为零的数值保护，也防止距离极近时粘性力发散。ALPHA 取 0.1 是经验值，调大可以让模拟更稳，但水花会显得偏黏。

3.4 时间积分与边界斥力：稳住的最后两环

时间积分采用 kick-drift-kick 的 leapfrog 格式，一次力计算可以推进半步，稳定性和二阶精度都够用：

def run(pos, vel, mass, steps, domain_w, domain_h, save_every=50): pairs = build_neighbors(pos) acc = compute_acceleration(pos, vel, mass, pairs, domain_w, domain_h) for step in range(steps): vel += 0.5 * DT * acc pos += DT * vel pairs = build_neighbors(pos) acc = compute_acceleration(pos, vel, mass, pairs, domain_w, domain_h) vel += 0.5 * DT * acc if step % save_every == 0: np.save(f"pos_{step:06d}.npy", pos)

边界处理是新手最容易翻车的位置。纯斥力函数虽然简单，但只适合在低压场景用；冲击壁面时压力可能把粒子从斥力区间里硬挤出去。我常用配合式：内壁外侧放两层固定边界粒子，让它们参与密度和压力计算，流体粒子靠近时自然感受到压力升高而被推开，而不是依赖解析斥力：

def compute_boundary_force(pos, domain_w, domain_h): f = np.zeros_like(pos) margin = BOUNDARY_MARGIN # 左墙：粒子离左侧内壁越近，向右的推力越强 dist = pos[:, 0] - 0.5 * DP m = dist < margin f[m, 0] += K_BOUNDARY * (1.0 - dist[m] / margin) # 右墙 dist = (domain_w - 0.5 * DP) - pos[:, 0] m = dist < margin f[m, 0] -= K_BOUNDARY * (1.0 - dist[m] / margin) # 底墙与顶墙同理，沿 y 方向 dist = pos[:, 1] - 0.5 * DP m = dist < margin f[m, 1] += K_BOUNDARY * (1.0 - dist[m] / margin) dist = (domain_h - 0.5 * DP) - pos[:, 1] m = dist < margin f[m, 1] -= K_BOUNDARY * (1.0 - dist[m] / margin) return f

把上面两个函数拼进 run，一个最小的二维溃坝就能跑起来。初期你会看到水柱先微微松弛，然后顶端开始翻卷，底部水流加速冲向右侧壁面——这就是 SPH 对强自由表面流动的直接响应。若出现粒子乱飞，先不要动核函数，按第 5 章的排查顺序走一遍。

4. 参数怎么定：光滑长度、粒子间距、声速与时间步长的匹配

SPH 调参没有“一套通吃”的黄金组合，但参数之间有强耦合关系。我把最影响结果的四个参数单独拎出来，讲清楚它们各自的职责，以及调一个参数连带影响了什么。

4.1 h/dp 比值：从 1.2 到 2.0，细节和稳定性的交换

光滑长度 h 与粒子间距 dp 的比值是所有参数里最核心的一个，它决定每个粒子支持域内的平均邻居数量。二维里邻居数约等于 (2h/dp)² 再乘一个密度因子，h/dp=1.5 时每个粒子大约有几十个有效邻居，这个数量足以让密度和压力场平滑。

h/dp 取 1.2 时，支持域内粒子偏少，压力场容易在局部出现空洞或尖峰；好处是细节更锐利，液滴脱离、薄水舌这类小尺度结构保留得更好。取 2.0 时数值更稳，粒子成团倾向减弱，但薄片结构会被过度平滑，看起来偏“黏”。我的习惯是先取 1.5 跑基准，再把 1.2 和 2.0 各跑一遍对比液面前沿位置曲线，视差值决定往哪边靠。这个对比实验成本很低，但能避免拍脑袋定 h。

还有一个容易忽略的点：h 是全局固定值，还是可以随密度自适应？自适应 h 能改善粒子稀疏区域的表现，但会让核函数失去对称性，动量守恒变差。除非粒子分布跨度极大，否则我不会在第一个项目里用自适应 h。

4.2 粒子间距与规模感：从 1 万到 100 万粒子怎么换算

粒子间距 dp 决定分辨率，也直接决定总粒子数。二维中粒子数 N 与面积成正比，三维中与体积成正比，dp 减半时三维粒子数变成 8 倍。拿一个 1m×0.6m 的二维水槽估算：dp=0.02m 时大约 1500 个粒子，演示没问题；dp=0.005m 时直接跳到 24000 个，计算量和内存同步上涨。

我的经验是先把规模分成三个档位来规划。教学和原理验证在几万粒子内，用 CPU 单线程配 NumPy 足够；十到几十万粒子的工程测试，需要上 C++/Rust 配合均匀网格近邻搜索；百万级场景基本绕不开 GPU 或大规模并行。下表给出我自己常用的规模参照：

这里的“内存压力”不是指粒子数组本身——一个二维粒子含位置、速度、密度、压力，约 5 个 float64 数组，每个粒子 40 字节（三维会更多）；压力主要来自近邻点对数组，它的长度能到粒子数的几十倍，才是内存消耗的大头。

4.3 人工声速与 CFL 条件：压缩性和时间步长的拉扯

弱可压缩 SPH 里，人工声速 c0 决定了流的“軟硬程度”。c0 取 10 倍最大流速，对应马赫数 0.1，密度振荡控制在 1% 量级。提高 c0 可以进一步逼近不可压缩，但声速越高，信息传播越快，显式时间积分被迫用更小的 dt，计算量随 c0 近似线性上涨，这笔账要把模拟时长一起算进去。

时间步长受 CFL 条件约束，常用判据是：

dt ≤ 0.25 * h / (c0 + v_max)

v_max 是当前最大流速。这个式子能解释大部分“模拟爆炸”事故：水柱刚倒下时流速快速增大，如果 dt 按初始流速算好就不再调整，一旦瞬时流速超过设计值，CFL 条件被破坏，局部分离会在一两个时间步内指数放大。

我把常用的起步参数整理成一张速查表，方便照抄：

参数之间存在联动：比如 c0 调大一倍，CFL 条件里的 vmax 贡献相对变小，但绝对限制更严，DT 往往要同步减半；ALPHA 调大后，高频噪声被压住，看似稳定了，但真实的细小涡旋也被抹掉。因此调参时要有“一次只动一个”的纪律，并且永远用同一个基准场景对比。

5. SPH 模拟的 4 个常见坑：从粒子爆炸到边界穿透

SPH 的坑很有辨识度：物理引擎报错会告诉你“穿模了”“飞出去了”，SPH 则是直接给出一堆 NaN 或者一帧比一帧离谱的粒子云。下面四条是我自己踩过也帮别人排查过的高频问题，按现象、原因、解决三步写。

5.1 粒子炸飞：密度尖峰后的压力雪崩

现象：模拟跑到几十步，某几个粒子速度突然变成几百甚至几千，随后整个粒子场在几步内全面发散，数值变成 NaN。断点调试时能看到爆炸前一个时间步的密度数组里出现异常大值。

原因：压力状态方程里 (ρ/ρ0)^γ 的 γ=7，密度一旦超过参考密度，压力以七次方放大。密度尖峰最常见的来源是粒子初始抖动过大导致支持域内粒子重叠，或者是 dt 太大导致相邻粒子在一步内穿过彼此。自由表面处粒子密度偏低产生负压，又把周围粒子吸附过来，形成局部坍缩。

解决：先把 dt 减半，看是否还炸；如果还炸，检查初始化抖动幅度，把 0.05dp 改成 0.01dp；再不行就检查核函数支持域内是否有粒子缺失，比如水柱角点处粒子天然稀疏，此时需要稍微加大 h。最后的手段是给压力做一个上下限截断，但那是治标，只能让模拟“不炸”，不能让你理解物理，不推荐新手优先用。

5.2 粒子扎堆：拉伸不稳定性在捣乱

现象：模拟中后期，粒子在流体内部或自由表面附近聚成小团，形成规则晶格状的分布，密度场出现与初始网格结构一致的花纹。这不是真实物理，是 SPH 的经典数值病——拉伸不稳定性。

原因：在受拉伸区域（负压环境），SPH 的力对粒子位置扰动呈正反馈，规则排列被不动点吸引，粒子趋于成团。早期三次样条核函数更容易触发这个现象，这也是我改用 Wendland 核的原因之一。粒子间距相对 h 过密时也更明显。

解决：优先把 h/dp 从 1.5 提到 1.8 或 2.0，让权重更平滑；再把人工粘性略微提高，抑制局部相对速度。如果还不行，检查 c0 是否过高——c0 越大，压力对密度误差越敏感，负压抖动越剧烈。换成更高阶的核函数或加密度修正项也能缓解，但代码复杂度会上去，建议作为第二步。

5.3 边界穿透：斥力挡不住高速冲击

现象：水柱坍塌后高速冲击底部或侧壁，一部分流体粒子直接穿出计算域，卡在边界外。用固定斥力函数的模拟尤其容易出现，因为斥力只在粒子进入边界的某个区间内生效，速度足够大时一个时间步内就能穿过整个区间。

原因：纯斥力的本质是“事后补救”，它不参与流体的密度和压力计算，相当于一道没有物理支撑的墙。粒子速度接近或超过 margin/dt 时，斥力还没来得及起作用，粒子已经越过边界。

解决：最可靠的做法是改用固定边界粒子。在域边界外侧布置二到三层不参与积分的虚拟粒子，让它们正常参与流体粒子的密度和压力求和。流体粒子接近边界时，密度真实升高，压力真实增大，自然被推回——这比任何解析排斥力都符合 SPH 的物理逻辑。边界粒子层的厚度至少要覆盖核函数支持域，三层是稳妥起步。

5.4 近邻搜索拖垮性能：O(N²) 的死胡同

现象：粒子数从 1 万翻到 4 万，运行时间不是涨 4 倍，而是涨将近 10 倍。用 KDTree 的教学代码在粒子数上升后每步重建树的开销越来越大，如果还为每个粒子单独查球，就更慢。

原因：近邻查询的复杂度与支持域半径、粒子密度相关，但实现不当很容易退化成全局搜索。每步重建 KDTree 本身也有常数开销，在十万粒子量级会明显拖后腿。

解决：换均匀网格空间哈希，网格尺寸直接取 h，每个粒子只需要查询自己所在网格及周围 3×3（二维）或 3×3×3（三维）的网格。网格大小和近邻搜索只遍历邻近网格，复杂度降到 O(N)，这是 SPH 工程实现的标配。静态边界粒子可以预先建好网格，每步只重建流体粒子部分。

6. 往工程走：GPU 加速、验证实验与参数标定的习惯

6.1 三个验证实验：静水压、溃坝、能量曲线

跑通模拟只是第一步，真正判断“做对了”要依赖验证手段。我习惯先做三个实验：静水压测试、溃坝对照、能量曲线检查。

静水压测试很简单，把水体放在容器里静置几千步，检查底部压力是否符合 ρgh 的理论值，且粒子速度长期保持接近零。这个测试能暴露核函数、状态方程和边界处理的系统性问题，比一上来就跑动态场景省心得多。溃坝对照是用公开实验数据里水柱前沿位置随时间的变化，跟自己的模拟结果画在同一张图上，看趋势是否一致，这是判断整体参数是否合理的金标准。能量曲线则是记录系统动能加势能：

E_kin = 0.5 * (mass * np.sum(vel ** 2, axis=1)).sum() E_pot = (mass * (-G[1]) * pos[:, 1]).sum() E_total = E_kin + E_pot

SPH 的总能量不严格守恒，因为人工粘性会耗散能量，但曲线应该平滑下降而不是上下剧跳。如果总能量出现明显回升，多半是时间步长或边界处理有问题，值得回头排查。

6.2 从 CPU 到 GPU：性能上台阶的正确姿势

粒子规模到了十万级，CPU 版本的速度会让人失去耐心。常见的加速路径是把近邻搜索和密度、压力、力的累加写成 GPU kernel，同一份数据在 GPU 上的计算速度通常比 CPU 高一个数量级。现有开源 SPH 框架大多提供 CPU 和 GPU 双后端，但我不建议直接跳到麻烦的 framework——先用第 3 章的最小实现理解每个阶段的数据流动，再迁移到并行版本，遇到瓶颈时你能清楚知道是近邻写入冲突还是内存带宽不足。

6.3 参数标定的第一个习惯

我把压箱底的习惯写在这里：每一个需要调试的项目，永远先建立一个“最小复现场景”，它要足够小到几十秒内能跑完，又足够典型到能复现你关心的物理现象。参数全部放进配置字典，每次只改一个，记录下它前后的能量曲线和液面前沿，别靠肉眼判断“看起来差不多”。我吃过不少亏，最典型的一次是同时调了 h 和 ALPHA，结果模拟从炸飞到死水都找不到原因，最后只能重跑分项排查。血泪经验换来的教训就是：参数实验是科学实验，要有对照和记录。希望帮到你。

本文还有配套的精品资源，点击获取