Matlab实现2D SPH流体模拟：从溃坝到粒子法入门

自从把2D SPH（光滑粒子流体动力学）用Matlab完整实现一遍之后，我最大的感受就一句话：流体模拟的入门门槛，真没有你想象中那么高。不需要动辄上千行的C++工程，也不用去啃完整的数值传热学教材，一套几十个函数的Matlab脚本，就能眼睁睁看着一堵水柱在重力作用下垮塌、飞溅、撞击容器壁然后反卷回来。这种直观的“即时满足”感，是我入门流体模拟以来最上头的体验。

这个项目本身是经典的2D SPH模拟器：大约两三千个粒子，每个粒子都携带位置、速度、密度和压力信息，通过核函数互相“感受”周围邻居的存在，从而计算压力梯度、粘性力和重力，再按时间步推进位置。因为是无网格的拉格朗日方法，自由表面、飞溅、破碎这些在网格法里很难追的事情，在这里天然地就出现了——你不需要专门的表面追踪算法，界面本身就是粒子群的几何外轮廓。

这个附带的Matlab代码特别适合三类人：正在学CFD但不想一上来就啃Fortran/C++的数值方法入门者；需要做毕业设计、想快速做出可视化效果的本科生或硕士生；以及做游戏开发、动画特效，想验证一下粒子流体思路的“手感党”。只要你的Matlab版本不是太老（建议R2019之后），随便一台普通笔记本都能跑起来，五千个粒子以内的规模完全不需要服务器。

下面我按自己的理解，把这个2D SPH模拟从原理到代码再到调参的心路历程全部拆开讲一遍。

1. 这套SPH代码到底在干什么

1.1 你看到的画面是什么样的

整套代码跑起来之后，默认是一个矩形水箱，水箱的左侧或者左下角有一个矩形的“水方块”。这堆水方块本质上是由几百上千个粒子按网格规则排列组成的，粒子与粒子之间保持一个均匀间距，比如0.01米。当你按下运行键，重力开始作用，水方块就会从静止状态垮塌下去，落地时粒子朝两侧飞溅，随后水体在箱底铺开、撞上右侧壁面、沿壁面爬升一段高度，再回落，最后在粘性耗散下慢慢接近静止水面。

这种场景在流体力学里有个专门名字叫溃坝问题（dam break），是SPH方法最常见的基准验证案例。它的好处是物理过程足够简单——不需要入口出口、不需要复杂的初始条件，重力、压力、粘性这三个力就能把整个流场的核心特征全部暴露出来。

代码里的可视化也很直白，每一帧都用 scatter 函数把粒子画成散点，颜色按粒子的速度大小或者压力大小来映射。速度大的区域用暖色，速度低的区域用冷色，一帧帧更新，看起来就像实实在在的水在动。我第一次跑通时盯着屏幕看了好几分钟，那种“一堆点竟然能呈现水的形态”的效果，比任何公式都更有说服力。

1.2 为什么选SPH而不是传统的网格法

你可能会问：做流体模拟，为什么不去用有限体积法、有限差分法这类经典CFD方法，偏偏用SPH？这里的关键差异在于坐标系视角。

传统网格法用的是欧拉视角。流体不动，网格固定，流体从网格里流过去。你要算出每个网格交点上的速度、压力，然后追踪自由表面的位置——这需要额外的界面处理方法，比如VOF（流体体积法）或者level set。一旦水面破碎、卷曲、某个液滴飞出去，网格法要处理的东西就非常麻烦，光是一个自由表面重构就够写好几篇论文。

SPH用的是拉格朗日视角。把流体离散成一个个粒子，这些粒子随身携带着速度、密度等物理量，在力的作用下各自运动，不需要网格，也不需要专门去追踪界面。水和空气的边界就是粒子群最外围那一圈，飞溅出来的一团粒子天然就是水滴。

我用一个类比来解释：网格法就像在一张固定坐标纸上画河流的流速分布，你得不断重新描绘河的边界；SPH则像是往水里撒一把锯末，锯末漂到哪里，水就流到哪里。对于有大量自由表面变形、飞溅、破碎的场景，SPH这种“跟着水走”的思路确实天生占优势。

所以这套代码的应用场景，从一开始就不是做那种高精度的工程管道流——它更适合水面波动、溃坝、液滴碰撞、液体搅拌这类“看得见的大变形”问题。Matlab做这些事情，恰好能兼顾算法实现和可视化验证。

2. SPH的核心原理：粒子如何“感受到”彼此

2.1 每个粒子身上都带什么“行李”

在SPH世界里，流体被拆成N个粒子，每个粒子不是一颗硬球，而是一个携带物理量的小“包裹”。具体来说，代码里的每个粒子都有这些属性：质量m、位置向量r、速度向量v、密度ρ、压力p。

质量在计算过程中通常保持不变，这是SPH的优势之一——质量守恒自动满足。位置和速度随时间更新。密度不是直接算出来的，而是通过对周围粒子的质量做加权求和得到的。压力则由密度通过状态方程换算。

那“光滑”体现在哪里呢？一个粒子的物理量不是只属于自己，而是通过一个核函数，像墨水滴到宣纸上一样“晕染”到周围一定范围内的空间去。这个范围由光滑长度h控制，一般取初始粒子间距的1.2到1.6倍。h太小，粒子只能感受到极少邻居，计算不稳定；h太大，粒子之间会“糊”在一起，丢失细节。

核函数本身在二维里用的是三次样条核（Cubic Spline Kernel）：

W(q) = (10 / (7πh²)) × { 1 - 1.5q² + 0.75q³, 0 ≤ q ≤ 1
0.25(2 - q)³, 1 < q ≤ 2
0, q > 2 }

这里的q就是粒子距离r除以光滑长度h。这个函数在q=0处最大，随距离增加先平缓下降、再加速衰减，到2h处刚好截断为零。好处是清晰的“紧支撑”——每个粒子只会影响2h范围内的邻居，不需要全局计算。

2.2 密度和压力：流体不可压缩是怎么逼近的

有了核函数，密度求和公式就顺理成章地写出来了：

ρ_i = Σ_j m_j W(r_ij, h)

意思就是，想知道第i个粒子的密度，就把周围每个邻居粒子j的质量乘以核函数权重，然后全加起来。直观理解是：你周围的粒子越多、越密，你的密度就越大。如果一个粒子附近几乎没邻居，比如飞溅出来的孤立液滴，它的密度会偏低，进而压力偏小，这符合物理直觉。

压力用Tait状态方程来算：

p_i = B × ((ρ_i / ρ_0)^γ - 1)

其中ρ_0是参考密度（初始设置的密度），γ一般取7，B是刚度系数。B越大，“反抗压缩”的能力越强，流体也就表现得越“硬”。为什么不用严格的不可压缩条件？因为在纯拉格朗日粒子法里，强制满足速度散度为零需要求解一个全局压力泊松方程，计算量和实现复杂度都暴增。弱可压缩SPH的妥协方法是：允许流体有微小可压缩性（通常密度变化不超过1%~3%），只要B取得足够大，视觉和工程精度上都足够以假乱真。

2.3 粒子受力：压力差、粘性、重力

粒子要动起来，就得受力。这套代码里，作用在每个粒子上的力分三块。

第一块是压力梯度力。流体从高压区向低压区流动，数学上就是压力梯度。粒子法里由离散求和实现：

F_pres_i = -Σ_j m_j × (p_i / ρ_i² + p_j / ρ_j²) × ∇W(r_ij, h)

这个式子看起来不直观，但要注意到压力项写成了“对称形式”——交换i和j之后符号相反，这保证了牛顿第三定律在离散层面也成立，一双粒子之间的压力相互抵消，总动量不会被自己内力破坏。这一点在写代码时极其重要，没有对称化的压力项，粒子系统自己就会慢慢产生非物理的净动量漂移。

第二块是粘性力，来自粒子间的相对运动。速度差越大、隔得越近，粘性力越强：

F_visc_i = Σ_j m_j × (μ × (v_j - v_i) / ρ_j) × ∇²W(r_ij, h)

这里μ是动力学粘性系数。粘性力相当于流体内部的“刹车系统”，能耗散掉多余动能，也是让水体最终能平静下来的关键。

第三块最简单，重力。每个粒子直接加速度g=9.81，方向向下。在2D模拟里，重力是驱动溃坝水柱倒塌的原始动力。

三个力算完，除以粒子自己的质量，就得到加速度a_i。然后就可以进入时间积分：更新速度、再更新位置。

2.4 时间推进：让画面动起来的节奏器

时间推进在代码里用的是类似蛙跳（leapfrog）的半隐式方法。基本循环是：

v_new = v_old + a × dt
r_new = r_old + v_new × dt

这个顺序看起来简单到无聊，但它是SPH里极少数能兼顾稳定性与实现简单的积分方式。你要是用前向欧拉把加速度先更新位置再更新速度，经常会看到粒子越飞越高、完全失控。

时间步长dt的选择是整个项目最关键、也最容易翻车的地方。CFL稳定条件告诉我们要满足：

dt ≤ 0.25 × h / c_max

其中c_max可以估算成最大声速或者最大粒子速度。初始粒子间距0.01米、声速取50到100米每秒量级时，dt量级在1e-4到5e-4秒之间。我知道很多人懒得做这个估计、随手填一个dt=0.001，结果跑20帧就看到粒子群“炸”成一团光斑。这不是代码错，是时间步长太大，数值稳定性没守住。

另外，为了让粒子不穿透水箱壁面，代码里还要加边界处理。最简单实用的方式是做墙体反射：当粒子坐标越过边界时，把位置拉回到边界内，同时把垂直于边界的速度分量反向并乘以一个恢复系数（比如0.5），模拟碰撞损耗。更精细的做法是放一排固定不动的墙体粒子，用它们之间的排斥力挡住内部粒子，但在入门代码里反射模型已经足够用了。

3. 实操：跑通代码，并让画面动起来

3.1 动手前先理清代码结构

这套Matlab代码不是一坨乱七八糟的脚本，而是按功能拆开的模块。我建议你拿到代码之后，先按下面这个结构在脑子里建立地图：

main.m是整个模拟的总调度。你打开它之后会发现第一段全是参数：粒子间距dp、粒子总数N、水箱长宽、参考密度rho0、光滑长度h、刚度系数B、粘性系数mu、重力加速度g、时间步长dt、总模拟时长tEnd。

我个人非常建议把初始化独立成一个函数，因为每次调整粒子规模或者水柱位置，只需要改这一处，其他地方完全不用动。如果你一次模拟想看多组对比，甚至可以把initialize.m写成一个接收参数的接口，非常方便。

3.2 初始参数怎么定

下面是一组我实测能稳定跑通的参数，供你直接复制使用：

初始粒子的分布我强烈建议用规则网格排列，不要用rand随机撒点。规则网格能保证初始时刻每个粒子周围邻居数量一致，密度均匀，压力为零。你如果随机撒粒子，初始密度就有高低起伏，程序一开始就背着巨大的压力修正工作，不仅慢，还容易在早期阶段产生虚假的密度波动，表现为粒子群边缘脱落、内部出现空洞。

具体做法很简单：取水柱区域矩形[xStart, xEnd] × [yStart, yEnd]，在两个方向上按dp间隔循环布置粒子，得到坐标矩阵后用meshgrid生成初始坐标。

3.3 邻居搜索与向量化：从暴力到网格

写SPH最直观的方式就是双重循环：对每个粒子i，遍历所有粒子j，判断是否满足r_ij < 2h，满足就算邻居。但当粒子数到4000以上时，双重循环意味着上千万次距离判断，Matlab里一帧就要算几秒钟，整个模拟跑下来会非常煎熬。

第一次实现我建议就用暴力法，因为代码简单、不容易写错，适合验证结果的正确性。等确认算法没问题了，再优化成网格邻居搜索。

网格搜索的核心思路是把模拟区域划分成边长约2h的方形格子。用坐标除以格子边长就能确定每个粒子落在哪个格子，然后只需要查当前格子及周围相邻格子里的粒子就行。在2D场景里，一个粒子最多只需要检查附近3×3共9个格子的粒子，计算量从O(N²)降到O(N×k)，k是平均邻居数——通常只有几十个。

Matlab里实现网格搜索有个很实用的函数accumarray，可以按格子编号把粒子索引分组。计算时的伪代码思路是：

% 将每个粒子放入格子编号 cellX = floor((x - xmin) / cellSize) + 1; cellY = floor((y - ymin) / cellSize) + 1; cellID = cellX + (cellY - 1) * nx; % 用accumarray或循环，把同格子粒子索引存进cell数组

然后主循环里对每个粒子，只遍历自身格子以及上下左右8个邻格里的粒子。实测中，即使只是简单写一下这个优化，同样的粒子规模速度能提升十几倍，动画流畅度完全不一样。

3.4 让动画流畅输出

Matlab里做实时动画最直观的写法是drawnow，但每帧都drawnow会带来很大的渲染开销。我建议你把绘图句柄提到主循环外，每帧只更新XData、YData和CData，而不是重新建一个scatter对象：

hScatter = scatter(x, y, 8, u, 'filled'); axis equal; xlim([0, Lx]); ylim([0, Ly]); for t = 1 : numSteps % 计算密度、压力、力、更新粒子 % ... set(hScatter, 'XData', x, 'YData', y, 'CData', speed); drawnow limitrate; end

使用drawnow limitrate可以避免画面更新过快拖垮CPU，还能让循环持续运行。如果你想把结果输出成GIF，就在循环里把当前帧的图片数据用getframe抓下来，再通过imwrite累加写入：

frame = getframe(gcf); [A, map] = rgb2ind(frame2im(frame), 256); if t == 1 imwrite(A, map, 'sph_result.gif', 'Loopcount', inf, 'DelayTime', 0.03); else imwrite(A, map, 'sph_result.gif', 'WriteMode', 'append', 'DelayTime', 0.03); end

不过注意，GIF写入本身比较慢，建议每隔10~20个时间步抓一帧，别每步都抓。

4. 常见问题与参数调试实录

4.1 粒子爆炸飞溅，像放烟花一样

这是SPH初学者第一个遇到的“灵异事件”。画面初始化看着没毛病，结果跑了没多少时间步，粒子就朝四面八方飞速飞出，速度大到离谱，甚至直接穿透水箱边界消失。

绝大多数情况下，原因是时间步长dt太大。SPH的力计算是非线性的，粒子一旦靠近，压力会急剧增大，如果dt不够小，积分就会像汽车在高速过弯时不减速——“甩尾”飞出。我的排查顺序是：先把dt缩小到原来的四分之一，比如从2e-4改成5e-5，看是否还爆炸；如果还不稳定，再查刚度系数B，B设得太大而dt又没相应缩小，也容易爆炸。

另外一个隐蔽原因：初始粒子有重叠。你如果手动修改初始区域时没调好循环范围，可能出现两个粒子距离极近，密度求出来异常大，压力巨高，直接互相弹射出去。检查方式很简单，初始化完成后画一下粒子位置，用矩形缩放看看有没有距离明显过近的点。

4.2 水体像“果冻”一样晃来晃去，波传播不自然

粒子不爆炸了，但水的行为怪怪的：整块水体像激光枪打中果冻，高频率微颤，或者波动传播得像慢动作，一点都不“水”。

这通常是因为刚度系数B太小。B本质上是等效声速的平方除以γ，B太小意味着流体“很软”，可以随意压缩，密度变化达到5%以上，压力波传播速度慢。把B调大到3到5倍，再配合合适的dt，水就会表现得“硬”起来。

但如果B已经很大还是像果冻，那多半是粘性系数mu太高了。尤其是很多人习惯把mu设成0.1以上，结果流体表现得像蜂蜜一样黏稠。水的实际动力粘性大约是1e-3量级，数值模拟中通常要加一点人工粘度来维持稳定，但不要超过0.05，否则你看到的就是一坨缓慢蠕动的胶水。

4.3 粒子从边界“渗漏”出去，水越跑越少

水箱是代码里用四个反射条件实现的，但如果处理得不好，粒子冲到墙边时会被困住，或者高速粒子直接“穿墙”。穿墙的根源是你的反射条件写在了速度更新之后，但是位置却更新在了速度前面——顺序弄反了。

正确的顺序是：先用蛙跳法得到新位置r_new，然后立刻检测边界，如果不满足水箱范围，就把坐标投影回去，再反转对应方向速度分量。而且反射检测要同时更新速度和位置，不能先算位置、再判断、却不改速度。调试时可以单独加一行输出，把每一帧穿出边界的粒子数量打印出来——这个数值应该始终为零。

还有一个更隐蔽的问题是“角落泄漏”。在角落处，粒子会同时越过x和y两个方向的边界。如果你用if-else而不是两个独立的if判断，角落粒子可能只处理了一个维度，另一个维度还留在墙外。务必写成两个独立的判断，一个对x方向、一个对y方向，各管各的反射。

4.4 跑得慢，等得慌

如果程序主循环跑一帧要好几十秒，说明还在暴力双重循环的坑里。我建议至少做两步优化：第一，把内层循环向量化，减少for循环次数；第二，实施网格邻居搜索。

先上向量化：在内层循环里，用数组运算一次算完该粒子到所有其他粒子的距离：

for recovery = 1:N dx = x_all - x_i; dy = y_all - y_i; r2 = dx.^2 + dy.^2; q = sqrt(r2) / h; idx = q < 2; % 只处理idx中为真的邻居 end

再把这种“每个粒子内层用全量数组”的方法改成网格邻域，速度就有明显提升。另外注意提前给所有数组分配内存，不要在时间循环里用不断增加元素的数组，避免Matlab反复扩容。用zeros预先分配好粒子数量N个行或者列的数组是基本素养。

4.5 每个症状的速查表

5. 从2D SPH出发还能往哪扩展

5.1 让流体更像真的：表面张力与多相流

基础SPH模拟里没有表面张力，所以水体看起来像“液体形状的软体”，不会自发地形成圆润液滴。想要更逼真的效果，可以引入表面张力模型：先通过粒子分布估算自由表面曲率，再用表面张力系数把沿法向的力加到表面粒子之上。这会让代码的效果有质的提升——水柱撞击地面后不再是直接散成一摊，而是能拉出细丝、甩出液滴，非常漂亮。

另一个自然的扩展是多相流。比如往水里放一个密度较小的气泡或者油滴。你只需要给不同粒子组设置不同的参考密度和粘性系数，再用相同的高斯核插值即可。2D代码里改起来很简单：给每个粒子加一个type或phase字段，在密度和受力计算时区分处理。

5.2 从2D到3D：难度在哪里

很多人跑完2D之后最想做的事就是直接改造3D。但3D不是简单地复制一遍公式就好，难度陡增三座大山：核函数归一化常数不同，2D是10/(7πh²)，3D变成了1/(πh³)，且连续导数形式也要跟着改；邻居数暴增，2D一个粒子五六十个邻居，3D直接两三百个起，网格搜索的格子也从9个变成27个，性能压力完全不是一个量级；可视化渲染不再能用scatter糊弄，得用particle plot、体积渲染或者ParaView导出来看，工作量很大。

我的建议是，除非毕业设计要求必须做3D，否则先用2D把SPH的物理和数值手感练熟，再考虑升维。2D里的大部分代码结构在3D中都能复用，但不用抱“改几个常量就完事”的幻想。

5.3 让结果更科学：数据导出与ParaView联动

Matlab的scatter画出来很直观，但真要写论文、做对比验证，还是需要专业后处理。最实用的路径是把粒子数据导出成VTK格式，然后用ParaView读取。

VTK导出很简单：把每个粒子的x、y、z（z全设为0）、速度、压力、密度写进固定的文本格式就行。ParaView打开之后，可以调整颜色映射、做切片、追踪某个粒子的轨迹，甚至录制动画。这套2D代码在讲解SPH到ParaView流水线时，比直接上3D反而更清爽，你不需要面对满屏粒子手忙脚乱。

如果要做定量验证，最经典的是对比水柱坍塌前沿位置随时间的变化曲线，这个结果与文献中的实验数据对照一下，能立刻看出自己的粘性参数设得对不对。

6. 我在实际调参过程中踩过的几个坑

最后分享一下我个人的几个小经验，这些都是常规文档里不会专门写给你的。

第一，要多打印状态量，不要只会画图。我在main循环里加了三个监控：总动能的均值、最大速度、穿出边界的粒子数。当模拟跑得不对劲时，这三个数字能快速定位问题——最大速度突然飙升，说明大约在爆炸边缘；穿界粒子数不为零，说明边界反射有问题；总动能单调增长，说明数值不守恒，多半是压力项不对称写错了。

第二，参数调节要有“复盘”的习惯。每次修改参数都记录一下效果，dt、B、h、mu这几个量对结果的影响是耦合的，只靠脑子记很容易混乱。我自己后来把参数组直接在代码注释里分版本保存，比如“v3：mu调低到0.01，dt缩小到1e-4，水柱倒塌轨迹明显更干净”。

第三，别焦虑粒子数量。很多人一上来就把粒子搞到一万两万，觉得“越多越真实”。但2D SPH里，从两千个加到八千个，效果提升有限，计算时间却翻好几倍。入门阶段，3000个粒子足够把溃坝的形态跑到像模像样。想要更光滑的水面，优先考虑把光滑长度h调大一点，或者换更高阶的核函数，而不是一味加粒子数。

我自己的体会是，SPH的迷人之处，在于你亲手把一堆数学公式变成动态的、有真实物理反馈的画面。刚开始跑那段代码时，我花了大半个晚上调dt、调边界、调配色，最后看着水柱砸向侧壁又反卷回来，那种成就感确实非常足。如果你也正在折腾这套2D SPH的Matlab代码，希望这篇分享能帮你少走一点弯路，早点看到自己那堵漂亮的“水墙”在屏幕上裂开。