说实话我一直觉得CFD离普通程序员挺远的满屏的偏微分方程、湍流模型、并行集群总让人有一种“这不是我能碰的东西”的错觉。直到我自己用Python把二维不可压缩气流和气压场跑出来并且在屏幕上看到风场箭头随着时间一点点卷成旋涡才明白CFD的内核并不神秘无非是把连续方程离散到网格上然后反复迭代而已。这篇文章会用一套完整、可以直接运行的Python代码带你实现二维气流气压模拟并绘制风场和气压场动画。整套代码核心不到200行依赖只有NumPy和Matplotlib普通笔记本就能跑。适合想亲手实践流体数值模拟的程序员、学生和工程师。1. 写程序之前CFD到底在求解什么1.1 控制方程与控制变量二维不可压缩流场的核心控制方程是Navier-Stokes方程。听起来吓人其实拆开看就三样东西在互相拉扯速度怎么变、压力怎么产生、以及流体体积怎么守恒。速度方程描述的是“动量守恒”流体的速度会由于自身惯性继续向前跑对流项会因为内部的粘性摩擦被拉住扩散项还会因为压力差异被推着走压力梯度项。而连续性方程则是一个约束条件既然流体不可压缩那么每个网格单元里流进来的体积必须等于流出去的体积数学上表现为速度场的散度等于零。这里最麻烦的地方在于压力和速度是耦合在一起的。你更新速度的时候需要压力压力又反过来依赖速度的散度直接联立求解非常痛苦。所以在教学和工程实践中我们通常用“投影法”把这个耦合关系拆开这也是后面程序实现的核心思路。1.2 有限差分把连续问题变成网格代数问题我们不在连续空间上解方程而是把一块矩形区域切成一堆小方格每个方格中心定义一个速度值、压力值然后用相邻格点数值的差来近似偏导数。举个例子对某人来说函数f在x方向的变化率可以近似为df/dx ≈ (f[i1] - f[i-1]) / (2 * dx)这就是中心差分精度是二阶足够教学使用。这样做的好处是把偏微分方程转换成了纯粹的数组运算而Python里的NumPy恰恰是处理数组的绝佳工具。1.3 为什么选顶盖驱动腔流作为第一个实验跑CFD程序第一步千万别上来就模拟风吹过高楼大厦那是自找苦吃。我强烈建议你先复现一个经典算例——顶盖驱动腔流。这个算例的背景是一个正方形封闭空腔四壁都不动只有顶盖以恒定速度向右平移。它的伟大之处在于边界条件极其简单不需要处理复杂的进出口和障碍物但流动形态却一点不简单顶盖拖动流体向右运动到右壁后被挡下来向下卷最终在腔内形成一个巨大的主旋涡。你能清清楚楚看到风场从一层平推慢慢卷成涡旋的全过程用来验证代码和可视化逻辑再合适不过。2. 程序骨架与关键参数设计2.1 依赖与安装环境准备只需要三个东西装起来最快五分钟依赖用途安装方式NumPy核心数组运算、向量化离散计算pip install numpyMatplotlib绘制气压场云图、风场箭矢/流线、生成动画pip install matplotlibFFmpeg把动画帧编码成MP4/GIF如果只保存PNG帧可跳过apt install ffmpeg或到官网下载并加入PATH我个人建议在命令行里先跑一下ffmpeg -version确认能输出版本信息再写动画代码。很多人在这一步卡住结果程序没问题纯粹是缺了编码器。2.2 计算域、网格与边界条件的设计我采用41×41的网格计算域是边长为1.0的正方形。网格数量看起来很小但这是刻意的教学场景下你的目标是看懂代码、跑通流程而不是追求网格无关性。如果一上来就设成256×256单步迭代虽然也不慢但动画帧一多生成时间会让人瞬间失去耐心。网格步长dx和dy都是1/(41-1)0.025。这里要注意网格数包含了边界点所以相邻网格点之间的物理间距是1除以间隔数而不是1除以网格数。边界条件是整个模拟的灵魂。顶盖驱动腔流的边界条件可以用一张表说清楚边界位置水平速度u垂直速度v左壁00右壁00底壁00顶盖1.00这四个壁都是无滑移壁面粘性流体贴着壁面时速度必须与壁面一致。顶盖给了一个向右的切向速度相当于用一个移动的板带动流场旋转。压力边界采用零梯度外推也就是边界上的压力值等于紧邻内点的压力值程序中用数组切片一步就能完成。2.3 时间步长、雷诺数与稳定性物理参数方面我取流体密度rho1.0运动粘性系数nu0.01。这个nu对应的雷诺数是多少呢雷诺数的定义是特征速度×特征长度/粘性系数。这里顶盖速度U1.0腔体边长L1.0所以Re1×1/0.01100。这是腔流算例里非常经典的工况主涡清晰、数值稳定非常适合初学者观察。时间步长dt取0.001这是经过稳定性评估后选出来的。显式格式求解对流和扩散项时有CFL条件限制在一个时间步内流体移动的距离不能超过一个网格宽度。网格步长0.025速度最大大概是1所以dt至少要小于0.025。取0.001相当于每个时间步里流体最多移动0.001个单位只占网格宽度的4%非常保守但绝对稳定。工程上这叫“安定优先”你可以后面慢慢试把dt调大观察什么时候开始发散。2.4 整体流程拆解整个程序的主循环可以拆成四步用上一时刻的速度场计算“中间速度”这一步暂时不直接使用当前压力只考虑对流和扩散计算中间速度的散度利用连续性方程约束得到一个关于压力的泊松方程用Jacobi迭代求解这个压力方程得到当前时刻的压力场用压力梯度修正中间速度让速度场重新满足不可压缩约束。这四步每走完一圈就推进了一个时间步。模拟的演化过程记录成一帧帧快照最终交给Matplotlib做动画。3. 二维气流与气压求解器的核心实现3.1 中间速度的推进对流与扩散中间速度公式的离散写法是代码里最长的部分但逻辑很直白新速度 旧速度 对流项贡献 扩散项贡献。注意对流项是借助上一时刻的速度un、vn来算的这就是显式格式的特征好处是每一步的数组运算都能完全向量化坏处是必须受时间步长稳定性限制。代码实现起来就是一堆数组切片的组合。用切片代替循环能充分利用NumPy的底层C语言速度。很多新手第一次接触可能不习惯这种密密麻麻的索引表达式我的经验是先在纸上画一个小网格标出i和j的方向再用难懂的部分逐项对照很快就能明白。3.2 压力泊松方程与Jacobi迭代速度中间场算完之后它的散度一般不等于零。连续性方程要求速度场无散度所以我们需要建立一个压力场让压力梯度去“顶”掉这个多余的散度。数学推导可以把速度修正公式两边取散度联立连续性约束最终得到一个二阶椭圆方程这就是压力泊松方程。代码里用Jacobi迭代来解原因是实现最简单、最透明。每次迭代时把每个网格点的新压力用周围四个点的旧压力加权平均再减掉一个由速度散度构造的源项。迭代若干次后压力分布就会收敛到满足泊松方程的解。对于41×41的小网格迭代50次已经能得到相当稳定的气压场。3.3 速度修正与边界处理有了压力后就可以用压力梯度修正速度了。修正公式也很简单新的速度为中间速度减去压力梯度乘以dt/ρ。修正完之后再重新施加边界条件。边界处理是整个程序最容易出错的地方。如果忘记把顶盖速度设回1.0顶盖的剪切驱动力就会消失整个腔流会慢慢停下来涡旋永远不会出现。我踩过这个坑当时盯着动画看了半天画面里一片死水最后检查了整整二十分钟才发现原因。3.4 可运行的完整求解器代码下面给出我整理好的完整核心代码。参数、边界、迭代和快照记录都打包在了一个脚本里建议你新建一个cavity.py把代码完整粘贴进去运行。import numpy as np import matplotlib.pyplot as plt # 网格与物理参数 nx, ny 41, 41 dx dy 1.0 / (nx - 1) nt 6000 nit 50 dt 0.001 rho 1.0 nu 0.01 save_step 50 # 初始化速度场、压力场 u np.zeros((ny, nx)) v np.zeros((ny, nx)) p np.zeros((ny, nx)) # 网格坐标用于绘图 X, Y np.meshgrid(np.linspace(0, 1, nx), np.linspace(0, 1, ny)) def set_velocity_boundaries(u, v): u[:, 0] 0 # 左壁 u[:, -1] 0 # 右壁 u[0, :] 0 # 底壁 u[-1, :] 1 # 顶盖向右移动 v[:, 0] 0 v[:, -1] 0 v[0, :] 0 v[-1, :] 0 def set_pressure_boundaries(p): p[:, 0] p[:, 1] p[:, -1] p[:, -2] p[0, :] p[1, :] p[-1, :] p[-2, :] set_velocity_boundaries(u, v) snapshots [] for n in range(nt): un u.copy() vn v.copy() # 中间速度对流 扩散 u[1:-1, 1:-1] ( un[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2]) - vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1]) nu * dt / dx**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] un[1:-1, 0:-2]) nu * dt / dy**2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] un[0:-2, 1:-1]) ) v[1:-1, 1:-1] ( vn[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) - vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) nu * dt / dx**2 * (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] vn[1:-1, 0:-2]) nu * dt / dy**2 * (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] vn[0:-2, 1:-1]) ) # 压力泊松方程Jacobi迭代 for it in range(nit): pn p.copy() div_u ( (u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy) ) p[1:-1, 1:-1] ( (pn[1:-1, 2:] pn[1:-1, 0:-2]) * dy**2 (pn[2:, 1:-1] pn[0:-2, 1:-1]) * dx**2 ) / (2 * (dx**2 dy**2)) - ( dx**2 * dy**2 / (2 * (dx**2 dy**2)) * rho / dt * div_u ) set_pressure_boundaries(p) # 用压力梯度修正速度 u[1:-1, 1:-1] - dt / rho * (p[1:-1, 2:] - p[1:-1, 0:-2]) / (2 * dx) v[1:-1, 1:-1] - dt / rho * (p[2:, 1:-1] - p[0:-2, 1:-1]) / (2 * dy) set_velocity_boundaries(u, v) if n % save_step 0: snapshots.append((n * dt, p.copy(), u.copy(), v.copy())) if n % 1000 0: print(fstep {n}/{nt}, mean |u| {np.abs(u).mean():.4f}) print(模拟完成共保存, len(snapshots), 帧数据)这段代码跑完后snapshots里就存了从初始时刻到结束时刻共120帧的流场数据。这里的p.copy()、u.copy()必须写因为Python的等号赋值只是引用传递不复制的话后面求解器修改数组会把你之前存的快照全部覆盖掉动画就变成只有最后一帧重复播放了。4. 把风场和气压场画成动画4.1 静态可视化气压云图与风场流线在生成动画之前先画一帧静态图检查结果。气压场用contourf画云图非常适合因为压力是连续标量场颜色等高图能清晰展示高压区和低压区。风场则有两个选择quiver画箭头图适合看速度方向streamplot画流线图适合看流动的整体走向。我一般把气压云图和流线叠加在同一张图上气压用背景颜色表示流线用黑色线条叠加这样一张图就能同时看到两个物理量的空间关系。如果只想画风场也可以用quiver每隔几个点采样一次避免箭头太密集。4.2 快照数据与动画生成的分离设计很多人写CFD动画时喜欢在求解循环里一边算一边画这其实不是好习惯。求解步很密集每步都重绘一次性能极差而且如果你想把动画帧率调成10fps又希望模拟跑6000步直接画图就得画6000帧完全没必要。我采用的办法是求解阶段只负责把每50步的压力场和速度场复制到snapshots列表里模拟全部跑完后再统一做可视化。这样模拟和绘图两个阶段完全分离你可以反复调节绘图参数而不必重新计算流场省时省力。4.3 用FuncAnimation生成MP4动画Matplotlib的FuncAnimation是专门用来做这种逐帧动画的。它的思路很简单你提供一个“画第i帧”的函数动画模块把这个函数反复调用并把结果按顺序输出成视频。记得在动画函数里调用ax.clear()清除上一帧内容否则旧图形会和当前帧重叠在一起。参数interval控制的是毫秒间隔fps则由保存时指定。dpi数值越大视频越清晰但文件体积也会越大。from matplotlib.animation import FuncAnimation fig, ax plt.subplots(1, 1, figsize(6.4, 6)) fig.tight_layout() def draw_frame(index): t, p_i, u_i, v_i snapshots[index] ax.clear() cf ax.contourf( X, Y, p_i, levels60, cmapRdBu_r, vmin-np.max(np.abs(p_i)), vmaxnp.max(np.abs(p_i)) ) ax.streamplot( X, Y, u_i, v_i, color0.25, density1.6, arrowstyle-, linewidth0.9 ) ax.set_title(ft {t:.3f} s, wind pressure) ax.set_xlabel(x) ax.set_ylabel(y) ax.set_aspect(equal) return cf anim FuncAnimation(fig, draw_frame, frameslen(snapshots), interval120) anim.save(wind_pressure_field.mp4, writerffmpeg, fps10, dpi120) print(动画已保存)这个代码里我用了vmin和vmax随当前帧的最大压力绝对值动态变化好处是颜色对比度始终饱满。但如果你希望不同帧之间的颜色系统完全可比就必须在循环前先确定全过程的压力最大值然后固定vmin、vmax。4.4 备选方案逐帧存PNG再合成如果你的机器上FuncAnimation保存时遇到编码器兼容问题还有一个土办法极其好用把每一帧用savefig保存成PNG图片然后用FFmpeg命令批量合成。这个方法没有任何神秘性还方便你做后期筛选。import os os.makedirs(frames, exist_okTrue) for i, (t, p_i, u_i, v_i) in enumerate(snapshots): fig2, ax2 plt.subplots(figsize(6.4, 6)) ax2.contourf(X, Y, p_i, levels60, cmapRdBu_r) ax2.streamplot(X, Y, u_i, v_i, color0.25, density1.6) ax2.set_title(ft {t:.3f} s) fig2.savefig(fframes/frame_{i:04d}.png, dpi120) plt.close(fig2)然后在命令行里执行ffmpeg -r 10 -i frames/frame_%04d.png -vcodec libx264 -pix_fmt yuv420p wind_pressure_field.mp4-r 10表示每秒播放10帧-i frames/frame_%04d.png表示按四位数字编号连续读取图片。这个方案的优点是每一步你都能看到进度而且某几帧画坏了可以直接单独改。4.5 让画面更直观的调节技巧画风场时很多新手抱怨streamplot出来的线条太稀疏或者太密其实density参数就是干这个的。数值越大线条越密集1.2到2.0是常见取值范围。流线的颜色统一设成深灰色0.25不会抢气压云的视觉焦点。气压场的色标选择上RdBu_r红蓝反向是CFD里的经典选择红色代表正压、蓝色代表负压。如果你的模拟场景里有明确的压力取值范围最好手动设置vmin、vmax不然Matplotlib会自动根据每一帧数据缩放颜色区间导致不同帧之间颜色含义不一致看动画时容易误判。5. 踩过的坑与排查思路5.1 常见问题速查表现象可能原因处理办法几步后速度出现NaN或极大值dt过大不满足CFL稳定性条件把dt降到0.0005或更小重试压力场呈棋盘状交替高低网格点定义与压力梯度离散不匹配改用交错网格或增加压力迭代次数观察是否缓解动画里流场一片静止顶盖边界速度被覆盖为零检查每次主循环末是否重新设置u[-1,:]1保存MP4报FFmpeg相关错误未安装ffmpeg或未加入PATH命令行运行ffmpeg -version验证动画播放速度过快或过慢fps与模拟时间尺度不匹配调节fps、interval参数压强颜色区间不断跳动vmin/vmax未固定计算全局压力范围后手动设定求解器太慢网格过大或nit过多减少网格到31×31或把nit降到305.2 压力棋盘振荡的根本原因这是所有CFD初学者的噩梦。当速度分量和压力定义在同一个网格点上非交错网格时离散后的压力方程和速度修正中会出现一种“失耦合”现象相邻网格点的压力呈棋盘状交替增大和减小而散度条件却依然被满足于是这种虚假振荡在迭代中不会被识别为错误。为什么顶盖腔流的教学代码里这个问题不明显因为均匀网格、简单边界条件和Jacobi迭代的阻尼特性给了它一定容错空间。但如果你把网格加密到128以上或者改变边界条件棋盘模式很可能立刻冒出来。专业CFD里最干净的办法是把u、v、p分别定义在交错的网格位置也就是所谓的MAC交错网格。我建议初学者先把当前版本跑通理解流程后再去研究交错网格一次学太多反而容易乱。5.3 一旦发散优先检查这五个位置模拟过程中出现数值爆炸时我有一套固定的排查顺序效率非常高检查dt是否太大。这是发散的第一嫌疑犯特别是把网格加密后CFL条件更严格原来的dt可能立刻失效。检查边界条件是否完整。有没有某个壁面的速度漏设置导致边界附近出现巨大梯度。检查压力迭代次数nit是否够多。压力没解准速度修正就会过冲表现为局部振荡逐步放大。检查中间速度计算里的切片索引有没有写反。比如某个数组维度顺序错误会让扩散贡献变成负向反而加热流场。检查每次迭代是否把上一时刻的速度复制给了un和vn。如果少了这步求解器用的是当前步已经更新过的速度会引入隐式耦合破坏稳定性。5.4 动画性能优化心得动画慢的根源通常是streamplot每帧都要从零开始计算插值路径。我实测下来41×41的网格已经很快但如果你加密到81×81生成一段10秒的MP4可能要等上几分钟这很正常。想提速有几个办法一是降低save_step对应的保存频率减少帧数二是把streamplot替换成quiverarrow图的计算量远小于流线图三是减小density参数。如果只是初略观察流场形态完全可以用快照里第几帧抽样来画成多图排布而不必非要做视频。最后再多说一句这套代码虽然小但已经覆盖了CDF求解器最基本的几个支柱——网格离散、显式推进、压力泊松迭代、无散度修正、边界处理。把它跑通之后你完全可以拿着这个骨架去改雷诺数、换边界条件、甚至加一个圆柱障碍物进去。比如把顶盖速度提高一个量级就能看到涡旋慢慢失稳、流动进入非定常状态的精彩过程。当你亲眼看着屏幕上原本整齐的流线开始扭动那种成就感跟跑通任何Web项目完全不一样。