简介：二维非定常 Navier-Stokes 方程的 MATLAB 实现示例，面向流体力学初学者、CFD 课程学生以及需要开展流场仿真的工程技术人员。资源以有限差分法与压力修正算法为主线，覆盖网格生成、边界条件施加、时间推进、泊松方程求解等典型环节，并通过多种辅助函数构建出清晰的模块化程序框架，便于逐段对照学习。压缩包共七十六个文件，包括四十四份 MATLAB 源程序、二十九张过程结果图、一份许可文本和一份超文本说明文档，整包大小约四百六十三 KB，轻巧易取。二十九张图像展示不同时刻的速度场与压力场演化，超文本文档辅助解释代码结构，帮助读者快速掌握非定常流动模拟的关键步骤。已有约一千九百四十一人学习浏览，适合用来巩固有限差分、压力修正等数值方法，并可在修改边界条件或参数后拓展至更多流动场景。

1. 一个能跑的MATLAB二维非定常Navier-Stokes示例：从哪来、能算啥、坑在哪

如果你在找的是那种“粘过去就能出图”的MATLAB二维非定常Navier-Stokes示例，那我先把结论说清楚：这份资源跑通之后，你能得到的是方腔顶盖驱动流从启动到定常的完整演化过程，以及一个自带解析解的衰减涡算例。它不算新，用的也是很多人觉得“过时”的有限差分加投影分步法，但恰恰是这套朴素的组合，把非定常不可压缩流动里最难啃的压力耦合问题拆成了三步，适合拿来当课程设计、毕业论文前期的验证基准。适合谁呢？正在学计算流体力学、需要在MATLAB里从零搭NS求解器、或者是想快速验证自己离散格式有没有写错的人。我见过不少人下载后跑出NaN就放弃了，绝大多数不是代码问题，而是CFL条件、边界索引、迭代容差这三个地方没对齐。这份示例真正值钱的不是那几行MATLAB代码，而是它替你标好了所有需要对齐的坐标。

2. 从连续方程到离散格式：交错网格与CFL条件怎么定

2.1 为什么选“不可压缩+二维+显式时间推进”这个组合

这份MATLAB示例求解的是二维不可压缩Navier-Stokes方程，连续方程是速度散度为零，动量方程是经典的对流项、压力梯度项和粘性项。二维意味着只有u、v两个速度分量，外加一个压力p。不可压缩的难点不在于未知数多，而在于压力没有独立的演化方程，它是以“散度为零”这个约束形式出现。直接同时解三个场会遇到鞍点问题，常规做法都是把压力处理成通过泊松方程求解的拉格朗日乘子。

这套示例选的是显式时间推进，也就是每个新时刻的速度都由上一时刻的已知量直接算出来。显式的好处是每步代价低、代码和直觉一致，坏处是时间步长受稳定性约束。二维低雷诺数流场，显式完全够用；如果你以后要跑高雷诺数湍流，再考虑隐式或压力投影的隐式版本也不迟。

方程形式如下，非定常体现在左边第一项时间导数：

$$\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)$$

代码里不会直接写成这种符号形式，而是拆成有限差分。这里的关键是：数值上满足连续方程的速度场，和满足动量方程的速度场往往不一致，所以需要投影分步法在每一步里“修正”一次。

2.2 交错网格：为什么压力和速度不能住在同一个格点上

如果让u、v、p都定义在同一套网格节点上，压力泊松方程的中心差分会允许一种棋盘状的振荡压力场存在：正负压力交替排列，但梯度为零，导致动量方程根本感受不到它，这种“伪模态”在迭代解泊松方程时特别容易被放大。使用交错网格是常见做法，也就是把速度分量放在网格单元边界上，压力放在单元中心。

这份示例里典型的数组尺寸是这样分配的，如果你看到u、v、p三个矩阵尺寸完全一样，那就要警惕棋盘振荡：

% 计算域与网格参数 Lx = 1; % x方向长度 Ly = 1; % y方向长度 nx = 64; % x方向控制体数量 ny = 64; % y方向控制体数量 dx = Lx / nx; dy = Ly / ny; % 交错网格数组 % p 定义在单元中心，尺寸为 nx x ny % u 定义在竖直单元边上，尺寸为 (nx+1) x ny % v 定义在水平单元边上，尺寸为 nx x (ny+1) u = zeros(nx+1, ny); % 水平速度分量 v = zeros(nx, ny+1); % 竖直速度分量 p = zeros(nx, ny); % 压力

这段初始化把网格数、计算域尺寸和三个场的存放位置全部框死了。网格坐标写法是x_c = (0.5:nx)*dx，意思是压力节点位于每个单元的中心，速度分量相应地错开半个网格。索引错位是初学交错网格最容易翻车的地方，建议在写迭代代码前先画一张小网格，把 u(i,j)、v(i,j)、p(i,j) 各自对应的物理位置标出来，再动手写循环。参数 Lx、Ly、nx、ny 决定了空间分辨率，Re 决定流动特性，这四个参数在后面的算例里会反复调整。

2.3 dt 不是随便取的：CFL 与粘性稳定性的合并约束

显式时间推进有两个稳定性约束必须同时满足。第一个是对流CFL条件，信息在一个时间步内不能跨过一个网格单元，量化后是u*dt/dx和v*dt/dy都必须小于某个安全系数，通常取0.25到0.5。第二个是显式扩散项的稳定性约束，这个容易被忽略：粘性项的显式处理要求2*dt/(Re*dx^2)小于1，也就是dt < Re*dx^2/2。很多人只盯着CFL，结果低雷诺数、细网格时照样发散。

综合判断时，我一般会写成一个动态计算dt的代码，而不是手工试：

% 动态时间步长：同时满足对流与粘性稳定性 umax = max(abs(u(:))); vmax = max(abs(v(:))); % 对流限制 dt_cfl = min(dx/umax, dy/vmax) * 0.25; % 粘性限制（显式扩散稳定约束） dt_diff = 0.5 * Re * min(dx^2, dy^2); % 实际采用两者中较小的 dt = min(dt_cfl, dt_diff);

这里的0.25是安全系数，如果流场中出现局部高速区，这个系数要再降低。dt_diff的公式来自一维扩散方程的von Neumann稳定性分析，二维时取两个方向里更严格的那个即可。算完之后最好把dt打印出来看一眼，如果比初始猜测值小一个数量级，那是正常的，别硬顶着大时间步跑，会直接爆掉。

3. 投影法三步走：压力泊松方程的MATLAB求解细节

3.1 第一步：忽略压力梯度，推进中间速度场

投影法的第一步是先去掉动量方程里的压力梯度项，用当前时刻的u、v计算一个新的中间速度场，通常记为u星。这一步物理上不完整，因为还没满足连续方程，但它已经包含了对流和粘性的贡献。

动量方程去掉压力项后，u的预测方程为：

$$u^* = u^n + dt \left[ -\left(u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y}\right) + \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \right]$$

MATLAB里的核心循环长这样：

% 预测步：计算中间速度 ustar、vstar % 以 u 分量为例，v 分量对称处理 % 粘性项：二阶中心差分 uxx = (u(3:nx+1, :) - 2*u(2:nx, :) + u(1:nx-1, :)) / dx^2; uyy = (u(:, 2:ny) - 2*u(:, 1:ny-1) + u(:, 1:ny-1)) / dy^2; % 注意边界索引 visc_u = (uxx + uyy) / Re; % 对流项：一阶迎风，先判断速度方向 % u*w 的乘积在交错网格上需要插值到同一位置 % 这里用简单的一阶迎风示意，完整代码里会做半格插值 conv_u = zeros(nx+1, ny); for i = 2:nx for j = 2:ny-1 ue = 0.5*(u(i+1,j) + u(i,j)); uw = 0.5*(u(i,j) + u(i-1,j)); un = 0.5*(u(i,j+1) + u(i,j)); us = 0.5*(u(i,j) + u(i,j-1)); conv_u(i,j) = ue*max(ue,0) - uw*max(uw,0) ... + un*max(un,0) - us*max(us,0); end end % 中间速度 ustar(2:nx, 2:ny-1) = u(2:nx, 2:ny-1) + dt*(-conv_u(2:nx,2:ny-1) + visc_u);

这段代码里最需要注意的就是交错网格上的乘积插值：对流项u*du/dx里的速度相乘，必须在同一位置取值，否则会产生数值不稳定。注释里“半格插值”的意思是相邻两个速度的平均。一阶迎风是低雷诺数算例的稳妥选择，虽然耗散偏大，但不会像中心差分那样在间断处产生振荡。如果你想要二阶精度，可以把迎风换成二阶迎风或者Quick格式，代价是实现复杂度上升。

3.2 第二步：解压力泊松方程，这一步的坑最多

中间速度u星不满足散度为零，我们要找一个压力场，用它的梯度去修正速度，使得修正后的速度场无散度。推导过程很经典：对修正公式两边取散度，要求div(u^{n+1})=0，就得到压力泊松方程：

$$\nabla^2 p = \frac{1}{dt} \nabla \cdot u^*$$

这个方程在代码里不会去组装大型稀疏矩阵，而是用Jacobi迭代逐点更新。原因有两个：一是迭代法在64乘64网格上收敛很快，代码也直观；二是避免在MATLAB里维护稀疏矩阵存储。代价是迭代次数和收敛容差需要手工控制。核心代码：

% 计算散度 RHS div_u = zeros(nx, ny); div_u(2:nx-1, 2:ny-1) = (ustar(3:nx, 2:ny-1) - ustar(2:nx-1, 2:ny-1)) / dx ... + (vstar(2:nx-1, 3:ny) - vstar(2:nx-1, 2:ny-1)) / dy; rhs = div_u / dt; % Jacobi 迭代求解压力泊松方程 p = zeros(nx, ny); tol = 1e-6; maxiter = 2000; for iter = 1:maxiter p_old = p; % 内部节点五点差分 p(2:nx-1, 2:ny-1) = 0.25 * ( ... p_old(1:nx-2, 2:ny-1) + p_old(3:nx, 2:ny-1) + ... p_old(2:nx-1, 1:ny-2) + p_old(2:nx-1, 3:ny) - ... dx^2 * rhs(2:nx-1, 2:ny-1) ); % 边界：Neumann 零梯度条件 p(1,:) = p(2,:); p(nx,:) = p(nx-1,:); p(:,1) = p(:,2); p(:,ny) = p(:,ny-1); if max(abs(p(:) - p_old(:))) < tol break; end end

rhs的计算是整个投影法的心脏，这里散度差分必须和后面速度修正里的压力梯度差分使用完全相同的网格位置。如果散度用(u(i+1)-u(i))/dx，梯度也必须用对应的交错位置，混用会导致修正后仍然有残余散度。Jacobi迭代的好处是天然满足Neumann边界条件，边界处直接让边界压力等于相邻内点压力。tol取1e-6在大多数算例够了，如果你观察流场有缓慢的周期波动，可以调到1e-8。

3.3 第三步：修正速度，让散度严格为零

得到压力场后，修正步用一个简单的梯度减法完成：

$$u^{n+1} = u^* - dt \frac{\partial p}{\partial x}, \quad v^{n+1} = v^* - dt \frac{\partial p}{\partial y}$$

这一步做完之后，理论上速度场就满足连续方程了。MATLAB实现很短：

% 修正步：压力梯度修正速度场 u(2:nx, 2:ny-1) = ustar(2:nx, 2:ny-1) - dt / dx * (p(2:nx, 2:ny-1) - p(1:nx-1, 2:ny-1)); v(2:nx-1, 2:ny) = vstar(2:nx-1, 2:ny) - dt / dy * (p(2:nx-1, 2:ny) - p(2:nx-1, 1:ny-1)); % 主时间循环骨架 for n = 1:nt % 1. 边界条件施加（固壁无滑移或解析解边界） apply_boundary_conditions(u, v); % 2. 预测步得到 ustar, vstar % 3. 求解压力泊松方程 % 4. 修正速度场 % 5. 记录输出（速度范数、残差、动画） end

每步循环顺序不能乱：先施加边界，再预测、投影、修正。如果边界条件在预测步之后才施加，边界上会积累误差。非定常的意义在于主循环里的每一步都对应真实的物理时间，每一个中间场都是物理上可解释的流动状态。衰减涡算例验证的正是这一点：同一时刻的数值解必须无限逼近解析解，而不是只求“最终定常结果正确”。

4. 边界条件与两个标准算例：把代码真正跑起来

4.1 顶盖驱动流：非定常推进到定常的经典路径

顶盖驱动流是所有NS代码的标配验证算例。计算域是单位方腔，顶盖沿x方向以恒定速度运动，其余三壁固定。雷诺数的定义是Re = U0 * L / nu，其中U0是顶盖速度，L是方腔边长。这份示例里就取U0等于1、L等于1，所以Re = 1/nu，代码里直接通过Re控制粘性系数大小。

边界条件的实现是一行行硬编码的：

% 顶盖驱动流边界 % 下壁、左壁、右壁：无滑移 u=0, v=0 u(:, 1) = 0; % 下壁 u(:, ny) = 1; % 顶盖：水平速度设为 1 v(:, 1) = 0; % 下壁 v(:, ny) = 0; % 顶盖：竖直速度为零 % 左壁与右壁 u(1, :) = 0; u(nx+1, :) = 0; v(1, :) = 0; v(nx+1, :) = 0;

跑顶盖驱动流时有个容易误解的点：题目叫非定常，但顶盖驱动流的经典验证目标是“定常解”。启动瞬间整个流场是一潭死水的，顶盖运动后涡逐渐发展，经历一个过渡过程，最终达到定常状态。所以这份示例既能展示非定常演化，又能通过定常终态验证收敛性。判断定常别用“跑够多少步”，而是看速度场的变化量：

% 定常收敛判据：连续两步速度场变化小于阈值 diff_u = max(abs(u(:) - u_prev(:))); diff_v = max(abs(v(:) - v_prev(:))); if max(diff_u, diff_v) < 1e-8 break; end

Re=100时，粗网格跑到定常大概需要几千步，细网格上万步。如果Re提高到1000，流场会呈现出更复杂的二次涡结构，收敛也更慢，这也是后续练手的好对象。

4.2 衰减涡算例：用解析解自校验时间精度

顶盖驱动流只能验证定常结果对不对，验证非定常时间积分的正确性需要另一个算例：衰减涡。它有一个解析解，在方腔内用正弦函数构造一个涡，其强度随时间指数衰减。这个算例的最大价值是能把每个时刻的数值解和精确解对比，直接算出误差阶。

衰减涡的解析速度场为：

% 衰减涡初始场：正弦涡 % 解析解：u = sin(pi*x) * cos(pi*y) * exp(-2*pi^2*t/Re) % v = -cos(pi*x) * sin(pi*y) * exp(-2*pi^2*t/Re) x_c = ((1:nx)-0.5) * dx; y_c = ((1:ny)-0.5) * dy; [Xc, Yc] = meshgrid(x_c, y_c); % 初始时刻 t=0 u = sin(pi*Xc) .* cos(pi*Yc); v = -cos(pi*Xc) .* sin(pi*Yc);

注意初始场给的是单元中心位置，而u、v数组在交错网格上分别位于不同位置，所以要对x、y坐标做适当的偏移，这又是一个索引细节。误差计算在任意时刻做：

% 在给定时刻对比数值解与解析解 t = n * dt; u_exact = sin(pi*Xc) .* cos(pi*Yc) .* exp(-2*pi^2*t/Re); err_u = sqrt(sum((u(2:nx+1, :) - u_exact(:,:)).^2) / numel(u_exact));

这里我习惯用相对L2误差来判断，而不是看最大误差，因为最大误差对个别点太敏感。跑衰减涡时，边界条件不能用固壁无滑移，而应该把解析解直接赋在边界上，这样误差只来自内部离散，不会被边界处理污染。

4.3 后处理：看见流场才算真正跑通

跑完计算后至少看三个图：速度矢量图、压力云图、流线图。MATLAB里没有特别多技巧，核心命令就这几个：

% 视觉化流场 figure quiver(x_u, y_u, u', v', 2); % 速度矢量，2是缩放比例 hold on contourf(x_c, y_c, p', 20); % 压力云图 axis equal title('二维非定常方腔流瞬时流场'); % 导出数据便于后续分析 save('cavity_result.mat', 'u', 'v', 'p', 't', 'Re', 'nx', 'ny', 'dx', 'dy', 'dt');

速度矢量图能直观看到主涡和角涡，压力云图能看出高低压区分布，流线图适合展示涡心位置。保存为后缀为mat的文件后，之后再想画中心线速度剖面、算涡量、做动画都不用重跑一遍计算，这点很重要，因为在细网格上重跑一次的代价不小。

5. 避坑指南：发散、不收敛与伪非定常的五个典型现场

5.1 前几步就爆NaN：十有八九是dt超限

现象：程序启动后前几个时间步速度场就出现Inf或NaN，绘图时报错，整个计算直接终止。

原因：时间步太长，超过了上一章说的对流或粘性稳定性边界。尤其是低雷诺数加细网格时，粘性约束项Re*dx^2/2往往比CFL约束更严格，只按CFL取dt必爆。另一个常见原因是初始速度场里存在不合理的局部高速，动态dt计算时被低估。

解决：把dt计算改成动态的，每步都重新扫描全场最大速度，并同时套用对流和粘性两个约束。如果仍然爆炸，把安全系数从0.25降到0.1试一次。这个问题的特征是“换个dt就从爆变稳”，不用怀疑代码逻辑。

5.2 压力场棋盘振荡：检查你的网格是不是“真交错”

现象：速度场看起来大体合理，但压力云图呈密集的红黑棋盘格交替，数值上相邻格点压力正负交替出现。

原因：压力棋盘振荡是泊松方程零能模态的表现，通常因为u、v、p数组尺寸相同，也就是名义上“交错”实际没交错；或者是散度和梯度算子用的差分方向不一致，导致压力修正没有作用到速度的正确位置。

解决：先确认数组尺寸满足size(p)=[nx,ny]、size(u)=[nx+1,ny]、size(v)=[nx,ny+1]。再检查修正步里的压力梯度索引，u(2:nx,...)对应p(2:nx,...)-p(1:nx-1,...)，如果写成了p(3:nx+1,...)-p(2:nx,...)，压力梯度会整体平移一格，相当于用错了离散模板。这两点检查完，棋盘振荡基本消失。

5.3 顶盖驱动流左上角压力爆表：角点奇点处理不当

现象：顶盖驱动流跑到几百步后，左上角和右上角附近的压力异常大，速度等值线在角点处剧烈扭曲。

原因：顶盖水平速度是1，左壁水平速度是0，这两个边界条件在角点直接冲突，形成了数学上的奇点。网格越细，角点处的速度梯度和压力越尖锐。

解决：在角点处对顶盖速度做一个平滑渐变，常见做法是把角点速度设为零，或者用一段余弦渐变过渡。代码里就是一行的事：

u(nx+1, ny) = 0; % 右上角速度置零，消除奇点 u(1, ny) = 0; % 左上角速度置零

加这一行之后，全局压力最大值会明显回落。这个改动只影响角点附近一个网格的精度，对整体流场影响很小。

5.4 衰减涡误差不降反升：边界条件施加错了

现象：网格从32加密到64再到128，误差几乎不变甚至在增大，完全不符合二阶格式应有的表现。

原因：衰减涡算例的边界上，如果把固壁条件当边界用，边界附近的数值解会和解析解完全脱节，误差由边界主导，加密网格也救不回来。另一个原因是时间步长没同步缩小，空间二阶但时间一阶，时间误差盖过了空间误差。

解决：衰减涡的所有边界都直接赋解析解的值，用当前位置和当前时刻的精确公式计算边界上的u、v。同时把dt减半和网格加密两件事分开做，先固定dt只加密网格测空间精度，再固定网格只减dt测时间精度，就能定位到底是哪一阶拖了后腿。

5.5 看着在动其实是伪非定常：收敛判据太乐观

现象：顶盖驱动流跑了预设步数，流场每一时刻都在缓慢变化，看起来好像有非定常特征，但速度变化量一直降不下去，始终在某个水平抖动。

原因：两种情况。一种是收敛阈值设得太宽，比如1e-4就停，流场还没真正平衡；另一种是压力泊松迭代的tol太松，每步的残余散度没有被完全压掉，导致速度场被反复注入微小扰动。

解决：把定常判据改成连续两步速度场最大差小于1e-8，并让压力泊松迭代的tol也降到1e-8，再看曲线是否平滑收敛。如果残差在某个数值附近震荡不降，切割压力泊松迭代次数通常能直接见效。记住一个习惯：判断定常要看历史曲线，不要当场拍脑袋。

6. 进阶：从复现走向验证的快速检查法

拿到这份MATLAB示例后，我建议你先别急着改物理参数，而是花半天时间做三件验证，都很便宜。

第一件是衰减涡的时间收敛阶测试。固定网格在64乘64，从dt=0.01开始逐步减半，记录每一档的L2误差，看误差是否大约按dt的第一次方下降。如果是，说明时间推进和投影法配合正确；如果误差不随dt变化，说明时间误差不是主导，问题可能在边界或初始场。

第二件是顶盖驱动流的速度剖面验证。在定常状态下取出方腔水平中线上u速度随y的分布，和同雷诺数的公开基准数据对比。Re=100时中心处回流的负速度峰值大约在-0.2附近，如果你的结果差出10%以上，优先怀疑网格分辨率或角点平滑的影响。第三件是质量守恒检查。计算任意时刻全场散度的平方和，这个数应当在压力泊松迭代的容差量级，如果明显偏大，说明投影法第三步的梯度索引写错了。

性能优化方面，整套代码最大的瓶颈在压力泊松方程的Jacobi迭代，可以用MATLAB矩阵操作代替双重for循环，加速比通常在5倍以上。如果网格超过128乘128，把迭代改成预条件共轭梯度法会更划算，但要注意边界条件的组装方式不同。这份示例包自带的算例控制在几十秒内跑完，没必要为了性能牺牲代码可读性。

从那以后我每次拿到任何NS求解代码，都会强制自己先跑一遍衰减涡误差测试再谈别的。误差阶不对，后面的流场图再漂亮也是废的。希望这份MATLAB二维非定常示例能帮你少走这几步弯路，早日把注意力放到真正要做的流动问题上。

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