简介：这份资源面向学习计算流体力学与数值方法的高校学生、科研人员和工程技术人员，提供二维非定常Navier-Stokes方程的MATLAB求解示例，帮助读者理解不可压缩流体在二维空间内运动的数值模拟流程。压缩包共76个文件，约463KB，以44个m脚本为核心，涵盖网格定义、形函数与高斯点计算、对流项与粘性项矩阵组装、压力修正及边界约束等模块，另有29张png结果图展示速度场与流线演化，并附1个html说明页和1个txt许可文件，便于按模块对照代码与计算结果。已有1947人学习下载，适合作为有限元或有限差分求解N-S方程的入门实践参考。读者可从中获得从方程离散、矩阵组装到时间推进与可视化输出的完整实现思路，并借助示例图直观理解非定常流动的演化过程，为空气动力学、水力学等方向的数值模拟打下基础。

1. 二维非定常 Navier-Stokes 的 MATLAB 示例：一份能跑起来的 CFD 入门资源

做流体仿真的人大概都有过这种体验：翻遍教材，N-S 方程的推导写得清清楚楚，可一到自己动手写代码，压力泊松方程怎么解、对流项用迎风还是中心差分、时间推进用显式还是隐式，全成了玄学。这份二维非定常 Navier-Stokes 的 MATLAB 示例，解决的正是这个断层——它把不可压 N-S 方程的数值求解流程拆成可执行的脚本，从网格生成、离散、时间推进到涡量场可视化，一条链路走通。适合刚接触 CFD 数值实现的研究生、需要快速验证算法思路的工程师，以及想把教材公式变成能跑代码的自学者。它不追求工业级求解器的鲁棒性，但胜在结构清晰、参数可调、每一步都能对应到方程里的具体项，拿来当数值方法的练手底稿正合适。

2. 方程怎么离散：从连续 N-S 到可迭代的矩阵形式

2.1 不可压 N-S 的控制方程与求解策略选型

二维不可压非定常 N-S 方程由动量方程和连续性方程组成。动量方程描述速度场的对流、扩散和压力梯度作用，连续性方程约束速度场散度为零。数值求解的核心难点在于压力与速度的耦合——压力没有独立的演化方程，它作为拉格朗日乘子存在，必须通过连续性约束间接求出。

常见的处理策略有三类：原始变量法中的投影法（分步法）、涡量-流函数法、以及基于交错网格的 SIMPLE 系列算法。这份示例走的是投影法路线，原因很实际：投影法把每个时间步拆成「先算中间速度、再解压力泊松方程修正速度」两步，逻辑线性、代码结构清楚，不需要组装大型稀疏矩阵求解耦合系统，对 MATLAB 这种解释型环境友好。涡量-流函数法虽然天然满足连续性，但边界条件处理绕，涡量边界条件在固壁处需要额外构造；SIMPLE 系列更适合稳态或大密度比问题，非定常推进时每步都要迭代到收敛，代码量和调试成本都上去了。

投影法的代价是压力泊松方程每步都要解一次。二维问题下，如果网格是结构化的，可以用快速泊松求解器（基于 DST/DCT）把复杂度压到接近 O(N log N)；示例里为了代码可读性，通常用稀疏矩阵直接求解或简单的迭代格式，网格规模控制在 128×128 以内时完全够用。

2.2 空间离散：交错网格与差分格式的取舍

空间离散采用交错网格（staggered grid），压力存在单元中心，u 速度存在竖直面上，v 速度存在水平面上。这样做的直接好处是压力梯度和散度算子天然紧致，不会出现棋盘格压力振荡——如果用同位网格配中心差分，压力和速度在同一套节点上，连续性方程对压力的约束会出现奇偶失联，压力场会长出高频振荡，这是新手最容易翻车的地方之一。

对流项默认用一阶迎风差分，扩散项用中心差分。迎风差分引入数值耗散，在低雷诺数下影响不大，但 Re 上到几千时涡结构会被抹平。示例里一般留了一个开关，可以切到二阶迎风或 QUICK 格式，代价是代码里要多写一层模板。我一般建议先用一阶迎风把流场跑稳、确认压力泊松方程收敛正常，再换高阶格式对比涡量场，否则一上来就调高阶格式，出了振荡根本分不清是格式问题还是边界条件写错了。

时间推进用显式 Euler 或二阶 Adams-Bashforth。显式格式受 CFL 条件约束，时间步长要满足 dt ≤ min(网格尺度/当地速度) 量级，粘性项还额外要求 dt ≤ Re·h²/4 左右。示例里通常给一个自适应步长的小函数，按当前最大速度自动缩 dt，避免手动调参调到崩溃。

2.3 压力泊松方程的组装与求解

压力泊松方程是投影法的核心。中间速度算出后，对其取散度，令修正后的速度散度为零，得到关于压力修正量的泊松方程。离散后每个单元中心一个方程，系数矩阵是五点拉普拉斯模板，用稀疏矩阵存储。

% 组装压力泊松方程的稀疏矩阵（五点差分，Neumann 边界） Nx = 64; Ny = 64; % 网格数 N = Nx * Ny; % 未知量个数 dx = 1.0 / Nx; dy = 1.0 / Ny; % 一维拉普拉斯模板， kron 拼成二维 e = ones(Nx, 1); L1d = spdiags([e -2*e e], -1:1, Nx, Nx) / dx^2; Ix = speye(Nx); Iy = speye(Ny); % 二维五点模板：x 方向 + y 方向 L = kron(Iy, L1d) + kron(spdiags([e -2*e e], -1:1, Ny, Ny) / dy^2, Ix); % 右端项：中间速度的散度 rhs = reshape(div_u_star, N, 1); % 直接求解（小网格够用），大网格换 pcg + 预条件 p = L \ rhs; p = reshape(p, Ny, Nx);

这段代码的关键点有三个。第一，kron拼装二维算子时要注意维度顺序，MATLAB 列优先存储，kron(Iy, L1d)对应的是沿 x 方向的二阶差分作用在按列展开的向量上，顺序搞反了压力场会整体错位。第二，纯 Neumann 边界下泊松方程的解不唯一（差一个常数），直接\求解可能报奇异或给出漂移的压力场，常见做法是固定一个参考点压力为零，或者在求解后减去全场均值。第三，网格超过 128×128 后直接求解开始吃力，换成pcg配不完全 Cholesky 预条件，收敛容差设 1e-6 量级即可，投影法对压力精度要求没那么苛刻。

2.4 时间推进主循环的骨架

主循环把上述模块串起来：算中间速度 → 解压力 → 修正速度 → 更新涡量 → 可视化。下面是一个精简骨架，省略了具体差分模板的展开。

% 非定常推进主循环 t = 0; tEnd = 2.0; dt = 1e-3; while t < tEnd % 1. 计算中间速度（对流 + 扩散，不含压力梯度） [u_star, v_star] = advect_diffuse(u, v, dx, dy, dt, Re); % 2. 组装并求解压力泊松方程 div_star = compute_divergence(u_star, v_star, dx, dy); p = solve_pressure(div_star, dx, dy); % 内部含稀疏求解 % 3. 用压力梯度修正速度 [u, v] = correct_velocity(u_star, v_star, p, dx, dy, dt); % 4. 更新涡量用于可视化 omega = compute_vorticity(u, v, dx, dy); % 5. 自适应步长（可选） dt = min(dt_max, 0.4 * min(dx, dy) / max(abs(u(:)) + 1e-8)); t = t + dt; if mod(round(t/dt), 50) == 0 visualize(omega, t); % 每 50 步画一次 end end

advect_diffuse里对流项用迎风、扩散项用中心差分，边界处固壁速度设零、出口用零梯度外推。correct_velocity里压力梯度用中心差分，注意交错网格上 u 和 v 的梯度位置不同，索引偏移要写对——这是调试时最容易出索引越界或错位的地方。自适应步长那行是保险丝，速度场突然发散时 dt 会被压小，给你留出观察窗口，不至于一步炸到 Inf。

3. 边界条件与初始场：让方腔驱动和圆柱绕流都能跑

3.1 顶盖驱动方腔的边界实现

顶盖驱动方腔是最经典的验证算例：四壁无滑移，顶盖以恒定速度 U=1 平移，Re 从 100 到 1000 对应不同的涡结构。边界条件在交错网格上的实现要分方向处理。

u 速度存在竖直面上，左右壁面处 u 直接设零；上下壁面处 u 的位置在壁面内侧半个网格，顶盖的 U 通过镜像虚拟节点施加，即 u(i, Ny+1) = 2U - u(i, Ny)。v 速度同理，上下壁面设零，左右壁面用镜像。压力边界用 Neumann 条件，法向梯度为零，离散时在边界单元上把泊松方程的模板截断即可。

% 顶盖驱动方腔边界条件（交错网格） U_lid = 1.0; % 左右壁：u = 0 u(1, :) = 0; u(end, :) = 0; % 上下壁：u 用镜像施加无滑移/驱动 u(:, 1) = -u(:, 2); % 底壁 u=0 u(:, end) = 2*U_lid - u(:, end-1); % 顶盖驱动 % 上下壁：v = 0 v(:, 1) = 0; v(:, end) = 0; % 左右壁：v 用镜像 v(1, :) = -v(2, :); v(end, :) = -v(end-1, :);

镜像公式里的系数 2 来自壁面到第一个内点的距离是半格，虚拟点距离壁面也是半格，线性插值后得到这个因子。写错成u(:,end) = U_lid会让顶盖附近速度剖面偏硬，涡心位置对不上基准解。Re=1000 时方腔会出现二级涡，网格至少 128×128 才能分辨，64×64 只能看到主涡，这是网格分辨率的硬约束，不是代码 bug。

3.2 圆柱绕流的来流与远场处理

圆柱绕流用来验证非定常特性：Re 超过约 47 后尾迹开始周期性脱落，升力系数做正弦振荡。几何上用浸没边界或贴体网格都行，示例里为了简单通常用阶梯近似——把圆柱占据的网格标记为固体，速度强制为零，压力在该处用相邻流体单元外推。

来流边界设在左侧，u=U∞、v=0；右侧出口用对流边界条件 ∂u/∂t + U∞∂u/∂x = 0，比零梯度更稳定，能减少出口反射。上下边界用滑移或远场条件。初始场给均匀来流加一个小扰动（比如在圆柱后方加一个高斯型速度脉冲），否则对称流态下数值舍入误差要很久才能激发出脱落，等得让人怀疑代码坏了。

% 圆柱区域标记与速度强制 [X, Y] = meshgrid(linspace(0, Lx, Nx), linspace(0, Ly, Ny)); solid = ((X - xc).^2 + (Y - yc).^2) <= r^2; % 圆柱内部 u(solid) = 0; v(solid) = 0; % 固体无滑移 % 出口对流边界（右边界） for j = 1:Ny u(end, j) = u(end-1, j) - U_inf * dt/dx * (u(end-1, j) - u(end-2, j)); end

出口对流条件里 dt/dx 乘的是流向梯度，系数超过 1 时会不稳定，所以 dt 要受这个条件额外约束。圆柱算例的 CFL 通常比方腔更紧，因为圆柱附近局部速度可以到 1.5U∞ 以上。阶梯近似在粗网格下会把圆柱「放大」成方块状，Strouhal 数偏高，想对基准数据得用浸没边界法或贴体网格，这是这份示例的边界所在，心里要有数。

4. 避坑与排查：五个让流场炸掉的常见问题

4.1 压力场出现棋盘格振荡

现象：涡量场看着还行，但压力云图出现高频明暗交替的格子状图案，相邻单元压力正负跳变。原因：用了同位网格配中心差分，压力和速度存在同一套节点上，连续性方程对压力的约束出现奇偶失联。解决：换交错网格，或者在同位网格上用 Rhie-Chow 动量插值构造面速度。示例默认交错网格，如果你改结构时把 u、v、p 放到了同一位置，这个坑必然出现。

4.2 时间步稍大就发散到 Inf

现象：跑几十步后速度指数增长，几秒内变 NaN。原因：显式格式违反 CFL 或粘性稳定性条件，尤其网格加密后 dx 变小，dt 上限按 dx² 缩，忘了同步调小。解决：先算dt_cfl = 0.4*min(dx,dy)/max(|u|)和dt_visc = 0.25*Re*min(dx,dy)^2，取两者较小值，再乘 0.8 安全系数。别迷信自适应步长函数，它只在速度已经偏大时才缩，发散初期可能来不及。

4.3 压力泊松方程求解报奇异或警告

现象：L \ rhs返回 NaN 或提示矩阵接近奇异。原因：纯 Neumann 边界下泊松方程解不唯一，系数矩阵秩亏一。解决：固定一个参考单元压力为零（把该行改成单位行、右端设零），或者求解后p = p - mean(p(:))。后者更省事，但要注意如果压力场本身有物理上的整体漂移（比如封闭方腔），减均值是合理的；有进出口的算例则要小心别把真实的压降抹掉。

4.4 涡量场对称性异常或涡心偏移

现象：方腔算例的涡心不在几何中心偏下的位置，或者左右不对称。原因：边界镜像公式系数写错（比如顶盖写成u(:,end)=U_lid而非2U-u(:,end-1)），或者交错网格索引偏移搞反，u 和 v 的梯度位置对不上。解决：先用 Re=100 的基准解对比涡心坐标，误差应在网格尺度量级；再检查correct_velocity里压力梯度的索引，u 在 x 面、v 在 y 面，梯度取法不同，写反了速度修正方向就错了。

4.5 长时间跑后质量不守恒、总动量漂移

现象：跑几千步后平均速度缓慢漂移，或者总动能不降反升。原因：离散散度算子与梯度算子不满足伴随关系（离散意义下的 ∇·∇p ≠ ∇²p），投影步残留散度；或者出口边界反射累积。解决：检查散度计算和梯度计算是否用了同一套面通量定义；出口改用对流边界；每隔若干步做一次散度清理（再投影一次）。这份示例规模下，前两条做到基本就稳了。

5. 进阶用法：用涡量场和基准数据验证你的实现

代码能跑不等于算对。验证分三层：网格收敛性、基准算例对比、守恒量监控。

网格收敛性最直接：同一算例跑 32×32、64×64、128×128 三档，看涡心位置、中心线速度剖面、Strouhal 数是否随网格加密收敛。方腔 Re=1000 的基准数据里，主涡心大约在 (0.53, 0.565) 附近（归一化坐标），二级涡在左下角和右下角。你的结果偏差超过 5% 就先查边界和分辨率，别急着调格式。

守恒量监控建议每步记录总动能和总涡量。无外力、无进出口的封闭方腔，总动能应单调衰减（粘性耗散），如果出现回升，说明数值格式在局部注入能量，多半是对流项中心差分在高速区失稳。圆柱绕流则监控升力系数的振荡幅值和频率，Strouhal 数 St = fD/U∞ 在 Re=100 时约 0.16～0.17，对不上先查圆柱直径是按网格数还是物理长度算的。

% 每步监控：动能、最大散度、涡量极值 KE = 0.5 * sum(u(:).^2 + v(:).^2) * dx * dy; divMax = max(abs(compute_divergence(u, v, dx, dy)), [], 'all'); omegaMax = max(abs(omega), [], 'all'); fprintf('t=%.3f KE=%.6e maxDiv=%.3e maxOmega=%.3f\n', t, KE, divMax, omegaMax);

maxDiv是投影法质量的直接指标，稳定运行时应在 1e-8 到 1e-6 量级（取决于求解容差），如果持续在 1e-3 以上，说明压力泊松方程没解到位或者散度/梯度算子不匹配。KE曲线拿来做收敛判断比看云图靠谱得多。

一个具体技巧：把涡量场存成矩阵序列，用VideoWriter导出 mp4，比逐帧截图省事，回放时涡的合并、脱落过程一目了然。我习惯在每次改完格式或边界后，先跑 Re=100 的短算例（2000 步以内）看maxDiv和KE曲线是否正常，再上高雷诺数长算例——这个习惯帮我省掉了无数次跑到一半发现边界写错、从头再来的时间。从那以后我每次动离散模板之前都强制走一遍这个短算例自检，希望帮到你。

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