简介这份资源是面向流体力学数值模拟学习者与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(i1,j)-u(i,j))/dx; else ux (u(i1,j) - u(i-1,j)) / (2*dx); end cu(i,j) -u(i,j)*ux; % 只展示 u 方向分量 end end endscheme参数控制格式选择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); endA(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、\都是基础功能。首次运行按这个顺序走把整个目录加入 MATLAB 路径右键文件夹 → 添加到路径 → 所选文件夹及子文件夹。打开main.m确认顶部参数区的网格数nx、ny雷诺数Re总步数nsteps时间步dt。先跑一个小算例nx64; ny64; nsteps200;观察命令行有没有奇异矩阵警告。运行结束后看plot_fields.m弹出的涡量云图确认涡结构在演化而不是一片噪声。第一次跑建议把网格压到 32×32、步数压到 50几十秒内出结果先确认整条链路通再往上加规模。直接上 256×256 跑几千步中途出错你都不知道是哪一步崩的。4.2 关键参数怎么调参数含义调整建议nx, ny网格点数从 64 起步翻倍前先确认 CFL 仍安全Re雷诺数涡街算例常见 1001000越高越需要细网格dt时间步长由 CFL 和扩散数反推不要凭感觉设nsteps总步数至少覆盖几个涡脱落周期才能看出非定常特征scheme对流格式先迎风跑稳再换中心差分看精度提升改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; enddt这里用了一个保守估计式实际跑起来再用 3.4 节的 CFL 监控微调。参数集中之后做参数扫描就是外层套一个循环改cfg.Re结果批量存盘比手动改脚本可靠得多。从那以后我每次拿到这类示例代码都强制先跑最小算例、打印散度和 CFL、再动任何一行物理相关的代码这个习惯让我少走了很多“改完不知道哪坏了”的弯路。希望帮到你。本文还有配套的精品资源点击获取