简介：这是面向流体力学初学者与MATLAB用户的二维非定常Navier-Stokes方程求解示例，聚焦不可压缩流体在二维空间内的运动模拟，适合学习有限差分、时间推进及压力修正等核心数值方法。资源共76个文件，含44个m脚本、29个png结果图、1个txt说明及1个html演示页，覆盖网格定义、边界条件、初始化、数值求解与可视化全流程，压缩包约463KB。文件组织清晰，主程序与各功能模块分离，便于按步骤理解代码逻辑。目前已有1941人浏览学习。通过该示例可掌握使用quiver绘制速度矢量场、contourf显示流线等可视化技巧，同时加深对N-S方程中压力梯度项、粘性扩散项及连续性约束的实际认识，为后续研究飞行器空气动力学、水力学等工程问题打下基础。

1. 二维非定常 N-S 方程：这个 MATLAB 示例能让你把“流动仿真”从玄学变成看得见的曲线

做过流体仿真的同学都有这种体验：学了一堆 N-S 方程推导，真要自己在 MATLAB 里跑一个非定常流场，却不知道第一行代码该写什么。网上的开源求解器动不动就是 C++ 几万行，装依赖就得折腾半天；教科书附录里的示例代码又往往只给定常版，一碰到时间项就集体沉默。这份“二维 Navier-Stokes（非定常）MATLAB 示例”把整个问题压缩到了一个可读、可改、能跑出流场图的工程里：它用流函数－涡量法离散不可压缩流动的控制方程，通过时间步进模拟顶盖驱动方腔流这一类经典算例，输出速度场、涡量场和流函数分布。你不用理解全部数值方法细节，照着改雷诺数、改网格密度，就能看到涡旋从生成到发展的全过程。适合正在选做课程设计、刚接手流体仿真任务、或者想找一个干净底子做二次开发的 MATLAB 使用者。我把它完整拆了一遍，从公式到代码行逐一讲清楚。

2. 数学模型与离散方案：为什么选流函数－涡量法而不是直接解压力

2.1 非定常不可压缩 N-S 方程的完整形式

在二维直角坐标系下，不可压缩非定常 Navier-Stokes 方程组的完整形式包含三个分量方程：两个动量方程和一个连续性方程。速度分量 u 和 v 分别沿 x 和 y 方向，p 为压力，ν 为运动黏度：

% 动量方程（x 方向） % du/dt + u*du/dx + v*du/dy = -1/rho * dp/dx + nu * (d2u/dx2 + d2u/dy2) % 动量方程（y 方向） % dv/dt + u*dv/dx + v*dv/dy = -1/rho * dp/dy + nu * (d2v/dx2 + d2v/dy2) % 连续性方程 % du/dx + dv/dy = 0

直接数值求解这套原始变量（u,v,p）的非定常方程会遇到一个经典麻烦：连续性方程没有显式的时间演化形式，压力场的求解需要反复迭代修正，程序写出来又复杂又容易发散。这个示例采用的处理方式是引入流函数 ψ 和涡量 ω，把问题从“三个方程耦合求解”降维成“两个方程分步求解”，稳定性也更好控制。这是上世纪七八十年代发展起来、至今仍广泛用于教学和工程评估的经典方案。

2.2 非定常涡量输运方程与运动学关系的离散

通过取旋度消去压力梯度项，得到非定常涡量输运方程：

% 涡量输运方程（非定常核心） % dω/dt + u*dω/dx + v*dω/dy = nu * (d2ω/dx2 + d2ω/dy2) % 其中 ω = dv/dx - du/dy % 而流函数满足泊松方程：d2ψ/dx2 + d2ψ/dy2 = -ω % 速度可由流函数恢复：u = dψ/dy，v = -dψ/dx

方程里每一项的物理含义值得说清楚：左边第一项是涡量随时间的变化率，即非定常项，这是本示例区别于定常求解器的核心；第二、三项是涡量的对流输运，它让流场中的涡旋结构随时间移动和变形；右边是分子黏性造成的涡量扩散。整个时间推进的循环就是：先根据当前时刻的速度场计算涡量的对流和扩散增量，更新得到下一时刻的涡量场；再由更新后的涡量场求解泊松方程，得到流函数场；最后从流函数差分出速度场，进入下一步。每一步之间环环相扣，是理解这段代码的主线。

空间离散采用二阶中心差分，时间推进采用显式 Euler——虽然精度只到一阶，但胜在实现简单、代码短，符合示例定位。在实际做项目时，我会把时间推进换成 Adams-Bashforth 二阶格式，这个在后面“改造”小节给出具体做法。对于显式格式，需要注意步长限制：时间步增量不能超过网格尺寸和对流速度共同决定的 CFL 上限，也不能超过黏性扩散的稳定极限，否则算到一半流场就炸了。

3. 程序模块拆解：从参数初始化到流场可视化的完整链路

3.1 工程目录与文件组织

这个示例代码设计得比较轻量，整体目录结构大致如下：

NS2D_Unsteady/ ├── main.m % 主程序：参数设置、时间循环、后处理 ├── initField.m % 初始化流函数、涡量、速度场 ├── solvePoisson.m % 泊松方程求解（流函数） ├── vorticityUpdate.m % 涡量输运方程的时间推进 ├── boundaryCondition.m % 边界涡量与流函数的约束 └── plotResult.m % 可视化输出

主程序负责装配整个计算流程，各子函数按职责拆分。这种结构最大的优点是把“数值方法”和“数据流程”分开了：改边界条件只需要动 boundaryCondition.m，改时间推进格式只需要动 vorticityUpdate.m，不需要在一大段顺序代码里寻找修改点。如果只是课程作业，把所有子函数并进一个 main.m 也能跑，但二次开发时你会后悔没有遵守这个模块边界。

3.2 参数初始化与网格生成

主程序开头是参数设置段。以顶盖驱动方腔流为例，方腔边长为 L，网格等分为 nx×ny 个单元，雷诺数 Re 是核心无量纲参数：

% 顶盖驱动方腔流参数设置 L = 1.0; % 方腔边长 nx = 64; ny = 64; % 网格点数 Re = 100; % 雷诺数（定义：Re = U_top * L / nu） nu = U_top * L / Re; % 运动黏度 U_top = 1.0; % 顶盖运动速度 dt = 0.001; % 时间步长 Tmax = 20; % 总模拟时长 dx = L / (nx-1); % 空间步长 dy = L / (ny-1); % 空间步长 x = linspace(0, L, nx); y = linspace(0, L, ny);

雷诺数直接控制流态：Re=100 时流场有一个稳定的主涡，Re=1000 时角落会出现次级涡，Re=5000 以上流动变得非定常甚至混乱。这个参数是研究流动状态最重要的旋钮，改它比改网格更直观地影响流场结构。初学者很容易把这里的运动黏度计算漏掉，直接写 nu=0.01 之类随意数值，这会让“雷诺数”名存实亡，后续和文献结果对比就失去了依据。

3.3 涡量输运方程的时间推进实现

涡量更新的核心代码不长，但集中了这个示例的数值方法精华。显式时间推进的意思很直接：用当前时刻的涡量空间导数去估计时间变化率，然后外推一个时间步长。这里给出 vorticityUpdate.m 的关键循环体：

% 涡量输运方程显式时间推进（以二维数组方式实现，循环速度远快于双重 for） function omega_new = vorticityUpdate(omega, psi, dx, dy, dt, nu, nx, ny) % 由流函数计算速度场：u = dpsi/dy，v = -dpsi/dx，中心差分 u = zeros(nx, ny); v = zeros(nx, ny); u(2:end-1, 2:end-1) = (psi(2:end-1, 3:end) - psi(2:end-1, 1:end-2)) / (2*dy); v(2:end-1, 2:end-1) = -(psi(3:end, 2:end-1) - psi(1:end-2, 2:end-1)) / (2*dx); % 涡量的空间导数（中心差分） dOmega_dx = (omega(2:end-1, 3:end) - omega(2:end-1, 1:end-2)) / (2*dx); dOmega_dy = (omega(3:end, 2:end-1) - omega(1:end-2, 2:end-1)) / (2*dy); % 拉普拉斯算子（二阶中心差分） lapOmega = (omega(2:end-1, 3:end) - 2*omega(2:end-1, 2:end-1) + omega(2:end-1, 1:end-2)) / dx^2 ... + (omega(3:end, 2:end-1) - 2*omega(2:end-1, 2:end-1) + omega(1:end-2, 2:end-1)) / dy^2; % 时间推进：dω/dt = -u*dω/dx - v*dω/dy + nu * lapOmega omega_new = zeros(nx, ny); omega_new(2:end-1, 2:end-1) = omega(2:end-1, 2:end-1) + dt * ( ... -u(2:end-1, 2:end-1) .* dOmega_dx ... -v(2:end-1, 2:end-1) .* dOmega_dy ... + nu * lapOmega); end

参数含义对应地解释一下：u 和 v 是当前时刻速度场，由流函数 ψ 差分得到；dOmega_dx、dOmega_dy 是对流项的一阶导数；lapOmega 是黏性扩散项的二阶导数。时间步进公式中，对流项是“涡量被流场携带运动”的数学表达，扩散项是“涡量从高浓度向低浓度抹平”的数学表达。代码里使用 MATLAB 矩阵整体运算而不是 for 循环嵌套，是一个很实用的性能技巧：同样的网格规模，向量化写法一般比嵌套循环快一个数量级，时间循环总共要跑上万步，这个差距直接决定你是喝一杯咖啡还是等一晚。

3.4 泊松方程求解与边界条件

每个时间步内，拿到新的涡量场后需要求解流函数泊松方程。示例直接调用 MATLAB 的反斜杠算子求解稀疏线性方程组：

% 泊松方程：d2(psi)/dx2 + d2(psi)/dy2 = -omega % 组装五点差分格式的系数矩阵 A N = (nx-2) * (ny-2); A = sparse(N, N); rhs = zeros(N, 1); for j = 2:ny-1 for i = 2:nx-1 idx = (j-2)*(nx-2) + (i-1); A(idx, idx) = -2/dx^2 - 2/dy^2; if i > 2, A(idx, idx-1) = 1/dx^2; end if i < nx-1, A(idx, idx+1) = 1/dx^2; end if j > 2, A(idx, idx-(nx-2)) = 1/dy^2; end if j < ny-1, A(idx, idx+(nx-2)) = 1/dy^2; end rhs(idx) = -omega(i, j); end end psi_vec = A \ rhs; % 稀疏矩阵求解 psi = zeros(nx, ny); psi(2:end-1, 2:end-1) = reshape(psi_vec, nx-2, ny-2)';

A 矩阵每一行对应方腔内部一个网格点的离散方程：对角线是 -2/dx²-2/dy²，四邻域各有权重。方程右侧 -ω 就是涡量对流的“源项”，物理意义是涡旋结构决定了流动的形态。用反斜杠求解的好处是 MATLAB 自动选用稀疏直接法或迭代法，对 64×64 这类网格在单时间步内达到毫秒级延迟。需要知道的一个局限是：矩阵组装用了两层 for 循环，只在初始化时执行一次，所以即使写得粗糙也不会成为性能瓶颈；真正决定程序快慢的是时间循环里的涡量更新段。

边界条件是这个示例中最容易写错、也最能体现数值功底的地方。顶盖驱动流对 ψ 和 ω 的约束如下：

% 边界条件设置 % 顶盖（y = L）：u = U_top，v = 0，dpsi/dy = U_top % 底面、左壁、右壁：无滑移，psi = 0（即不穿透） % 边界涡量由流函数导数反推： % 左壁：omega(1,j) = -2 * psi(2,j) / dx^2 % 右壁：omega(nx,j) = -2 * psi(nx-1,j) / dx^2 % 底面：omega(i,1) = -2 * psi(i,2) / dy^2 % 顶盖：omega(i,ny) = -2 * psi(i,ny-1) / dy^2 - 2*U_top/dy

边界涡量的公式来自“无滑移边界下单侧导数近似”：壁面涡量正比于紧邻壁面的流函数值。顶盖那一项里多出来的 -2U_top/dy，是因为顶盖在运动，剪切层内部顶着有限速度差。很多人在这一步忽略顶盖项，结果就是流场整体形状对，但顶盖附近的涡量峰值明显偏小，与文献基准解一比就露馅了。

3.5 可视化输出

示例自带绘图段，输出三个子图：涡量云图、流函数等值线、速度矢量图。绘图的代码比较常规：

% 可视化当前时刻流场 figure(1) subplot(1, 3, 1) contourf(x, y, omega', 20); colorbar; title('涡量场') subplot(1, 3, 2) contour(x, y, psi', 20); title('流函数') subplot(1, 3, 3) quiver(x(1:4:end), y(1:4:end), u(1:4:end, 1:4:end)', v(1:4:end, 1:4:end)') title('速度场') drawnow

可视化不仅可以用来“看个大概”，还有实际价值：观察涡量云图是否对称、流函数等值线是否光滑，能快速判断边界条件有没有写对。绘图本身只消耗很少的时间比例，建议每几百步输出一次，而不要每一步都画图——否则 matplotlib 式的地帧重绘会让时间循环显著变慢。实际做项目时，我会把特定时刻的完整流场数据存成 .mat 文件，后处理单独做，避免反复重算。

4. 运行指南与算例验证：把命令跑通，再把结果跑对

4.1 从零到第一个流场图的完整步骤

拿到源码后的标准操作流程如下：在 MATLAB 里打开 main.m，直接 F5 运行。无需安装任何工具箱，稀疏矩阵反斜杠、contourf、quiver 这些都是 MATLAB 基础功能。文件里如果包含子函数，把它们全部加入当前工作目录或统一放在一个文件夹中。运行过程中命令行窗口会输出当前时间步和模拟时间，几秒后弹出流场图。

% 在 MATLAB 命令行中运行 >> main % 或者直接点击编辑器中的“运行”按钮

建议第一次跑用默认参数：64×64 网格、Re=100、总模拟时长 Tmax=20。Re=100 属于层流稳态工况，流场在 5~10 个无量纲时间单位后就基本稳定，你会看到一个位于方腔中央偏右的大涡。这时把 Re 提高到 1000，缩短 Tmax 到 10 左右，观察涡心偏移和角涡出现的规律。如果出现 NaN 或者数值发散，优先减小 dt。

4.2 如何验证结果与经典基准解对比

跑通不等于跑对。几乎所有流体数值代码都需要和经典基准解对照，否则改了一个参数后结果逐渐失真你却浑然不知。有两类最常用的对比方式。第一类是流场特征对照：方腔中心剖面上 u 速度沿竖直方向的分布曲线，与文献中的基准数据对比，Re=100 时曲线应当重合；第二类是涡量极值对照，顶盖处的涡量最大值随雷诺数的变化规律。后者尤其适合检查边界涡量公式是否实现正确。

% 提取竖直中线上的 u 速度分布用于结果对比 midLine = round(nx/2); uVertical = u(midLine, :); plot(uVertical, y); % 应近似关于 y 对称的 S 形曲线 hold on; % 此处可加载文献基准数据 plot(ref_u, ref_y, 'ro')

需要特别说明的是：不要期望曲线与文献基准解完全重合。网格分辨率 64×64 加上一阶时间精度，中心线上尾部的误差在 1% 到 3% 之间属于正常范围；误差偏大且曲线形状不一，才需要考虑数值问题。

4.3 参数扫描方法与实验设计

研究流动的一个常见需求是“看不同雷诺数下的流场演化”。写一个外层循环自动扫描 Re 值，并把各个算例的结果保存成独立文件即可。另一个实用的实验设计是“固定 Re、加密网格”，观察结果是否收敛——如果 Re=100 下 32×32 和 128×128 网格得到的速度曲线差异很大，说明你的网格还没达到网格无关性，后续所有结果都要打一个问号。这个“网格无关性验证”是论文和报告中必须呈现的核心内容，示例代码同样支持这种用法：直接改 nx、ny 重新运行即可。

5. 常见问题排查：关于发散、边界涡量和可视化异常的排雷日志

5.1 时间步进几十步后流场出现 NaN

现象：程序运行 20 到 50 个时间步后，涡量场突然出现 NaN，随后整个流场崩溃，绘图窗口空白或显示错误。原因分析：显式格式的时间步长超出了稳定极限。涡量输运方程包含对流项和扩散项，对流项的稳定性要求 dt ≤ dx/|u|max，扩散项要求 dt ≤ dx²/(4ν)。实际程序中报错前的几万步通常都是“潜伏期”，网格越细黏度越低，潜伏期越短。解决办法也很直接：把 dt 从 0.001 调整到 0.0005 或更低；或者更系统地，写一个自动计算最大稳定步长的函数。有经验之后我习惯把 dt 取三个限制量中的最小值再乘 0.5 的保险系数，彻底告别这类问题。

5.2 边界附近涡量云图出现锯齿状等值线

现象：涡量云图整体轮廓正常，但紧贴壁面处等值线不再光滑，出现锯齿纹理。原因分析：边界涡量公式中的单侧差分只有一阶精度，内部点使用二阶中心差分，这种精度不匹配会在边界引起虚假振荡，用高分辨率云图看得特别清楚。解决办法：用二阶精度的单侧公式替代一阶公式。例如左壁的涡量计算公式从 -2ψ(2,j)/dx² 改写为 (-3ψ(1,j) + 4ψ(2,j) - ψ(3,j)) / (2dx²)。改完边界涡量，等值线明显光滑。这里需要提醒一个反直觉的事实：边界涡量的精度直接决定整体求解精度，不是边界内部各占一半，而是边界负责给内部提供“锚点”，锚歪了内部再准也无济于事。

5.3 速度矢量图看起来方向和预期相反

现象：顶盖驱动流中，靠近顶盖的水平速度方向与顶盖运动方向相反，流场存在一个明显的大涡但旋转方向反了。原因分析：流函数与速度的关系符号约定不一致，常见的是 ψ 的 y 向导数符号弄反。解决办法：检查速度恢复公式，正确关系是 u = ∂ψ/∂y、v = -∂ψ/∂x。这里有一个查验方法：方腔内一个大涡对应 ψ 呈山峰状分布，沿 x 方向看 ψ 上升、沿 y 方向看 ψ 下降，就能推断出正确的速度分量方向。这类错误与雷洛数、网格无关，是纯粹的符号约定问题，思想方法最有效。

5.4 高位 Re 下程序越跑越慢、内存占用飙升

现象：把 Re 从 100 改到 10000 后，程序快要完成时突然变得极慢，甚至报告内存不足。原因分析：高频流动出现小尺度涡结构，在固定网格下 CFL 要求时间步大幅缩小，总步数翻倍。而 MATLAB 每次循环如果意外产生了完整矩阵的拼接或动态增长，会出现内存碎片化。解决办法：提前用 zeros 初始化所有大型数组；如果确实需要高频非定常计算，则要考虑加密网格、减少总模拟时长、或将结果降采样保存。需要说明的是，显式方法的扩散项对时间步长的限制与 dt ≤ dx²/(2ν) 相关，即细网格、低黏度下时间步会急剧压缩，这是显式方法的固有局限，改用隐式处理扩散项可以绕过这个坎，但实现复杂度也上来了。

6. 二次开发实战：把“能跑的示例”改造成“能研究问题的小型求解器”

先说第一个建议：动手改之前，先跑一遍原版代码并对照标准算例验证，确认你手上的版本是“基线正确”的。任何后续改造都要基于这个基线，否则会遭遇“改了十处代码和跑出一个不知道哪里的错”的混乱局面。基线验证通过后，按由浅到深的顺序选择改造点：

第一步，把时间推进从显式 Euler 升级为二阶 Adams-Bashforth。做法是保存前一步的涡量增量，当前增量与上一步增量做加权组合：ω(n+1) = ω(n) + dt×(1.5×F(n) - 0.5×F(n-1))，其中 F 表示涡量输运方程的右端。这个改动只涉及循环内的几行代码，但时间精度从一阶升到二阶，长期演化更稳，而且为后续加源项提供了干净入口。这是“性价比”最高的一处改造。

第二步，改造几何模型。原示例默认方腔流，修改壁面形状需要更换网格生成和边界条件两处。比如改成圆柱绕流，抛弃笛卡尔网格是更好的选择，整体框架才换得更彻底；但改成长宽比不同的矩形通道流则只需调整 L 和 nx、ny，几乎零成本。这类保守改造恰好适合练手。

第三步，把压力场恢复出来。流函数-涡量法算完不直接给压力场，研究受力就需要恢复压力。做法是对动量方程取散度得到压力泊松方程，边界用诺伊曼条件，这与原示例的泊松求解套路一致，可以直接复用稀疏矩阵框架。恢复压力后，可以计算圆柱或翼型的升阻力系数，后续再研究降噪、受力分析等工程问题才有根基。

第四步，把算例参数文件独立出来。将 Re、网格尺寸、求解时长、输出频率全部放进一个结构体，配合循环脚本自动批量跑参数工况，批量后处理时效率极高。我现在做多工况对比都是这套流程：先把单例结果验证到和文献数据对得上，再放开参数扫描，最后统一出图和分析。从那以后，我每次改造求解器都强制走一遍“基线复现到确认没问题再动刀”的流程，改一个参数就重新观察一次流场演化，杜绝原地大改。希望你也能把求解器的每一次改动变成一次可以回滚、可以对照的实验记录，这比灵光一现的“改完就跑”高效得多。希望帮到你。

本文还有配套的精品资源，点击获取