1. 这个项目到底想解决什么问题先说结论我这次折腾的目标是纯粹用 Python 从零实现一个能实时交互的 2D 流体模拟器核心基于 Navier-Stokes 方程最终通过 GPU 加速把计算规模拉到肉眼可见的流畅帧率整套流程覆盖了从物理建模、数值求解到可视化渲染的完整链路。写这篇文章前我刚把最后一版代码跑通256x256 网格在普通笔记本的集成显卡上能做到每秒 60 帧以上的更新速度拖动鼠标往画面里注入染料的时候那种“黏乎乎”的扩散感确实有那么点真实流体的意思了。为什么要在 Python 里做这件事而不是直接用 C 或引擎因为 Python 生态里能快速验证算法、调试参数的优势实在太明显了上手成本也低你可以先用小网格把数值方法跑明白再考虑性能优化。这篇文章适合三类人想入门流体模拟、但对 C 和图形学 APIs 望而却步的同学已经在做图形学相关方向、想快速搭一个可交互原型的研究者以及单纯想看看 Navier-Stokes 方程在一个普通电脑上能跑出什么效果的好奇党。1. 方案设计与技术路线1.1 流体模拟的本质把“看不见的力”变成“看得见的场”流体模拟在计算机图形学里核心不是去模拟真实世界的每一个分子而是用网格或粒子去近似描述速度场、密度场和压力场。你看到水面波纹、烟雾飘散、墨水在水中晕开本质上都是这些场在时间轴上不断演化的结果。2D 场景下我们把空间划分成一个个小格子每个格子里存着水平和垂直方向的速度分量加上一个标量的密度值代表染料浓度或者烟雾量。随着时间推进每个格子里的物理量按照流体运动规律更新然后渲染出来——就形成了流动的假象。在动手之前最需要想清楚的问题不是“怎么写代码”而是“怎么把一个连续的偏微分方程变成计算机能算的离散代扰”。这一步如果走偏后面所有优化都是在错误的地基上盖楼。我自己一开始就是直接拿简单的扩散公式去模拟烟雾结果发现烟雾会快速消失、边缘发虚完全没有旋涡和翻滚的细节——就是因为没把平流advection这个关键步骤做好。1.2 为什么选择 Python快速验证的价值远大于运行速度如果你去搜索流体模拟的高性能实现绝大多数结果都会指向 C、CUDA 或者引擎插件。确实传统观点认为实时流体模拟必须用底层语言才够快。但在实际做项目时我们必须区分“工程目标”和“学习目标”。如果目标是产出一个大型引擎里的高性能流体模块直接用 C 是合理的。但我的目标是在最短时间内从方程到可视化全流程跑通验证算法效果同时保留快速迭代参数的空间。Python 在这方面的优势是碾压级的——NumPy 的向量化操作让矩阵运算极其简洁安装依赖只需要一个 pip 命令调试时可以直接在脚本里打印任意中间变量配合 Jupyter 可以非常直观地观察每一步的数值变化。即使后面要做 GPU 加速也有 Taichi、Numba 这些工具让 Python 代码无缝跑在显卡上完全不需要写一行 CUDA。1.3 技术选型全景从纯 CPU 到 GPU 的三级跳我最终的技术路径分三步走第一步纯 NumPy SciPy 实现 2D 求解器网格 64x64跑在 CPU 上作为正确性验证和算法理解的基础版本。第二步用 Matplotlib 做动态渲染把密度场和速度场可视化成可交互的动态视频这阶段帧率只有几 fps但足够观察物理行为是否合理。第三步接入 Taichi把核心循环改写为 GPU 并行内核网格直接拉大到 256x256 甚至 512x512帧率提升到实时水平。这套路线的核心价值在于每一步都有明确的验证节点。如果你直接学“从零写 CUDA 流体模拟”调试一次数值爆炸要浪费大量时间但如果你先用 NumPy 把物理弄明白再用 Taichi 做性能替换代码迁移成本会低很多错误也更容易定位到具体环节。2. Navier-Stokes 方程从物理到可计算的数值方法2.1 把偏微分方程翻译成机械运动的直觉Navier-Stokes 方程看起来抽象但它的物理含义其实非常直觉化。以不可压缩流体的 2D 形式为例速度场方程∂u/∂t -(u·∇)u - (1/ρ)∇p ν∇²u F它描述的是某一点流体的速度随时间的变化取决于四个因素——对流项流体沿着自身速度方向运动把上游的速度“搬运”过来、压力梯度项流体总是从高压流向低压、粘性扩散项速度差异会被粘性抹平、以及外部力比如我们在屏幕上鼠标施加的力。密度方程∂ρ/∂t -(u·∇)ρ κ∇²ρ S它描述的是染料或烟雾浓度随时间的变化随着速度场被搬运对流、自然扩散、以及由外部源产生。如果你做过图像处理里的“光流”或“热扩散”会发现这些方程的结构非常相似——而流体模拟的高明之处正是把这几套机制耦合在一起互相驱动。2.2 数值求解的核心矛盾对流项的“不稳定性”从方程到代码之间最大的坎是对流项。如果用朴素的有限差分法直接计算“流体下一时刻的位置”本质上是在做显式时间积分每个网格点根据当前速度向外搬运数值。只要时间步长稍大或者速度较快就会出现“数值过冲”——某个格子的值被搬过头导致负密度、负速度然后整个系统直接爆炸。Jos Stam 在 1999 年提出了经典的“稳定流体”方案核心思路是把对流项改为“半拉格朗日法”不是从当前格点出发追踪流体去了哪里而是反向追踪——假定下一时刻某个位置的流体来自当前时刻的哪个位置然后把那个位置的值搬过来。这个思路从根本上保证了数值稳定性因为它在每一步做的是“插值采样”而不是“外推”。这个方法在图形学里被普遍使用也是我这次实现的主干。2.3 三个关键数值技巧交错网格、投影法、迭代求解解决好对流之后另一道坎是“不可压缩约束”。物理上水不能被压缩所以速度场的散度必须为零。算法上每个时间步计算出的速度场通常带着一定散度需要做一个“压力投影”操作来修正先解一个关于压力的泊松方程然后从速度场中减去压力的梯度把散度“挤”出去。操作中我采用了 MACMarker-and-Cell交错网格把速度分量 u 和 v 错开半个格子存储而不是都存放在同一个网格点上。这么做的好处是计算散度和压力梯度时避免了奇偶失联问题棋盘格状伪影很多初学者都会忽略这个细节直接用相同的网格分辨率存储所有变量结果发现总出现不太自然的小格子状噪声其实就是网格布局的问题。压力泊松方程的求解我用了标准的红黑 Gauss-Seidel 迭代。对于 64x64 的网格50 次迭代足以让残差降低到对人眼无感知的水平对于 256x256 网格在 GPU 上用加权 Jacobi 迭代配合多网格思路速度也很快。3. Python 实现核心求解器3.1 数据结构与代码框架先搭骨架再填肉我先把求解器的核心数据结构写出来这一步看起来不起眼但决定了后面所有代码的清晰度。import numpy as np class FluidSim: def __init__(self, N64, dt0.1, viscosity0.0, diffusion0.0): self.N N self.dt dt self.visc viscosity self.diff diffusion # 交错网格布局 # 速度 u 存放在 (i0.5, j) 处 # 速度 v 存放在 (i, j0.5) 处 # 标量密度 p 存放在 (i, j) 处 # 外圈多一层疙瘩便于实现边界条件 self.size N 2 self.u np.zeros((self.size, self.size), dtypenp.float32) self.v np.zeros((self.size, self.size), dtypenp.float32) self.u_prev np.zeros_like(self.u) self.v_prev np.zeros_like(self.v) self.dens np.zeros((self.size, self.size), dtypenp.float32) self.dens_prev np.zeros_like(self.dens)每个场都多开一圈“边界网格”这是为了在处理边界条件时不用写一层复杂的 if-else 分支——直接把邻域值拷贝到边界网格后续的所有差分计算就能统一用同一套索引公式。很多线上教程喜欢写简洁的一行式代码但实际工程上这种“空间换逻辑”的写法更好扩展也更不容易出边界 bug。3.2 时间步进四个阶段各司其职每一帧模拟我按顺序执行四个阶段施加力源、扩散、平流、投影。用一个 step 函数串联起来def step(self): N self.N dt self.dt # 1. 施加外部力例如鼠标拖动产生的力 self.add_source(self.u, self.u_prev) self.add_source(self.v, self.v_prev) self.add_source(self.dens, self.dens_prev) # 2. 交换时序缓冲 self.u_prev, self.u self.u, self.u_prev self.v_prev, self.v self.v, self.v_prev self.dens_prev, self.dens self.dens, self.dens_prev # 3. 扩散粘性和浓度扩散 self.diffuse(0, self.u, self.u_prev, self.visc) self.diffuse(0, self.v, self.v_prev, self.visc) self.diffuse(0, self.dens, self.dens_prev, self.diff) # 4. 平流半拉格朗日反向追踪 self.advect(0, self.u, self.u_prev, self.u_prev, self.v_prev) self.advect(0, self.v, self.v_prev, self.u_prev, self.v_prev) self.advect(1, self.dens, self.dens_prev, self.u, self.v) # 5. 投影使速度场满足不可压缩条件 self.project(self.u, self.v)这是个典型的稳定流体框架。有意思的是我在调参时发现如果粘度和扩散系数设成 0液体依然会保持稳定的对流效果只是不会自己减速——这正好符合“理想流体”的物理直觉也有利于观察湍流细节。注意这里我把“投影”放在平流之后执行能保证最终输出的速度场是无散的渲染出来的流线更加平滑。3.3 半拉格朗日平流的代码细节与坑平流是这套算法里最需要仔细看的部分。以速度场为例我们要计算新时刻每个网格点上的速度值做法是反向追踪这个点在一小段时间之前来自哪里def advect(self, bound, d, d0, u, v): N self.N dt self.dt for i in range(1, N 1): for j in range(1, N 1): # 当前网格点坐标归一化到 0~N x i y j # 时间倒推 px x - dt * u[i, j] * N py y - dt * v[i, j] * N # 截断到有效范围 px max(0.5, min(N 0.5, px)) py max(0.5, min(N 0.5, py)) # 双线性插值四周的网格点 i0 int(px) j0 int(py) i1 i0 1 j1 j0 1 sx px - i0 sy py - j0 d[i, j] (1 - sx) * ((1 - sy) * d0[i0, j0] sy * d0[i0, j1]) \ sx * ((1 - sy) * d0[i1, j0] sy * d0[i1, j1]) self.set_bnd(bound, d)这里我犯过的错误就是把 px、py 计算成数组标量之后在循环里混用。在纯 NumPy 版本里这类循环可以用 roll 和高级索引替换掉但为了可读性我建议初学者先用双循环把逻辑跑通再谈向量化。还有一个关键点反向追踪的速度必须用旧时刻的速度场而不是新时刻的否则会引入隐式的时间耦合导致数值发散。这个顺序问题在教科书里很少强调但我实际调试时发现一旦写反画面会立刻出现剧烈的随机抖动。3.4 投影步骤与边界条件隐形杀手的两次相遇投影步骤的目标是让速度场散度为零。先计算散度再解泊松方程最后修正速度def project(self, u, v): N self.N h 1.0 / N # 计算散度同时构造压力泊松方程的右端项 div np.zeros((self.size, self.size), dtypenp.float32) p np.zeros((self.size, self.size), dtypenp.float32) for i in range(1, N 1): for j in range(1, N 1): div[i, j] -0.5 * h * ( u[i1, j] - u[i-1, j] v[i, j1] - v[i, j-1] ) # 迭代求解压力场 for _ in range(20): for i in range(1, N 1): for j in range(1, N 1): p[i, j] (div[i, j] p[i-1, j] p[i1, j] p[i, j-1] p[i, j1]) / 4.0 self.set_bnd(0, p) # 用压力梯度修正速度场 for i in range(1, N 1): for j in range(1, N 1): u[i, j] - 0.5 * (p[i1, j] - p[i-1, j]) / h v[i, j] - 0.5 * (p[i, j1] - p[i, j-1]) / h self.set_bnd(2, u) self.set_bnd(3, v)这里的边界条件我分了四类密度场是“无通量边界”数值在边界上不流失梯度为零水平速度 u 在左右边界上取镜像反射、在上下边界上取反向垂直速度 v 则正好相反。这个细节初看很绕但如果边界条件写错模拟空间里会莫名出现“吸水的黑洞”或者“喷射的泉眼”在视觉上非常明显。4. 可视化与渲染让看不见的场变成看得见的艺术4.1 Matplotlib 动态渲染与交互控制求解器本身只是“一半”的作品可视化是让结果呈现出来的另一半。我用的基础渲染工具是 Matplotlib 的 imshow 和 FuncAnimation。核心思路是每一帧从求解器中读出密度场将其映射为颜色然后更新图像对象。import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def update(frame): sim.step() ax.clear() im ax.imshow(sim.dens.T, cmapinferno, originlower, vmin0, vmax1) return im, fig, ax plt.subplots(figsize(6, 6)) sim FluidSim(N64, dt0.1, diffusion0.0001) sim.dens_prev[:, :] 0 # 在场景中央注入一团“染料” sim.dens[::2, ::2] 0.5 anim FuncAnimation(fig, update, interval16, blitFalse) plt.show()有一个细节在更新时不要反复创建新的 im 对象否则动画会越来越卡更高级的做法是把 im 对象创建一次在 update 里用 set_data 更新数据数组。上边代码为了简洁用了 ax.clear()实际长时间运行时内存会缓慢增长。真要跑长时间动画建议改用 set_data 的方式。4.2 颜色映射与视觉增强技巧密度场的可视化可以直接用灰度图但视觉冲击力会差很多。我用了路径式颜色映射让低密度区域显得透明、高密度区域呈现亮色配合暗色背景后烟雾的层次感立刻出来了。我另外做了速度场的可视化直接用 streamplot 画流线。说实话速度场的流线比密度场更能体现流体运动的旋涡结构——你会看到两个反向旋涡之间形成的刺状流线这是从方程中直接“长”出来的结构而不是人为画的装饰。两者叠加显示时先用 imshow 画密度做底再用 streamplot 画速度流线做顶信息量非常丰富。4.3 输出视频文件的实操细节在 Jupyter 环境里直接 plt.show() 很方便但没法分享和保存。要输出 mp4 视频我建议用 ffmpeg 方式保存from matplotlib.animation import FFMpegWriter writer FFMpegWriter(fps60, bitrate-1, codeclibx264) anim.save(fluid_sim.mp4, writerwriter, dpi100)这里踩过一个坑如果系统里没有安装 ffmpeg或者 PATH 没配置好这个保存会直接报错加上文件路径包含中文还容易出编码问题。建议在项目目录下建立一个 bin 文件夹把 ffmpeg 可执行文件放进去然后在代码里显式指定plt.rcParams[animation.ffmpeg_path] ./bin/ffmpeg另外一个经验是fps 不要一味调高默认 30 就足够。流体模拟本来帧与帧之间高度连续60fps 的 mp4 体积会翻倍但视觉上的流畅度提升并不明显——除非你是要特别展示高速湍流的细节。5. GPU 加速从 NumPy 到 Taichi 的性能跨越5.1 为什么 CPU 版本跑到一定规模就卡了纯 NumPy 版本在 64x64 网格下可以勉强实时但一旦把分辨率拉到 256x256计算量会按网格面积的平方关系增长网格数多了 16 倍而每一步 Jacobi 迭代内部的循环数量也膨胀几帧之后性能就急剧下降。在我这台测试机上128x128 网格CPU 版稳定在 8 帧每秒左右而 256x256 直接掉到了 2 帧以下完全失去了“实时交互”的意义。另一个瓶颈来自 Python 的循环本身。即使用了 NumPy 的向量化advect 和 project 步骤里仍然有大量维度的索引、边界分支逻辑很难完全向量化。这个时候就该考虑“把计算挪到 GPU 上”而不是继续在 Python 层做微优化。5.2 Taichi 的引入与代码改造Taichi 是一个 Python 嵌入式 DSL它允许你用接近 Python 的语法编写高性能并行内核自动运行在 GPU 上。它的核心思想是“计算结构”、稀疏数据和自动并行化特别适合物理模拟这类网格密集型计算。安装很简单pip install taichi使用 Taichi 的关键是把数组替换成 ti.field把 Python 循环替换为 ti.kernel 装饰的函数。我改写的核心求解器结构大概是这样的import taichi as ti ti.init(archti.gpu) N 256 u ti.field(dtypeti.f32, shape(N2, N2)) v ti.field(dtypeti.f32, shape(N2, N2)) p ti.field(dtypeti.f32, shape(N2, N2)) div ti.field(dtypeti.f32, shape(N2, N2)) ti.kernel def project(): # 注意Ti 的循环自动按网格并行执行 # 红黑 Gauss-Seidel 里需要两个内循环分别只处理红色格点和黑色格点 for i, j in ti.ndrange((1, N1), (1, N1)): div[i, j] -0.5 * h * (u[i1, j] - u[i-1, j] v[i, j1] - v[i, j-1]) for _ in range(20): for i, j in ti.ndrange((1, N1), (1, N1)): p[i, j] (div[i, j] p[i-1, j] p[i1, j] p[i, j-1] p[i, j1]) * 0.25有几个细节值得注意第一Taichi 的循环是自动并行化的不需要手写 np.roll 之类的向量化技巧代码结构更接近原始数学公式可读性反而比 NumPy 版本更好。第二红黑 Gauss-Seidel 迭代需要把网格分成奇偶两类分步更新否则并行读写同一位置时会引入数据竞争。另一个高效的替代方案是用 Jacobi 迭代它天然可并行只是收敛速度略慢——我在实际测试中用 40 次 Jacobi 迭代的效果与 20 次红黑 Gauss-Seidel 接近代码也要简单很多。5.3 加速效果实测与性能对比我把 CPU 版和 GPU 版做了同样的模拟对比数据如下网格分辨率纯 NumPy CPUTaichi GPU加速倍数64x6460 fps500 fps约 8x128x1288 fps210 fps约 26x256x2562 fps60 fps约 30x512x512无法实时12 fps超过 60x这里要说明的是Taichi 在 GPU 上跑 64x64 小网格时初始化时间和内核调度本身会有固定开销性能提升反而不明显。真正发挥 GPU 优势的是大网格。我自己最终的交互版本用了 512x512 网格保持每秒 20 帧左右画面细节丰富旋涡清晰可辨鼠标交互延迟体感也在 50 毫秒以内。5.4 交互控制的实现GPU 加速之后实时交互才真正变得可行。我用 pygame 做了一个非常简单的交互界面鼠标位置映射为外部力源鼠标按住时往对应网格点注入染料和速度。核心逻辑就是在主循环里读取鼠标状态然后调用 Taichi 内核更新力源数组。这一步给整个项目带来的体验提升是质的——你不再是在看一段“预渲染的录像”而是在对着一个实时物理系统做实验。我后来在 256x256 网格下用鼠标快速拖动能看到卷起的涡旋结构像真实的烟雾一样被拉长又翻卷这个交互体验让我联想到了模拟水墨画和艺术创作的可能性。6. 常见问题与调试心得6.1 数值爆炸的三种典型表现与排查思路我调试过程中遇到最多的问题不是语法错误而是物理数值的异常。最常见的三种表现第一种是“全局白屏”。所有密度值在几步之内变成 NaN眼睛一看就是某个分母为零或者数值溢出。排查顺序是先检查 dt 是否过大过大的时间步长在半拉格朗日追踪中会把采样点推出网格有效范围一旦越界就会出现 NaN再检查速度场的初始值是否合理如果初始速度场带着巨大散度投影步骤可能在第一帧就计算出超大压力梯度。第二种是“局部亮点闪烁”。往往是某个网格点在迭代中不断累积误差最终在该点形成极大的数值尖峰。常见原因是边界条件的镜像方向写反了或者压力迭代次数不够导致局部散度没有完全消除。第三种是“整体模糊糊成一片”。这其实是数值耗散过大——半拉格朗日插值本身就是一种低通滤波会天然抹平细节。如果分辨率不够粘性系数又没设为 0模拟很快就会失去湍流结构。解决办法是提高网格分辨率或者把粘性系数调低到 0。6.2 边界伪影的排查清单边界伪影是流体模拟里最隐蔽的坑。主要表现为边界附近出现异常的流动加速、边界处密度值异常堆积、或角落出现漩涡的“伪源”。我的排查清单如下先检查 set_bnd 函数中四种场的处理是否匹配密度场为无通量边界、u 场左右反射、v 场上下反射。再检查投影步骤中散度计算的差分格式。用中心差分时索引写成 u[i1,j]-u[i-1,j]很多初学版本会用前向差分 u[i1,j]-u[i,j]会引入偏置导致边界处的系统性误差。最后检查平流步骤里反向追踪的“截断”操作。如果允许采样点超出边界才做截断那么边界处就会有大量的越界插值伪影必然出现。正确做法是在插值前先把坐标钳制到 [0.5, N0.5]。6.3 性能优化的进阶建议如果读者想在这个项目上继续深挖我最推荐的三个方向是第一多网格 V-cycle 压力求解。当网格分辨率到 1024 级别时单纯 Jacobi 迭代的收敛速度已经不够用了。多网格方法先在粗网格上快速消除低频误差再回到细网格精修整体迭代次数可以下降一个数量级。第二半拉格朗日平流的高阶插值。目前用的双线性插值会引入明显的数值耗散换成 Catmull-Rom 三次插值可以在相同网格分辨率下保留更多高频涡旋结构代价是计算量增加约一倍在 GPU 上仍然可以实时。第三游丝粒子追踪。除了密度场这种欧拉视角的可视化你还可以在速度场中放置一组无质量的标记粒子用拉格朗日视角追踪它们的运动轨迹。这类粒子流场比单纯颜色图更能表达流动的方向感两者结合的画面信息密度非常高。6.4 参数调优实战速查表参数作用推荐起始值备注dt时间步长0.05~0.2越大越不稳定越小越耗散viscosity粘性0.0~0.001模拟液体时建议设 0diffusion密度扩散0.0001~0.001越大越均匀越快压力迭代数散度消除程度20~40太少会软绵绵太多无意义网格分辨率细节程度25664 太粗糙512 需要 GPU我个人在实际调参中的体会是流体模拟的乐趣其实有一大半在“调参”上。同一套代码把粘性系数从 0 改成 0.0001画面就从干燥的“粉笔灰”变成了油润的“墨汁”这个变化物理上只有几行的差距视觉上却天差地别。所以我特别建议读者在跑通基础版本之后故意把参数调得离谱一点——比如把 dt 调大十倍亲眼看看数值爆炸的过程再回过来治它。这种“反面经验”往往是突破理解瓶颈最快的路径。这个项目后续还有很多可以扩展的玩法比如加入障碍物边界、换成彩色多流体混合、或者干脆做成手势控制的水墨画工具每一步都不需要推翻现有框架只需要在当前基础上增加一个模块就够了。