简介：这份资源是面向流体力学数值模拟学习者与MATLAB使用者的二维非定常Navier-Stokes方程求解示例，聚焦不可压缩流体在二维空间内的运动描述，适合具备一定数值方法与编程基础、希望理解有限元或差分实现思路的中高级读者参考。压缩包共76个文件，以44个m脚本为核心，涵盖单元矩阵组装、形函数与高斯点计算、梯度与对流项处理、约束施加及初始化等模块，另有29个png结果图、1个html说明页和1个txt许可文件，整体约463KB，结构清晰便于按功能检索。已有1947人学习下载。读者可从中获得从网格定义、边界条件设置到时间推进与压力修正的完整代码框架，并借助速度矢量图、流线图等可视化结果对照理解方程各项的物理含义，为空气动力学、水力学等方向的数值模拟实践提供可复用的脚本基础与排错参考。

1. 从一张“跑不起来”的二维非定常流场图说起

很多人第一次在 MATLAB 里写二维非定常 Navier-Stokes 求解器，卡住的地方不是公式，而是时间推进那一步：压力泊松方程解出来的场和速度场对不上，跑几十步就发散，或者涡量莫名其妙地衰减。这份示例资源就是围绕“二维、非定常、不可压”这三个限定词搭起来的一套可运行骨架，把投影法（Projection Method）从数学推导落到 MATLAB 的矩阵运算里，包含网格生成、对流项与扩散项离散、压力修正、时间步进和流场可视化几个部分。它适合两类人：一类是刚接触 CFD 数值实现、想看懂“方程怎么变成代码”的学生；另一类是要在 MATLAB 环境里快速搭一个二维流场原型、后续再往里塞湍流模型或边界条件的研究者。整套代码不依赖工具箱之外的重量级库，拿到手改几个参数就能观察涡街、剪切层这类非定常现象。

2. 项目骨架与技术选型：为什么是 MATLAB 加投影法

2.1 目录结构与文件职责

这类二维非定常求解器在工程上通常拆成“驱动脚本 + 若干函数文件”的形式，主脚本负责参数、时间循环和出图，函数文件各管一块物理或数值逻辑。一个合格的目录组织大致如下，读者拿到源码后可以按这个映射去对号入座：

ns2d_unsteady/ ├── main.m % 驱动脚本：参数、初始化、时间推进、可视化 ├── init_fields.m % 流场初始化：速度、压力、涡量 ├── build_grid.m % 交错网格/同位网格坐标与算子 ├── assemble_convection.m % 对流项离散（迎风或中心差分） ├── assemble_diffusion.m % 扩散项离散（拉普拉斯算子） ├── solve_pressure.m % 压力泊松方程求解 ├── apply_bc.m % 边界条件施加 ├── compute_vorticity.m % 涡量后处理 └── plot_fields.m % 流场、涡量、压力云图绘制

主脚本main.m是唯一需要手动运行的入口，其余都是被调用的函数。这种拆法的好处是：调对流格式只动assemble_convection.m，换压力求解器只动solve_pressure.m，不会牵一发动全身。我一般会建议读者先把main.m从头到尾读一遍，把时间循环里的调用顺序画出来，再逐个进函数看实现，这样比一上来啃离散公式效率高得多。

2.2 为什么选投影法而不是直接解耦合方程

不可压 Navier-Stokes 的难点在于压力没有独立的演化方程，它通过连续性约束隐式地影响速度。直接离散耦合方程组会得到一个大型鞍点系统，MATLAB 里解起来又慢又难调。投影法把一步拆成两步：先忽略压力梯度算一个中间速度，再用压力泊松方程把中间速度投影到无散空间。这样每步只需要解一个标量泊松方程，矩阵规模小、条件数可控，特别适合 MATLAB 这种以稠密/稀疏矩阵运算见长的环境。

常见做法是采用 Chorin 型一阶投影，时间离散写成：

% 一步时间推进的骨架（示意） [u_star, v_star] = assemble_convection(u, v, dx, dy, dt); % 对流项显式 [u_star, v_star] = assemble_diffusion(u_star, v_star, nu, dx, dy, dt); % 扩散项 rhs = (u_star(2:end-1,2:end-1) - u_star(1:end-2,2:end-1))/dx ... + (v_star(2:end-1,2:end-1) - v_star(2:end-1,1:end-2))/dy; % 速度散度 p = solve_pressure(rhs, dx, dy); % 解泊松方程 [u, v] = correct_velocity(u_star, v_star, p, dx, dy, dt); % 压力修正

这段代码里dt受 CFL 条件约束，nu是运动黏度，dx、dy是网格间距。对流项显式、扩散项可以显式也可以隐式，示例里为了代码清晰通常都走显式，代价是时间步不能太大。读者如果发现跑着跑着数值爆炸，第一反应应该是把dt减半再看，而不是怀疑公式写错。

2.3 网格与变量存放方式

二维结构化网格上，速度分量和压力的存放位置决定了离散格式的写法。同位网格（collocated）实现简单，但容易出现棋盘格压力振荡；交错网格（staggered）把 u 存在竖直面上、v 存在水平面上、p 存在格心，天然抑制振荡，是教学示例里更常见的选择。示例代码如果采用交错网格，build_grid.m里会分别生成三套坐标，读者改边界条件时要特别注意 u、v 的索引范围比 p 各少一格，这是最容易翻车的地方。

提示：拿到源码先确认是同位还是交错网格，这决定了后面所有差分公式的索引偏移，搞错了会出现“压力场看着正常但速度场对不上”的玄学现象。

3. 核心模块拆解：从对流项到压力泊松方程

3.1 对流项离散：迎风还是中心差分

对流项是非定常流场里产生非线性行为的根源，也是数值不稳定的主要来源。中心差分精度高但在高雷诺数下容易振荡，一阶迎风稳定但数值耗散大、会抹平涡结构。示例里通常提供一个可切换的开关，读者可以两种都跑一遍对比涡量场。

function [cu, cv] = assemble_convection(u, v, dx, dy, scheme) % 输入：u,v 为速度场，dx,dy 网格间距，scheme 取 'upwind' 或 'central' % 输出：cu,cv 为对流项贡献 [nx, ny] = size(u); cu = zeros(nx, ny); cv = zeros(nx, ny); for i = 2:nx-1 for j = 2:ny-1 if strcmp(scheme, 'upwind') ux = (u(i,j) > 0) * (u(i,j)-u(i-1,j))/dx ... + (u(i,j) < 0) * (u(i+1,j)-u(i,j))/dx; else ux = (u(i+1,j) - u(i-1,j)) / (2*dx); end cu(i,j) = -u(i,j)*ux; % 只展示 u 方向分量 end end end

scheme参数控制格式选择，u(i,j) > 0这种写法是迎风方向的判断，正速度取后向差分、负速度取前向差分。实际工程代码里会用向量化代替双重循环，但教学示例保留循环是为了让索引关系一目了然。读者二次开发时如果换成向量化，务必先用小网格和循环版对拍，确认结果一致再替换，否则一个转置错误就能让你调一整天。

3.2 扩散项与黏性系数

扩散项就是拉普拉斯算子乘黏度，离散上比较直接，五点在每个内点做二阶中心差分即可。真正需要留意的是黏性系数nu和雷诺数的关系：示例里通常给定Re，再由特征长度和特征速度反算nu = U*L/Re。改雷诺数时不要只改一个数就完事，要确认特征尺度是否同步调整，否则你观察到的“流态变化”可能只是无量纲化没对齐。

function [du, dv] = assemble_diffusion(u, v, nu, dx, dy) % 五点差分拉普拉斯，nu 为运动黏度 du = nu * ( (u(3:end,2:end-1) - 2*u(2:end-1,2:end-1) + u(1:end-2,2:end-1))/dx^2 ... + (u(2:end-1,3:end) - 2*u(2:end-1,2:end-1) + u(2:end-1,1:end-2))/dy^2 ); dv = nu * ( (v(3:end,2:end-1) - 2*v(2:end-1,2:end-1) + v(1:end-2,2:end-1))/dx^2 ... + (v(2:end-1,3:end) - 2*v(2:end-1,2:end-1) + v(2:end-1,1:end-2))/dy^2 ); end

这里用切片一次性算完所有内点，返回的du、dv尺寸比原场小一圈，调用方要负责把边界值补回去。nu越大扩散越强、涡衰减越快，跑涡街算例时如果涡很快就消失，先检查nu是不是设大了。

3.3 压力泊松方程：求解器选择与边界处理

压力泊松方程是每步计算量最大的部分，形式为 ∇²p = ∇·u*/dt。离散后得到一个稀疏线性系统，示例里一般用mldivide（反斜杠）直接解，小网格下够用；网格一大就得换迭代法或预条件。

function p = solve_pressure(div_u, dx, dy, bc) % div_u 为中间速度散度，bc 指定压力边界类型 [nx, ny] = size(div_u); N = nx * ny; % 组装五点拉普拉斯稀疏矩阵 e = ones(N,1); Ax = spdiags([e -2*e e], -1:1, N, N) / dx^2; Ay = spdiags([e -2*e e], -ny:ny, N, N) / dy^2; A = Ax + Ay; b = div_u(:); % 诺伊曼压力边界需固定一个参考点消除零空间 A(1,:) = 0; A(1,1) = 1; b(1) = 0; p = reshape(A \ b, nx, ny); end

A(1,:) = 0; A(1,1) = 1这两行是处理纯诺伊曼边界的关键：压力只定到常数，不固定参考点矩阵就奇异，反斜杠会给出警告或垃圾解。这是新手最常踩的坑之一，现象是压力场整体漂移、速度场却看着还行，很容易被忽略。bc参数在完整实现里会区分固壁、入口、出口，示例里如果只做了固壁，读者扩展到开口边界时要自己补对应行。

3.4 时间推进与 CFL 监控

主循环里每一步都要检查 CFL 数，超过阈值就报警或自动缩步。示例代码通常把dt设成固定值，读者最好自己加一段监控：

cfl = max(max(abs(u)))*dt/dx + max(max(abs(v)))*dt/dy; if cfl > 0.5 warning('CFL=%.2f 偏大，考虑减小 dt', cfl); end

显式对流加显式扩散的稳定性条件更严，扩散项还额外要求nu*dt/dx^2 < 0.5。两个条件取交集才是安全步长。跑长时程算例前先做几步试算、打印 CFL 和扩散数，比跑崩了再回头找原因省事得多。

4. 跑起来：参数配置、运行流程与结果验证

4.1 环境准备与首次运行

MATLAB 版本上没有特殊要求，R2016b 之后的版本对隐式扩展支持较好，稀疏矩阵接口也稳定。不需要额外工具箱，spdiags、sparse、\都是基础功能。首次运行按这个顺序走：

1. 把整个目录加入 MATLAB 路径（右键文件夹 → 添加到路径 → 所选文件夹及子文件夹）。
2. 打开main.m，确认顶部参数区的网格数nx、ny，雷诺数Re，总步数nsteps，时间步dt。
3. 先跑一个小算例：nx=64; ny=64; nsteps=200;，观察命令行有没有奇异矩阵警告。
4. 运行结束后看plot_fields.m弹出的涡量云图，确认涡结构在演化而不是一片噪声。

第一次跑建议把网格压到 32×32、步数压到 50，几十秒内出结果，先确认整条链路通，再往上加规模。直接上 256×256 跑几千步，中途出错你都不知道是哪一步崩的。

4.2 关键参数怎么调

改Re时同步检查nu的换算，改nx时同步检查dt是否还满足稳定性条件，这两个联动是最容易被忽略的。

4.3 结果对不对：三个自检手段

数值解没有解析解对照时，用这几个办法判断实现是否正确。第一，看散度：每步修正后计算max(abs(div(u)))，应该在机器精度量级，如果一直在 1e-2 以上说明压力修正没生效。第二，看动能随时间的变化：无强迫的衰减湍流里动能应单调下降，如果出现非物理增长，多半是对流项或时间步有问题。第三，做网格收敛：同一算例 64×64 和 128×128 的涡量场大尺度结构应相似，差异只在小尺度，如果大结构都变了，说明离散或边界条件有系统性错误。

注意：投影法一阶格式在时间上精度有限，做定量对比时不要拿它和二阶格式的结果逐点比，看趋势和统计量更靠谱。

5. 避坑与排查：那些让流场“炸掉”的细节

5.1 现象：跑几十步后速度场出现高频振荡

原因通常是对流项用了中心差分而网格又不够密，局部 Peclet 数超过 2，中心差分失去稳定性。解决办法是切回迎风格式先跑稳，或者加密网格、降低雷诺数，让网格 Peclet 数回到安全区。不要试图靠减小dt解决，时间步再小也压不住空间离散的振荡。

5.2 现象：压力场整体漂移、每步加一个常数

这是纯诺伊曼压力边界没固定参考点的典型表现，矩阵奇异，反斜杠给出的解带任意常数。解决就是在组装矩阵后把第一个未知量钉死，如 3.3 节代码所示。如果用的是迭代求解器，还要确保右端项和矩阵都做了同样的处理，否则迭代不收敛。

5.3 现象：涡量场看起来正常，但速度散度不为零

多半是压力修正那一步的索引错位，交错网格下 u、v、p 的尺寸不同，修正时切片范围写错一格就会让散度残留。排查办法是把中间速度、压力、修正后速度三个场的尺寸打印出来，逐个核对，再在单个内点上手算一遍修正公式和代码结果对比。

5.4 现象：换了大网格后内存爆掉或求解极慢

反斜杠解稀疏系统在二维几万未知量时还能接受，再大就该换迭代法（如pcg配不完全 Cholesky 预条件）或直接用多重网格。示例代码如果只写了直接解，读者扩规模前先评估内存，A的稀疏存储虽然省空间，但分解产生的填充可能远超预期。

5.5 现象：边界条件改了之后入口附近出现非物理回流

固壁无滑移和入口给定速度的施加顺序有讲究，先设入口再设固壁，或者反过来，结果可能不同。正确做法是先给所有边界赋默认值，再按边界类型逐个覆盖，最后统一处理角点。角点归属不明确时，回流往往就从那里冒出来。

6. 二次开发：把示例改成你自己的算例

这套代码最大的价值是骨架清晰，改造时按“物理 → 离散 → 求解 → 后处理”四层动刀，不要一上来就重写主循环。我一般会先做一件事：把main.m里的参数区抽成一个config.m返回结构体，这样换算例只改配置文件，主脚本不动，版本管理也干净。

第一个常见改造是加力项。比如加一个体积力驱动或浮力项，只需要在投影之前给中间速度叠加force*dt，其余流程不变。第二个是换边界条件，把固壁换成周期边界，要在apply_bc.m里同时处理左右和上下两对边界，并且注意压力泊松方程在周期边界下的可解性条件（右端项均值为零），否则又回到奇异矩阵的老问题。第三个是提升时间精度，把一阶投影换成二阶的 Adams-Bashforth 对流加 Crank-Nicolson 扩散，改动集中在时间循环和对流、扩散两个函数的调用方式上，压力步不变。

验证改造是否正确的习惯做法是：每改一处，先跑一个已知有解析解或标准基准的退化算例（比如无对流纯扩散、或低雷诺数定常解），确认新代码在退化情形下和旧结果一致，再上完整非定常算例。这个“先退化再完整”的流程帮我省过很多次通宵排查。

% config.m 示例：把参数集中管理 function cfg = config() cfg.nx = 128; cfg.ny = 128; cfg.Re = 200; cfg.U = 1.0; cfg.L = 1.0; cfg.nu = cfg.U * cfg.L / cfg.Re; cfg.dt = 0.2 * min(1/(cfg.U/cfg.nx), cfg.nu); % 保守估计 cfg.nsteps = 2000; cfg.scheme = 'upwind'; end

dt这里用了一个保守估计式，实际跑起来再用 3.4 节的 CFL 监控微调。参数集中之后，做参数扫描就是外层套一个循环改cfg.Re，结果批量存盘，比手动改脚本可靠得多。

从那以后我每次拿到这类示例代码，都强制先跑最小算例、打印散度和 CFL、再动任何一行物理相关的代码，这个习惯让我少走了很多“改完不知道哪坏了”的弯路。希望帮到你。

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