流体模拟进阶实战用Python实现基于粒子方法的二维Navier-Stokes方程求解在计算机图形学、物理引擎和科学计算中流体模拟一直是极具挑战性和应用价值的方向。本文将带你从零开始构建一个基于粒子方法SPH, Smoothed Particle Hydrodynamics的二维流体模拟器使用Python NumPy Matplotlib实现核心算法并通过可视化展示流动过程。 核心思想SPH方法简介SPH是一种无网格数值方法适用于自由表面流体模拟。它将连续介质离散为一系列粒子每个粒子携带密度、速度、压力等信息通过核函数进行局部插值来近似偏微分方程。关键公式简化版密度估计$$\rho_i \sum_j m_j W(|\mathbf{r}_i - \mathbf{r}_j|, h)$$压力力计算内聚力$$\mathbf{F}{\text{pressure}, i} -\sum_j m_j \left( \frac{p_i}{\rho_i^2} \frac{p_j}{\rho_j^2} \right) \nabla W{ij}$$粘性力项人工粘性$$\mathbf{F}_{\text{visc}, i} \mu \sum_j m_j \frac{\mathbf{v}_j - \mathbf{v}i}{\rho_j} \nabla W{ij}$$其中WWW是平滑核函数常用的是Gaussian或Cubic Splinehhh是光滑长度μ\muμ是粘度系数。 代码实现基础框架搭建importnumpyasnpimportmatplotlib.pyplotaspltfrommatplotlib.animationimportfuncAnimation# 参数设置N1000# 粒子数dt0.01# 时间步长h0.05# 光滑长度mu0.1# 粘度系数gravitynp.array([0,-9.8])# 重力加速度# 初始化粒子位置与属性posnp.random.rand(N,2)*0.80.1# 初始均匀分布于[0.1, 0.9]²velnp.zeros_like(pos)mass0.001densitynp.ones(N)*1000# 初始密度一致pressurenp.zeros(N)defcompute_density(pos,mass,h):计算每个粒子的密度Nlen(pos)densitynp.zeros(N)foriinrange(N):rnp.linalg.norm(pos[i]-pos,axis1)weightsnp.where(rh,kernel(r,h),0)density[i]np.sum(mass*weights)returndensitydefkernel(r,h):Cubic spline核函数ifr0:return1/(np.pi*h**2)elifrh:return15/(7*np.pi*h**2)*(1-1.5*(r/h)**20.75*(r/h)**3)else:return0defcompute_pressure_force(pos,density,pressure,mass,h):计算压力力Nlen(pos)forcenp.zeros_like(pos)foriinrange(N):rnp.linalg.norm(pos[i]-pos,axis1)grad_kernelgradient_kernel(r,pos[i],pos,h)p_term-(pressure[i]/density[i]**2pressure/density**2)force[i]np.sum(mass*p_term[:,None]*grad_kernel,axis0)returnforcedefgradient_kernel(r,xi,xj,h):梯度核函数drxj-xi normnp.linalg.norm9dr,axis1)masknormh gradnp.zeros_like(dr)grad[mask]-15/(7*np.pi*h**3)*dr[mask]8(1-1.5*(norm[mask]/h)**2)returngraddefupdate_velocity(pos,vel,force,dt,mass):更新速度accforce/mass velacc*dtreturnveldefintegrate(pos,vel,dt):积分更新位置posvel*dt# 边界处理简单反射pos[pos0]-pos[pos0]pos[pos1]2-pos[pos1]# 主循环fig,axplt.subplots(figsize(8,8))scatterax.scatter([],[],cblue,s10)defanimate(frame):globalpos,vel,density,pressure# 步骤1: 更新密度densitycompute_density(pos,mass,h)# 步骤2: 计算压力力pressure100*density# 压力正比于密度简化模型pressure_forcecompute_pressure_force(pos,density,pressure,mass,h)# 步骤3: 添加粘性力简略版visc_forcenp.zeros_like(pos)foriinrange(len(pos)):neighborsnp.where(np.linalg.norm(pos-pos[i],axis1)h)[0]iflen(neighbors)1:avg_velnp.mean(vel[neighbors],axis0)visc_force[i]mu*(avg_vel-vel[i])# 总力 压力 粘性 重力total_forcepressure_forcevisc_forcegravity*mass# 更新速度和位置velupdate-velocity(pos,vel,total_force,dt,mass0 integrate(pos,vel,dt0 scatter.set_offsets(pos)ax.set_xlim(0,1)ax.set_ylim(0,1)ax.set_aspect(equal)ax.set_title(fFrame{frame})returnscatter,aniFuncAnimation(fig,animate,frames1000,interval50,blitTrue)plt.show()3## 可视化效果说明运行上述脚本后你会看到一个动态的粒子云在重力作用下下落并相互挤压、扩散的过程这就是最基础的流体行为模拟扩展建议加入边界碰撞检测矩形墙引入表面张力增加粒子间吸引力使用更复杂的核函数如 Wendland C2转换为GPU加速版本可用PyCUDA或JAX 模拟流程图伪代码逻辑示意初始化粒子状态 → 计算当前密度 → 计算压力力 粘性力 → 更新速度 → 积分位置 → 边界检查 → 绘制帧 → 循环直到结束✅ 图中标注清晰适合嵌入到文章中作为技术路线图。 小结为何选择这种方案相比传统的格点法如Finite Difference/VolumeSPH具有以下优势自适应分辨率无需固定网格自然处理自由表面比如水滴飞溅易于并行化尤其适合现代GPU架构开源生态丰富可结合Blender、unity做动画这正是我们在游戏开发、影视特效甚至气象模拟中广泛采用的原因。 如果你是刚接触流体力学的程序员不妨先跑通这个例子在此基础上逐步加入复杂特性。你会发现编程不仅是工具更是探索物理世界的新语言