✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。
🔥 内容介绍
一、引言：流体模拟的 “粒子革命”——2D SPH 为何出圈？
1.1 从天体物理到工程实验室：SPH 的跨界之旅
回溯到 1977 年，当科学家们试图揭开宇宙中星体碰撞的神秘面纱时，光滑粒子流体动力学（Smoothed Particle Hydrodynamics，简称 SPH）方法应运而生 。在传统的流体模拟中，网格就像是一张密织的渔网，将流体空间划分成规整的小块，通过在这些网格节点上求解流体力学方程来模拟流体行为。但在面对星体这种尺度巨大、形态多变且运动复杂的对象时，传统网格法显得力不从心。SPH 方法则另辟蹊径，它摒弃了固定网格的束缚，把流体看作是由无数个携带质量、速度等物理量的粒子组成。每个粒子就像是一个微小的 “使者”，通过与周围粒子的相互作用，来传递和反映流体的物理特性，如同在浩瀚星空中，每颗星星都遵循着引力法则相互影响，共同塑造出壮观的宇宙动态。
随着技术的不断演进，云境智仿与 Peridyno 引擎的耦合，让 SPH 从晦涩的学术论文中走进了大众视野。曾经只有科研机构才能触及的高端模拟技术，如今在云端就能轻松实现，工程师们可以在虚拟环境中利用 SPH 模拟汽车碰撞时的流体飞溅、船舶航行时的水流阻力 ；而爱好者们也能借助它创造出美轮美奂的流体艺术效果，比如梦幻的水花绽放、灵动的液体流动。在众多 SPH 应用场景中，2D SPH 以其较低的计算量和易于可视化的特点，成为了人们踏入流体模拟领域的最佳起点，就像一把小巧的钥匙，开启了通往奇妙流体世界的大门。
1.2 传统网格法的 “痛点”：SPH 无网格技术的核心优势
传统网格法在处理简单流体问题时，就像一位训练有素的工匠，能够有条不紊地完成任务。它通过在规则网格上离散化流体方程，精确地计算每个网格点上的物理量，对于那些边界规则、流体运动平稳的场景，比如管道内的水流，传统网格法可以给出相当准确的模拟结果。但当面对复杂边界和大变形流体时，传统网格法的局限性便暴露无遗。想象一下，当模拟海浪冲击崎岖的海岸线，或者液滴在高速撞击下破碎飞溅的场景，网格就像被强行扭曲的拼图，难以准确贴合流体的复杂形态。随着流体的剧烈运动，网格可能会严重扭曲甚至撕裂，导致计算精度大幅下降，就如同用一张破损的地图去导航，最终迷失方向。
而 SPH 作为纯 Lagrange 方法，就像是一群自由灵动的舞者，不受固定框架的约束。在 SPH 的世界里，粒子才是主角，它们自由地穿梭、相互作用，每个粒子都带着独特的 “使命”—— 质量、速度、压力等物理属性，通过核函数的 “桥梁”，与周围粒子交换信息，共同演绎出流体的运动轨迹。尤其是在 2D 场景下，这种优势更加凸显。较低的维度使得计算复杂度降低，我们可以更快速地搭建模型，验证各种假设。同时，2D 画面的直观性让我们能清晰地观察到流体自由表面的流动细节，比如涟漪的扩散、液膜的破裂，这些细微而美妙的变化在 2D SPH 模拟中都能被生动地呈现出来。
1.3 本文导航：从原理吃透到亲手做 2D SPH 模拟
在接下来的内容中，我们将深入剖析 2D SPH 的核心原理，从基础的数学公式到物理意义的解读，一步一步揭开它神秘的面纱，让你不仅知其然，更知其所以然。随后，通过经典案例分析，带你领略 2D SPH 在不同场景下的精彩应用，感受它强大的模拟能力。还会为你带来开源工具的实操教程，手把手教你如何利用这些免费又强大的资源，搭建自己的 2D SPH 模拟环境，实现从理论到实践的跨越。最后，分享一些进阶技巧，帮助你优化模拟效果，提升计算效率，让你的 2D SPH 模拟更上一层楼。无论你是对流体模拟充满好奇的小白，还是渴望深入探索 SPH 技术的进阶者，本文都将成为你掌握 2D SPH 模拟的得力指南 ，带你开启一段充满惊喜与挑战的流体模拟之旅。
二、核心解密：2D SPH 模拟的底层逻辑
2.1 无网格的本质：流体的 “粒子化离散”
在 2D SPH 模拟中，我们彻底告别了传统网格的限制，将流体看作是由大量离散粒子组成的集合 。就好比把一片平静的湖面想象成无数颗闪闪发光的水珠，每颗水珠都是一个独立的粒子，它们各自携带质量、密度、速度等关键物理属性，在二维平面上自由分布、相互作用，共同演绎出湖水的流动、涟漪的产生等各种奇妙现象。
当风吹过湖面，这些粒子就开始运动，它们的速度和方向不断改变，而这种改变正是通过与周围粒子的相互作用实现的。我们通过求解粒子动力学方程，就像给每颗水珠设定了独特的 “行动指南”，来精确跟踪每个粒子的运动轨迹，从而在宏观上还原出流体的整体行为。这种粒子化离散的方式，让我们能够轻松应对流体的复杂变形和自由表面运动，就像一位技艺高超的舞者，可以在舞台上自由地旋转、跳跃，展现出无限的灵动之美。
在 2D 场景下，粒子仅在平面内分布，无需考虑 Z 轴维度，这就好比我们只在一张纸上作画，不用考虑画面的立体深度。这种简化使得计算复杂度大幅降低，就像从解一道复杂的多元方程变成了解简单的一元方程，计算资源的消耗也大大减少。同时，与传统网格法相比，SPH 对质点排列没有严格要求，粒子们可以自由地 “散落” 在平面上，更加贴合流体的自然特性，也为模拟带来了更大的灵活性。
2.2 关键参数：核函数与光滑长度的 “魔法”
2.2.1 核函数：粒子间相互作用的 “纽带”
核函数在 2D SPH 模拟中扮演着至关重要的角色，它就像是粒子间相互作用的 “纽带”，定义了粒子间相互作用的权重 。简单来说，当两个粒子距离较近时，它们之间的相互作用就越强，权重也就越高；反之，距离较远时，相互作用就越弱，权重越低。
以常见的 Poly6 核函数为例，它的形式在 2D 场景下与 3D 有所不同，但其核心思想都是通过数学公式来描述这种距离与权重的关系 。Poly6 核函数常用于密度计算，它能够根据周围粒子的分布情况，准确地计算出每个粒子处的密度值，就像一个精准的 “密度探测器”。而 Spiky 函数则主要用于压力计算，它通过巧妙的数学构造，使得压力在粒子间的传递更加合理，保证了模拟中压力分布的稳定性，就像给流体的压力传递搭建了一座坚固的 “桥梁”。
不同的核函数对模拟精度有着显著的影响。在压力计算中，如果选择的核函数不合适，可能会导致压力震荡，使得模拟结果出现不稳定的情况，就像一座摇晃的桥梁，无法承载流体的真实运动。而合适的核函数则能让压力计算更加稳定，模拟结果更加接近真实的流体行为，让我们看到的流体动画更加流畅、自然。
2.2.2 光滑长度 h：控制粒子影响范围的 “开关”
光滑长度 h 是核函数的关键参数，它就像是一个控制粒子影响范围的 “开关”，决定了粒子间相互作用的有效范围 。在 2D 模拟中，我们可以把每个粒子想象成一个拥有 “势力范围” 的小团体，这个势力范围的大小就是由光滑长度 h 决定的。当一个粒子在运动时，它只会与处于自己光滑长度范围内的其他粒子发生相互作用，交换物理信息，就像在一个社交圈子里，人们只会与自己身边的人交流互动。
光滑长度的选取至关重要。如果 h 选取过大，粒子的影响范围就会过大，导致模拟结果变得模糊，丢失很多细节信息，就像用一个超大号的画笔去描绘精细的图案，只能得到一个粗糙的轮廓；反之，如果 h 选取过小，每个粒子周围的邻居粒子数量就会不足，这会导致数值计算出现震荡，模拟结果变得不稳定，就像一个团队成员太少，无法有效地完成任务。
为了提升模拟精度，我们可以根据粒子密度动态调整光滑长度 。在粒子密度较高的区域，适当减小光滑长度，让粒子间的相互作用更加精细；在粒子密度较低的区域，增大光滑长度，保证每个粒子都能与足够的邻居粒子相互作用，就像在热闹的集市和冷清的小巷中，采用不同的社交策略，以达到最佳的交流效果。
2.3 2D SPH 的简化方程：从连续到离散的数学转换
2D SPH 模拟的背后，是一系列从连续流体力学方程到离散化形式的巧妙数学转换 。我们从流体力学的基本方程出发，这些方程描述了流体在连续介质中的运动规律，就像描述一条河流在大地上的流淌。但在 SPH 中，我们要将其转化为适合粒子描述的形式。
以连续性方程为例，在连续介质中，它描述了流体密度随时间和空间的变化关系。而在 2D SPH 中，我们通过对周围粒子的质量和分布进行求和，将其离散化，用来更新每个粒子的密度，就像把连续的水流分割成一个个小水滴，通过计算每个小水滴周围的水滴分布来确定它的密度变化。
动量方程则负责计算粒子所受的力，包括压力梯度力和粘性力等 。在 2D 场景下，压力梯度力通过对周围粒子的压力和位置进行计算得到，粘性力则用于模拟流体的内摩擦力，它们共同决定了粒子的加速度和运动方向，就像一群舞者在舞台上，根据周围同伴的位置和力量，调整自己的动作和方向。
与 3D SPH 方程相比，2D SPH 方程的积分维度降低，这大大简化了计算过程，就像从三维空间的复杂舞蹈变成了二维平面的简单舞步；边界条件也相对简化，我们无需考虑复杂的三维边界情况，让模拟更加容易实现 。这些方程虽然看起来有些复杂，但它们背后的物理意义却很直观，它们是我们理解和模拟流体运动的有力工具，就像地图是我们探索未知世界的指南。
2.4 2D vs 3D SPH：入门首选的核心原因
对于初学者来说，2D SPH 无疑是踏入流体模拟领域的最佳选择 。首先，它的计算资源消耗低，普通的个人电脑就能轻松运行模拟，不像 3D SPH 需要强大的计算集群，这就像骑自行车和开跑车，前者更容易上手。2D 模拟的可视化效果非常直观，我们可以在平面上清晰地看到流体的运动轨迹、形态变化，比如水波的荡漾、液滴的融合，就像看一场平面动画，所有细节尽收眼底。
2D SPH 的模型调试周期短 。在开发模拟程序时，我们可以快速地修改参数、调整模型，然后迅速得到结果反馈，不断优化模拟效果，就像在一个小实验室里做实验，能够快速验证各种想法。当然，2D SPH 也有其局限性，它无法模拟三维空间中复杂的流场，比如龙卷风的三维旋转、深海中复杂的水流运动。但这并不影响它作为入门工具的地位，通过学习 2D SPH，我们可以深入理解 SPH 的基本原理和方法，为进一步探索 3D SPH 打下坚实的基础，就像先学会在浅水区游泳，再挑战深水区的巨浪。
⛳️ 运行结果
📣 部分代码
talNumParticles) = 1000;
%Loops through all of the fluid particles
for i = 1:numParticles
xAcc_i = 0; % Reset values for every particle
yAcc_i = 0;
rho_i = Rho_RhoHalf_dRho(1,i);
drho_i = 0;
p_i = Pressures(1, i);
x_i = Pos(1,i);
y_i = Pos(2,i);
vx_i = Vel(1,i);
vy_i = Vel(2,i);
Neighbor = NL{i,1};
[t, m] = size(Neighbor);
for k = 2:t
j = Neighbor(k);
%Loops through all of the neighbors to the fluid
% for j = 1:totalNumParticles
x_j = Pos(1,j);
y_j = Pos(2,j);
dx = x_i - x_j;
dy = y_i - y_j;
r_ij = [dx; dy];
norm_r_ij = norm(r_ij);
q = norm_r_ij/h;
% if ((i ~= j) && (q<=2) && (q>=0))
rho_j = Rho_RhoHalf_dRho(1,j);
p_j = Pressures(1, j);
vx_j = Vel(1,j);
vy_j = Vel(2,j);
dvx = vx_i - vx_j;
dvy = vy_i - vy_j;
v_ij = [dvx; dvy];
dw_ij = dW(q, h);
grad_a_wab = (dw_ij/(h2*q)) * r_ij;
rho_bar = (rho_i+rho_j)/2;
% Everything is a scalar except unit_r_ij
pressureTerm = (particleMass*((p_i/(rho_i*rho_i)) + (p_j/(rho_j*rho_j))))...
.* grad_a_wab;
muNumerator = particleMass * 2*mu .* (r_ij' * grad_a_wab) .* v_ij;
muDenominator = (rho_bar*rho_bar)*(norm_r_ij*norm_r_ij + h2*epsilon);
muTerm = (muNumerator./muDenominator);
xAcc_i = xAcc_i - pressureTerm(1) + muTerm(1);
yAcc_i = yAcc_i - pressureTerm(2) + muTerm(2);
drho_i = drho_i + (rho_i/rho_j)*particleMass*v_ij'*grad_a_wab;
% end
end
% Leap frog stuff:
% stepTerm = 0.5 for first step. Otherwise stepTerm==1.
% Acceleration
Acc(1,i) = xAcc_i;
Acc(2,i) = yAcc_i + g;
% Velocity at t+0.5*dt
VelHalf(1,i) = VelHalf(1,i) + stepTerm*dt*Acc(1,i);
VelHalf(2,i) = VelHalf(2,i) + stepTerm*dt*Acc(2,i);
% rho at t+0.5*dt
Rho_RhoHalf_dRhoNext(2,i) = Rho_RhoHalf_dRho(2,i) + stepTerm*dt*drho_i;
% dRho
Rho_RhoHalf_dRhoNext(3,i) = drho_i;
% Temporary Position
PosNext(1,i) = Pos(1,i) + dt * VelHalf(1,i);
PosNext(2,i) = Pos(2,i) + dt * VelHalf(2,i);
end
Pos = PosNext;
Rho_RhoHalf_dRho = Rho_RhoHalf_dRhoNext;
end
🔗 参考文献
🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：
🌈 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱调度、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划（2E-VRP）、充电车辆路径规划（EVRP）、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题、港口调度、港口岸桥调度、停机位分配、机场航班调度、泄漏源定位
🌈 机器学习和深度学习时序、回归、分类、聚类和降维
2.1 bp时序、回归预测和分类
2.2 ENS声神经网络时序、回归预测和分类
2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类
2.4 CNN|TCN|GCN卷积神经网络系列时序、回归预测和分类
2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类
2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类
2.7 ELMAN递归神经网络时序、回归\预测和分类
2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类
2.9 RBF径向基神经网络时序、回归预测和分类
