企业官网建设流程全解析

说实话，最开始我抱着"顺便玩玩"的心态去写这个二维气流模拟，真正开始之后才发现，把一个二维不可压缩流场从零算出来，比想象中有意思得多：整个过程不依赖任何商业软件，只用 Python 和 NumPy 就能在笔记本上跑出一个完整的风场和气压场动画。这种感觉就像你终于知道"流场图是怎么画出来的"，而不是每次都在等别人算好的结果。

这篇文章我尽量写得像一次手把手的实操记录，适合那种会一点 Python 基本语法、但对计算流体力学还比较陌生的读者。我也会把那些代码注释里往往不写的坑讲清楚：为什么有些参数不能乱调，为什么边界上的涡量要额外处理，为什么别人跑的动画不抖而你的却发散。如果看完你也能复现出一个二维顶盖驱动流的模拟动画，那这篇教程的目的就达到了。

1. 为什么选"涡量–流函数"而不是压力–速度直接求解

1.1 两种方程形式的核心差异

在做二维不可压缩流动的数值模拟时，第一反应往往是把原始的纳维–斯托克斯方程直接拿来离散。二维不可压缩 N-S 方程写成无量纲形式是这样的：

[ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 ]

[ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} =-\frac{\partial p}{\partial x}+\frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right) ]

[ \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} =-\frac{\partial p}{\partial y}+\frac{1}{Re}\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) ]

这套方程看起来很完整，但实际写代码时会立刻碰到一个麻烦：连续性方程里没有压力项，可是动量方程里压力是未知量，它并不存在一个独立的"压力演化方程"。如果直接做时间和空间离散，压力和速度会像捉迷藏一样互相依赖，需要额外引入投影法或压力 Poisson 方程求解步骤。对初学者来说，这一整套正确写出来很容易卡住。

所以，在二维问题上我更推荐避开原始变量，改用涡量–流函数形式。定义涡量：

[ \omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} ]

因为流动不可压缩，二维速度场一定可以表示成一个流函数 ψ 的导数：

[ u=\frac{\partial \psi}{\partial y}, \quad v=-\frac{\partial \psi}{\partial x} ]

把涡量的定义代进去，就得到了著名的流函数 Poisson 方程：

[ \frac{\partial^2 \psi}{\partial x^2}+\frac{\partial^2 \psi}{\partial y^2}=-\omega ]

再把涡量定义代入 N-S 方程的旋度形式，得到涡量输运方程：

[ \frac{\partial \omega}{\partial t}+u\frac{\partial \omega}{\partial x}+v\frac{\partial \omega}{\partial y} =\frac{1}{Re}\left(\frac{\partial^2 \omega}{\partial x^2}+\frac{\partial^2 \omega}{\partial y^2}\right) ]

这两条方程一个管"涡量怎么随时间演化"，一个管"根据当前涡量怎么求出流函数，再反推出速度"。压力项在这套表述里被彻底消掉了，你不必再去求解繁琐的压力 Poisson 方程，主程序的结构一下子清晰很多。

1.2 无量纲化为什么只留下一个雷诺数

标题里要同时展示风场和气压场，而对可压缩气流来说，气压场本来应该是热力学状态量。为了让问题既接近标题又方便入门，我建议先处理二维不可压缩流动，这里的"气压场"实际可以对应为压力场。在一套无量纲方程里，真正决定流场形态的只有一个参数：雷诺数 Re。

无量纲化的思路并不复杂。取特征长度 L 为方腔边长，特征速度 U 为顶盖驱动速度，特征时间就是 L/U，然后用这些量把方程里的变量全部换成无量纲量。做完之后，原来方程里的运动粘度 ν 会合并成：

[ Re=\frac{UL}{\nu} ]

这个 Re 一旦定了，不管你是用厘米还是米做单位，流场的涡结构形态是同一套。这是所有 CFD 问题里特别值得先意识到的一点：先无量纲化，后面调参数会方便很多。

关于压力的显示，我后面会额外用泊松方程从速度场重建压力场，画成云图当作"气压场"展示。这其实是很多教学程序里会忽略的一步，但标题既然提到气压场，我觉得值得把它补上，而且实现起来也不复杂。

1.3 两条方程的数值角色分配

把两条方程放在一起看，计算节奏很明确：

• 涡量输运方程是时间推进方程，知道当前时刻的 ω、u、v，就能算出下一时刻的 ω。
• 流函数 Poisson 方程是椭圆型方程，每个时间步都要重新解一遍，得到当前涡量对应的 ψ。
• 有了 ψ，用导数关系马上得到新的 u 和 v，然后再进入下一个时间步。

也就是说，整个程序的骨架不是一条方程走到底，而是"涡量推进一次 → 解一次泊松方程 → 更新速度 → 再推进涡量"这样循环。后面第 2 章我会把这套循环的代码骨架完整写出来。

2. 程序整体骨架：从数据结构到主循环

2.1 数组怎么组织才不容易乱

第一版程序不需要太复杂的类结构，用函数加全局数组完全够用。关键是数组的维度顺序要一致。我是这样安排的：

• 计算域是 1×1 的方腔，网格数 nx=101、ny=101（含边界点）；
• 网格间距均匀，dx=dy；
• 所有物理场都用二维 NumPy 数组保存，第一维是 y 方向索引，第二维是 x 方向索引。

为什么要坚持"先 y 后 x"？因为 matplotlib 画图时，数组的第 0 维对应纵坐标方向，第 1 维对应横坐标方向。你如果一开始按 (nx, ny) 组织数组，后面画图时要不停转置，非常容易把自己绕晕。不如第一次就统一成 (ny, nx)。

初始化代码大概长这样：

import numpy as np nx, ny = 101, 101 # 网格点数 Lx, Ly = 1.0, 1.0 # 计算域尺寸 dx = Lx / (nx - 1) dy = Ly / (ny - 1) Re = 100.0 # 雷诺数 u_top = 1.0 # 顶盖驱动速度 omega = np.zeros((ny, nx)) # 涡量场 psi = np.zeros((ny, nx)) # 流函数场 u = np.zeros((ny, nx)) # x 方向速度 v = np.zeros((ny, nx)) # y 方向速度

顺便说一句，u 表示 x 方向速度，这个符号约定全篇统一，后面代码里看着也不会迷糊。

2.2 时间推进主循环的顺序

整个模拟最核心的地方，就是这个主循环。每一步做的事情，顺序一定不能乱：

1. 根据当前流函数 ψ 计算速度场 u、v；
2. 根据当前速度和涡量，离散涡量输运方程，得到下一时刻的涡量；
3. 用新涡量作为右边源项，解流函数 Poisson 方程，得到新的 ψ；
4. 更新边界上的涡量；
5. 记录当前步数、残差、时间，然后进入下一步。

写成伪代码是这样：

for step in range(max_steps): u = dpsi_dy(psi) v = -dpsi_dx(psi) omega_new = omega + dt * compute_omega_rhs(omega, u, v, Re, dx, dy) omega = update_boundary_omega(omega_new, psi, u_top) psi = solve_poisson_sor(omega, psi) if step % save_interval == 0: save_snapshot(...)

注意第 4 步边界涡量的更新，务必放在"解完泊松方程之后、下一次计算涡量输运之前"。因为边界涡量依赖最新的流函数，如果在错误的时间更新，边界条件会滞后一个时间步，高雷诺数下容易出现边界附近的数值振荡。

2.3 为什么先选顶盖驱动流作为验证算例

流体的数值模拟不是从零随便造一个场就开始跑，那样很难发现程序哪里错了。标准做法是先找一个有大量公开结果的经典算例来验证程序。我选择的是顶盖驱动流（lid-driven cavity flow）：一个正方形方腔，底部和左右壁面静止，顶盖以恒定速度向右滑动。

这个算例优点非常明显：

• 几何超级简单，就是一个正方形，不用处理进出口边界；
• 不需要设计进口速度剖面、出口边界，省掉一堆麻烦；
• 雷诺数从 100 到 1000，流场会出现清晰的主涡和角涡，视觉效果很适合做动画；
• 文献里的结果非常多，拿自己的流场和标准结果对比，很快能判断程序对不对。

我在读研那会儿第一次写 CFD 程序，导师给的建议就是"先别惦记圆柱绕流，先把方腔算清楚"。这个忠告很值钱，因为方腔流把边界条件、迭代求解、时间步长控制这些基础能力全部练了一遍，后面再去算复杂几何，其实只是在这个框架上加网格和边界处理而已。

3. 核心求解器：涡量输运方程与泊松方程的实现

3.1 对流项差分格式的取舍

涡量输运方程里最需要小心的是对流项 (u\frac{\partial \omega}{\partial x}+v\frac{\partial \omega}{\partial y})。如果用最简单的中心差分：

[ \left.\frac{\partial \omega}{\partial x}\right|{i,j} \approx \frac{\omega{i,j+1}-\omega_{i,j-1}}{2\Delta x} ]

精度确实是二阶，但显式时间推进下，这种格式在网格雷诺数较大时会产生强烈的非物理振荡，气流场会像被揉皱的纸一样到处是锯齿。对于第一版程序，我更推荐迎风差分：根据速度方向选择上游一侧的差分。

以 x 方向为例：

[ u\frac{\partial \omega}{\partial x} \approx \begin{cases} u_{i,j}\frac{\omega_{i,j}-\omega_{i,j-1}}{\Delta x}, & u_{i,j}>0 \[4pt] u_{i,j}\frac{\omega_{i,j+1}-\omega_{i,j}}{\Delta x}, & u_{i,j}<0 \end{cases} ]

迎风格式牺牲了一点精度，但换来了很好的稳定性，对教学用的一阶程序来说，这是性价比极高的选择。等你跑通整个流程之后，可以再把对流项换成二阶中心差分或者 Adams-Bashforth，那时候再跟一阶结果对比，体会会更深刻。

3.2 涡量方程的离散与推进代码

把扩散项和二阶导数离散放在一起，整个涡量导数函数写成这样：

def compute_domega_dt(omega, u, v, Re, dx, dy): domega = np.zeros_like(omega) # 内部格点（避开边界） for j in range(1, ny - 1): for i in range(1, nx - 1): # 对流项迎风差分 if u[j, i] > 0: u_dw = u[j, i] * (omega[j, i] - omega[j, i - 1]) / dx else: u_dw = u[j, i] * (omega[j, i + 1] - omega[j, i]) / dx if v[j, i] > 0: v_dw = v[j, i] * (omega[j, i] - omega[j - 1, i]) / dy else: v_dw = v[j, i] * (omega[j, i + 1] - omega[j, i]) / dy # 二阶导（扩散项） lap = (omega[j, i + 1] - 2.0 * omega[j, i] + omega[j, i - 1]) / dx**2 lap += (omega[j + 1, i] - 2.0 * omega[j, i] + omega[j - 1, i]) / dy**2 domega[j, i] = -u_dw - v_dw + (1.0 / Re) * lap return domega

然后时间推进用最简单的显式欧拉：

dt = 0.001 # 后面会讲怎么根据 CFL 条件取 domega = compute_domega_dt(omega, u, v, Re, dx, dy) omega_new = omega.copy() omega_new[1:-1, 1:-1] += dt * domega[1:-1, 1:-1]

注意一定用omega.copy()生成新的数组。我在初期吃过亏：直接写omega += ...，结果边界更新时使用了已经被污染的内部值，整个流场在几百步后变得乱七八糟。

3.3 流函数泊松方程与 SOR 迭代

每个时间步结束后，都要解一次：

[ \frac{\partial^2\psi}{\partial x^2}+\frac{\partial^2\psi}{\partial y^2}=-\omega ]

这是一条椭圆型方程，不会随时间演化。求解方式有很多，第一版程序我用 SOR 迭代，也就是逐次超松弛迭代。基本原理先回忆一下：普通的高斯–赛德尔迭代，走完一圈更新之后，把新值和旧值再做一个加权混合，权重就是松弛因子 ρ。写成公式：

[ \psi_{i,j}^{(k+1)}= (1-\rho)\psi_{i,j}^{(k)} +\frac{\rho}{4}\left[ \psi_{i+1,j}^{(k+1?)}+\psi_{i-1,j}^{(k)}+\psi_{i,j+1}^{(k?)}+\psi_{i,j-1}^{(k)} +\omega_{i,j}\Delta x^2 \right] ]

实际编写时，SOR 天然是串行更新的，用 Python 三重循环最直观。这里先给出内部点的基本写法：

def poisson_sor(omega, psi, rho=1.5, max_iter=500, tol=1e-6): res = 1.0 it = 0 while res > tol and it < max_iter: res = 0.0 for j in range(1, ny - 1): for i in range(1, nx - 1): psi_old = psi[j, i] psi_new = 0.25 * ( psi[j + 1, i] + psi[j - 1, i] + psi[j, i + 1] + psi[j, i - 1] + omega[j, i] * dx * dy ) psi[j, i] = (1 - rho) * psi_old + rho * psi_new res = max(res, abs(psi[j, i] - psi_old)) it += 1 return psi, res, it

注意我这里系数写成0.25的前提是 dx=dy。如果你后面用了不均匀网格或长方形计算域，就必须换成更一般的五点差分系数，否则会得到错误结果。这里正是新手最容易忽略的地方。

松弛因子 ρ 的选择也值得说。ρ 太小，收敛慢；ρ 接近 2，会过冲甚至发散。在方腔网格 101×101 的情况下，ρ=1.5 是一个非常稳妥的选择，一般几百次迭代就能把残差压到 1e-6 以下。

3.4 时间步长和 CFL 条件

显式时间推进必然有时间步长上限，否则程序会发散。CFL 条件通常在两个地方同时起作用：

• 对流稳定性：(dt \le cfl \cdot \Delta x / U_{max})
• 扩散稳定性：(dt \le cfl \cdot \Delta x^2 Re / 2)

实际写代码时，我建议每次都动态计算：

Umax = max(np.max(np.abs(u)), np.max(np.abs(v)), 1e-10) dt_cfl = 0.5 * dx / Umax dt_diff = 0.5 * dx * dx * Re / 2.0 dt = min(dt_cfl, dt_diff)

这个方法既安全又可以适应不同网格。比如 101×101 网格下，dx=0.01，若 Re=100，扩散限制大约是 0.0025；但如果 Re 提高到 1000，极限时间步会显著变小，必须动态调整，不能写死。

初版程序我建议用dt = 1e-4~1e-3这种保守取值先跑通流程，然后再用上面的公式自动选择。能看到动画稳定推进，比省那几步计算时间重要得多。

4. 边界条件的处理：最容易出错的位置

4.1 固壁上的涡量如何由流函数反推

涡量输运方程里并没有边界涡量的演化公式，壁面上的涡量必须靠流函数反推。这个环节网上资料写得最乱，我直接给出我验证过的写法。

对于静止壁面，用流函数的二阶导数近似：

[ \omega_{wall}=-\frac{2(\psi_{inner}-\psi_{wall})}{(\Delta n)^2} ]

其中 Δn 是到壁面的垂直网格间距。比如底边：

omega[0, :] = -2.0 * (psi[1, :] - psi[0, :]) / dy**2

左右边类似。对顶盖来说，因为顶盖本身有速度 (u_{top})，壁面切向速度不再是 0，所以公式后面要多一项驱动项：

[ \omega_{top}=-\frac{2(\psi_{inner}-\psi_{top})}{dy^2} -\frac{2u_{top}}{dy} ]

对应的代码就是：

omega[-1, :] = -2.0 * (psi[-2, :] - psi[-1, :]) / dy**2 - 2.0 * u_top / dy

这几个公式我建议直接用，不要试图自己改符号。初版程序只要这一点写错，流场要么完全退化，要么在壁面附近出现严重锯齿。如果你想要更高精度，可以考虑加入单侧三阶近似，但对第一版来说，上面公式已经够用。

4.2 初场怎么给，要不要加扰动

顶盖驱动流的好处是初场可以全设为 0，也就是说一开始腔体内流体完全静止。驱动力完全来自顶盖壁面的剪切，这一项会不断把涡量从边界带进内部流场，逐渐形成主涡。对教学程序来说，这比"人为加一个扰动"更干净，因为你可以确信流动结构是方程自己长出来的，不是初始条件预设的。

唯一需要注意的是，因为从零开始演化需要时间，前面几百步动画看起来会很无聊：流场几乎不动。解决办法有两个：一是先跑几百步再开始保存动画；二是把保存帧间隔设大一点，只保存每 20 步或 50 步的结果。

4.3 常见翻车现场排查表

我把实际调试中遇到的典型问题整理了一下，遇到动画变成雪花纸片的时候，对着表排查会快很多。

这个表是我个人的调试清单，思路就是"先看残差、再看边界、最后才怀疑格式"。很多人一上来就怀疑差分格式精度不够，其实绝大多数初版程序的崩溃，都是时间步和边界条件的问题。

5. 可视化：从静态场到风场气压场动画

5.1 先用静态图确认流场长对了

动态之前，我建议先画两张静态图：一张流函数等值线（相当于流线），一张速度矢量图（风场）。这样能快速判断程序是否基本正确。

import matplotlib.pyplot as plt X, Y = np.meshgrid( np.linspace(0, Lx, nx), np.linspace(0, Ly, ny) ) fig, axes = plt.subplots(1, 2, figsize=(12, 5)) # 流函数等值线 ax = axes[0] cs = ax.contourf(X, Y, psi, levels=50, cmap='RdBu_r') fig.colorbar(cs, ax=ax) ax.set_title('Stream Function') # 速度矢量 ax = axes[1] ax.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2]) ax.set_title('Velocity Field') plt.tight_layout() plt.show()

这里有一个小坑：X是 x 方向坐标，Y是 y 方向坐标，但我们的数组是psi[j, i]，所以如果直接contourf(X, Y, psi)可能会发现坐标反了。正确做法是保证psi的索引第一维对应Y的第二维对应X，这正是 2.1 节强调(ny, nx)组织数组的原因。

5.2 用 FuncAnimation 制作演化动画

静态图确认无误后，就可以把求解一步 → 更新绘图对象这个过程做成动画。用 matplotlib 的FuncAnimation最方便，基本结构是：

from matplotlib.animation import FuncAnimation, PillowWriter fig, axes = plt.subplots(1, 2, figsize=(12, 5)) ax1 = axes[0] ax2 = axes[1] cf = ax1.contourf(X, Y, psi, levels=50, cmap='RdBu_r') qv = ax2.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2]) def update(frame): global omega, psi, u, v for _ in range(steps_per_frame): # 这里放主循环里的推进逻辑 u = dpsi_dy(psi, dx) v = -dpsi_dx(psi, dy) domega = compute_domega_dt(omega, u, v, Re, dx, dy) omega_new = omega.copy() omega_new[1:-1, 1:-1] += dt * domega[1:-1, 1:-1] omega = update_boundary_omega(omega_new, psi, u_top) psi = poisson_sor(omega, psi)[0] # 更新绘图对象 ax1.clear() ax1.contourf(X, Y, psi, levels=50, cmap='RdBu_r') ax2.clear() ax2.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2]) return ax1, ax2 anim = FuncAnimation(fig, update, frames=200, interval=30) anim.save('cavity.gif', writer=PillowWriter(fps=20))

这里每帧之间多跑几步steps_per_frame，这样即使时间步很小，动画也能在有限帧数里看到明显的流动演化。保存 gif 用PillowWriter，保存 mp4 则需要系统里装了 ffmpeg，writer 换成FfmpegWriter即可。

5.3 气压场：用泊松方程从速度场重建

标题里提到气压场动画，有些读者会觉得难以直接呈现。对于不可压缩流而言，压力并不是从状态方程解出来的，而是满足一个压力泊松方程。二维不可压缩流的压力方程形式为：

[ \nabla^2 p = -2\left(\frac{\partial u}{\partial x}\frac{\partial v}{\partial y} -\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}\right) ]

可以把它当作又一个泊松方程来解，源项是速度梯度组合。压力和流函数一样用 SOR 迭代求解，只是右侧源项换成上面的表达式。边界条件相对复杂，初版可以简单采用 Neumann 边界，也就是法向压力梯度由动量方程给出。这一步做出来后，就能把压力场画成云图，和风场放一起展示。

如果你觉得边界处理太繁琐，也有更省事的办法：先把主涡跑出来，然后用流函数等值线做定性参照，把压力场作为辅助图展示。但对想认真理解气压场动画的读者，我建议还是把压力泊松方程完整加上去，它会让整套模拟看起来更接近真实 CFD 后处理。

5.4 动画文件太大的处理技巧

跑过几次之后你会发现，数值模拟的动画文件非常大，尤其是一帧里画了 200 个箭头加上 50 层等值线。这里有几个实际可操作的小技巧：

• 矢量化箭头每隔 2 个或 3 个点画一个，画多了肉眼看不出更密，只会让文件膨胀；
• 等值线层数从 50 降到 20~30，动画流畅度提升明显；
• gif 颜色用cmap的细分数量少一点，也能显著缩小体积；
• 保存帧间隔控制在 100~300 帧，不要每步都存，否则文件动辄几十兆。

6. 调试验证、调参与一点性能心得

6.1 怎么判断算出来的流场是对的

程序跑出动画后，第一件事不是发朋友圈，而是对比参考数据。顶盖驱动流的标准结果非常多，以 Re=100 为例，文献中主涡中心大约在 (x\approx0.6, y\approx0.76) 附近；Re=400 时主涡中心会移到腔体上方偏右的位置；Re=1000 时还会出现明显的左下角和右下角角涡。你可以用自己的流函数最小值点位置和这些值对比，误差在网格间距范围内就说明程序没问题。

如果你画的流线图里只有一个巨大的中心涡，没有右上方角涡，那大概率是计算时间不够，流场还没完全发展，多跑几百步再看。

6.2 网格分辨率、雷诺数与时间步长的经验组合

不带性能优化的纯 Python 三重循环版本，网格数别一上来就选 257×257，否则每次泊松迭代的耗时就会让你怀疑人生。我建议按照下表来推进：

上表是我在普通笔记本电脑上的实测感受，不是绝对基准，但能帮你建立直觉：越细的网格，时间步必须压得越小，总计算量增长非常快。这也是为什么很多人第一次玩 CFD 都觉得"为什么别人动画那么快，我的却像放幻灯片"，本质是网格和 dt 没匹配好。

6.3 用 NumPy 向量化替代部分循环

第一版代码里，为了可读性我故意用了三层循环。但如果你想让网格数跑到 201×201 以上，建议把涡量导数的内部循环改成 NumPy 切片运算。比如扩散项整个计算域可以这样写：

lap = ( omega[1:-1, 2:] + omega[1:-1, :-2] + omega[2:, 1:-1] + omega[:-2, 1:-1] - 4.0 * omega[1:-1, 1:-1] ) / dx**2

对流项的迎风差分也可以根据速度场符号，用np.where一次性选择左右差分。这样改写之后，计算速度通常会快一个数量级。SOR 迭代因为天然串行，比较难整体向量化，但可以逐行更新或者用更少迭代次数的多重网格法，初学者暂时不用太纠结。

动画展示环节对性能更敏感，建议绘图数据只取降采样的网格，比如u[::2, ::2]。这样箭头数量变成原来的四分之一，绘图耗时大幅下降，而流场趋势依然清晰。

6.4 我调试几天下来最想提醒的三件事

最后分享三条实际调试中最有价值的经验。

第一，一定要保留每一步的残差和最大涡量作为"生命体征"。我每次迭代都会打印：

step=1200 time=1.200 res_psi=3.2e-7 max_omega=4.20

一旦最大涡量开始持续上涨，说明大概率快发散了，这时候暂停程序去查参数，比等它变成 NaN 再回头找原因高效得多。

第二，初次跑通后，不要着急调高雷诺数，先把 101×101、Re=100 的结果彻底看明白。真正理解主涡、角涡、边界层这些结构之后，再去挑战 Re=1000 会更有底。很多教程把 Re=1000 当作默认展示，但初学者直接复现时，边界条件、迭代格式和 dt 三处只要有一点点问题，流场就会异常花哨，容易误以为是"物理现象"。

第三，动画和计算解耦。我强烈建议先算出结果、保存数据，再单独做可视化。如果你想边算边画，update函数里每帧做几十步计算还可以接受，但如果是几十万步的计算，实时绘图会成为瓶颈。先存数据，后处理动画，这套工作流在真实 CFD 项目里也更常见。

如果你按照这篇文章把顶盖驱动流跑出来了，下一步可以试试两个方向：一是把顶盖驱动改成圆柱绕流，引入更复杂的几何边界和进出口条件；二是把显式欧拉时间推进换成二阶 Runge-Kutta，你会发现同样的网格下稳定性和精度都能进一步提升。祝玩得开心。