简介面向流体力学数值模拟学习者的二维非定常纳维斯托克斯方程完整示例非常适合正在学习计算流体力学或MATLAB编程的读者。资源以实现不可压缩流体在二维空间的非定常运动模拟为目标清晰演示了空间离散、时间推进、速度压强耦合、边界条件处理与压力修正等完整流程。压缩包内共有76个文件主体为44个MATLAB脚本与函数并配有29张流场演化结果图、1份说明文档及1份许可文本整个包体积仅463KB。已有1932人学习使用。通过分模块阅读代码可深入理解速度场更新、压力泊松方程求解及涡量演化等核心环节配合不同时刻的流场图像还能直观验证数值格式的正确性为后续开展流体工程仿真或算法改进提供扎实的实践基础。 先聊一个常见现象计算流体的课上了NS方程也推导过无数遍但真要自己动手在MATLAB里写一个二维非定常Navier-Stokes求解器大多数人卡在第一步——不知道从哪儿下手。二维不可压缩NS方程本身只有两三行可一旦落到离散网格上速度放哪、压力放哪、时间步取多大、边界怎么给、压力怎么解这些问题每一个都能劝退新手。这篇文章用最经典的方腔驱动流算例把一个能从MATLAB里直接跑起来的二维非定常NS求解链路完整拆开。适合两类人刚学完流体力学基础、想第一次动手写CFD代码的人以及想快速验证数值格式、需要稳定测试床的研究生和工程师。看完你会得到一个能出图、能对比文献数据的MATLAB算例也会理解每一步背后的取舍。1. 为什么拿二维非定常NS方程开刀方程与计算目标二维不可压缩Navier-Stokes方程的无量纲形式其实很紧凑动量方程加上连续性方程就是全部∂u/∂t u∂u/∂x v∂u/∂y -∂p/∂x (1/Re)(∂²u/∂x² ∂²u/∂y²) ∂v/∂t u∂v/∂x v∂v/∂y -∂p/∂y (1/Re)(∂²v/∂x² ∂²v/∂y²) ∂u/∂x ∂v/∂y 0无量纲化之后整个问题只剩下一个参数雷诺数Re。它的物理意义是惯性力与粘性力的比值Re低的时候流动被扩散主导Re高的时候对流主导流场结构会越来越丰富。这也是我推荐用它练手的原因——一个参数就能覆盖从层流到复杂涡结构的一大片物理现象。1.1 为什么是二维而不是直接上三维三维NS的问题在于变量直接变成u、v、w、p四个计算量按网格点数的三次方增长个人电脑跑一个像样的算例经常要等一小时起步。二维保留了NS的核心困难——对流项的非线性、不可压缩约束的处理、压力速度耦合——又把计算规模控制在笔记本能接受的范围。而且二维最大的红利是可以直接画流线和涡量云图流场结构一眼就能看懂这对验证代码正确性太重要了。1.2 为什么是非定常很多初学者以为定常问题更简单实际上定常NS的数值处理反而更绕。非定常把时间当作一个明确的推进方向每个时间步做一次预测-修正物理图像非常清楚。定常问题通常也要用伪时间推进来迭代本质还是在解非定常问题。从非定常入手还有一个额外的好处能看到启动阶段的涡生成和演化比一上来就追稳态解有趣得多。所以这个算例的目标就很明确了给定初始速度场和边界条件每一时间步求出一组满足连续性方程的速度场和对应的压力场。方腔驱动流跑到足够久之后会收敛到文献里的稳态解那就到了验收时刻。2. 先定方案投影法而不是涡量-流函数法二维不可压缩NS的主流数值方案有两条路线一个是本文采用的投影法也叫分步法另一个是涡量-流函数法。很多学校课程里先教涡量-流函数法因为它把变量从u、v、p三个消成ω、ψ两个连续性约束被流函数自动满足方程数量少数值上也稳。但这个方法有个硬伤边界涡量需要额外推导而且基本没办法推广到三维。再说投影法。它的核心思路是把速度更新拆成三步先忽略压力梯度只用对流项和扩散项算一个中间速度u*解压力泊松方程找一个能让中间速度散度归零的压力场用压力梯度修正速度得到满足连续性方程的速度场。数学上写出来就是u* u^n dt * ( -(u^n·∇)u^n (1/Re)∇²u^n ) ∇²p^{n1} (1/dt) * ∇·u* u^{n1} u* - dt * ∇p^{n1}第二步的方程习惯上直接写为∇²p (∇·u*)/dt来源是让修正后的速度散度等于零。这是整个算例里最需要理解清楚的一步后面单独展开。两种方案放在一起对比选型逻辑就很清楚了对比维度投影法速度-压力涡量-流函数法求解变量u, v, pω, ψ连续性约束投影修正强制满足流函数自动满足压力求解需要解泊松方程不需要显式压力三维推广直接扩展基本不可用边界条件速度边界直观涡量边界需推导我实际用下来的体会是投影法虽然每一步多一个压力泊松方程要解但整个求解逻辑与现代CFD的主流框架完全一致后面想往LES、DNS、或者三维方向走现在打的地基不会浪费。3. 离散落地的关键选择网格布置、差分格式与边界条件方程选定了下一步就是把连续方程搬到离散网格上。这里每一步选择都会影响后面的稳定性和收敛性我按顺序拆开讲。3.1 为什么要提交错网格不可压缩CFD里有一个著名的坑叫棋盘压力振荡如果速度、压力都放在同一个网格节点上离散后的压力修正方程会出现相邻节点解耦压力场看起来像棋盘一样一高一低交替实际却是错的。交错网格是标准解法——u放在x方向的半网格点上v放在y方向的半网格点上压力放在网格中心。这样每个压力点周围都有一圈真实的速度分量在撑场信息不会各玩各的。用生活里的类比这就像两把错开半个齿距的梳子速度插在压力的齿缝里谁也躲不开谁。本文为了代码可读性采用非交错网格的教学版本低雷诺数下跑方腔没有大问题压力场会有一点轻微棋盘但不影响速度场主结构。真要上高Re或者复杂边界请务必换成交错网格或者引入Rhie-Chow动量插值。3.2 空间差分和时间推进的选择对流项和扩散项我都用二阶中心差分。以x方向导数为例∂u/∂x ≈ (u(i1,j) - u(i-1,j)) / (2*dx)扩散项的拉普拉斯算子就是标准的五点格式。中心差分的优点是二阶精度、实现简单代价是在Re较高时对流项容易产生数值振荡。如果之后要跑Re2000的方腔可以考虑换成迎风格式或者QUICK格式网格也要相应加密。时间推进这里先用一阶显式Euler稳定区间窄但写起来最直观方便排查问题。代码跑通之后再升级到Adams-Bashforth二阶或者三阶Runge-Kutta都不难核心结构完全不用动。3.3 方腔驱动流的边界条件方腔驱动流是最经典的验证算例一个单位正方形腔体顶盖以水平速度向右运动其余三个壁面静止。具体来说就是顶盖u 1, v 0左、右、底壁u 0, v 0压力边界法向梯度∂p/∂n 0这里有一个容易踩的坑顶盖速度从0突然跳到1和两侧静止壁面之间形成了数学上的间断启动阶段速度场会被这个间断搅得天翻地覆前几步就可能让数值解直接爆掉。我的做法是对顶盖速度做一个线性升速utop min(n * dt / 2, 1); % 前2秒从0线性升到1这个小技巧成本极低但对稳定性的改善非常明显强烈建议保留。4. 压力泊松方程整个示例里最微妙的一环投影法能不能跑稳九成取决于压力泊松方程这一步。它本身不复杂但细节非常多。4.1 方程怎么来的回顾投影法的修正步u^{n1} u* - dt * ∇p^{n1}对两边取散度∇·u^{n1} ∇·u* - dt * ∇²p^{n1}我们希望下一时刻的速度场满足连续性方程也就是∇·u^{n1} 0于是∇²p^{n1} (1/dt) * ∇·u*这个方程里右端的散度完全由中间速度u*决定所以每步都要先算出u*、再求散度、再解这个泊松方程最后才能做速度修正。4.2 五点离散与SOR迭代在均匀网格上二维泊松方程用五点格式离散(p(i1,j) - 2p(i,j) p(i-1,j))/dx² (p(i,j1) - 2p(i,j) p(i,j-1))/dy² b(i,j)直接解这个线性方程组可以用MATLAB的反斜杠但要提前组装一个(Nx-1)×(Ny-1)维的大型稀疏矩阵代码会很长。教学算例里更常见的是SOR迭代公式是p(i,j) (1-ω)p(i,j) ω/(2/dx²2/dy²) * ( (p(i1,j)p(i-1,j))/dx² (p(i,j1)p(i,j-1))/dy² - b(i,j) )松弛因子ω一般取1.5到1.7我习惯取1.5稳。迭代停止条件用最大残差小于1e-640×40网格通常几百步就能收敛。下面是一个可以直接复制进脚本的SOR求解函数function p sorPoisson(p0, b, dx, dy, maxIter, tol) % 教学简化版SOR求解压力泊松方程 % 边界点不参与迭代等效于零法向梯度的近似处理 p p0; omega 1.5; dx2 dx^2; dy2 dy^2; denom 2/dx2 2/dy2; [Nx, Ny] size(p); for k 1:maxIter pOld p; for j 2:Ny-1 for i 2:Nx-1 p(i,j) (1-omega)*pOld(i,j) omega/denom * ( ... (p(i1,j)p(i-1,j))/dx2 ... (p(i,j1)p(i,j-1))/dy2 - b(i,j) ); end end if max(abs(p(:) - pOld(:))) tol break; end end end4.3 压力边界的处理方式严格CFD做法是在边界上离散∂p/∂n 0工程代码里常用p(1,:)p(2,:)这类一阶外推。教学版本经常直接让边界压力不更新也就是代码里SOR只扫内点。这带来的误差在低Re方腔算例里对速度场影响很小因为压力的绝对值本身不重要真正起作用的是压力梯度。这个简化可以接受但心里要清楚它的局限。如果发现SOR迭代不收敛优先检查两件事第一右端项b是不是只在内部点有值、边界上是不是也被填入了非零散度第二迭代内循环里用的是更新后的p还是旧pSOR公式要求边用新值边覆盖写错这一行就会从SOR退化成Jacobi收敛速度立刻掉一个量级。5. 方腔驱动流验证与MATLAB代码骨架算例写完之后怎么判断对不对不能只盯着屏幕看流场漂不漂亮。方腔驱动流的好处是文献数据极其丰富最常用的是Ghia等人在1982年发表的结果。以主涡涡心位置为例Re主涡涡心x主涡涡心y1000.6170.7344000.5550.60610000.5310.563跑稳态后在速度云图里找速度模量最小的区域坐标误差在3%以内就可以认为代码正确。注意Re100时大约需要t20的物理时间才能到近似稳态别跑几百步就下结论。5.1 MATLAB主循环代码下面是完整的主循环我尽量保持代码简短可读。直接复制到脚本里配合上面的sorPoisson函数就能跑% 参数 Nx 40; Ny 40; % 网格数 Re 100; nu 1/Re; % 无量纲粘性系数 Lx 1; Ly 1; dx Lx/Nx; dy Ly/Ny; dt 0.002; T 10; nt round(T/dt); x linspace(0, Lx, Nx1); y linspace(0, Ly, Ny1); u zeros(Nx1, Ny1); v zeros(Nx1, Ny1); p zeros(Nx1, Ny1); % 时间推进 for n 1:nt % 顶盖升速避免启动震荡 utop min(n * dt / 2, 1); % 2秒内线性升到1 % 固定边界 u(:, end) utop; u(:, 1) 0; u(1, :) 0; u(end, :) 0; v(:, end) 0; v(:, 1) 0; v(1, :) 0; v(end, :) 0; uo u; vo v; % 动量预测显式Euler暂时不处理压力 u(2:end-1,2:end-1) uo(2:end-1,2:end-1) - dt * ( ... uo(2:end-1,2:end-1) .* (uo(3:end,2:end-1)-uo(1:end-2,2:end-1))/(2*dx) ... vo(2:end-1,2:end-1) .* (uo(2:end-1,3:end)-uo(2:end-1,1:end-2))/(2*dy) ) ... dt*nu * ( (uo(3:end,2:end-1)-2*uo(2:end-1,2:end-1)uo(1:end-2,2:end-1))/dx^2 ... (uo(2:end-1,3:end)-2*uo(2:end-1,2:end-1)uo(2:end-1,1:end-2))/dy^2 ); v(2:end-1,2:end-1) vo(2:end-1,2:end-1) - dt * ( ... uo(2:end-1,2:end-1) .* (vo(3:end,2:end-1)-vo(1:end-2,2:end-1))/(2*dx) ... vo(2:end-1,2:end-1) .* (vo(2:end-1,3:end)-vo(2:end-1,1:end-2))/(2*dy) ) ... dt*nu * ( (vo(3:end,2:end-1)-2*vo(2:end-1,2:end-1)vo(1:end-2,2:end-1))/dx^2 ... (vo(2:end-1,3:end)-2*vo(2:end-1,2:end-1)vo(2:end-1,1:end-2))/dy^2 ); % 固定中间速度边界 u(:, end) utop; u(:, 1) 0; u(1, :) 0; u(end, :) 0; v(:, end) 0; v(:, 1) 0; v(1, :) 0; v(end, :) 0; % 压力泊松方程右端项 b zeros(Nx1, Ny1); b(2:end-1,2:end-1) ... ( (u(3:end,2:end-1) - u(1:end-2,2:end-1))/(2*dx) ... (v(2:end-1,3:end) - v(2:end-1,1:end-2))/(2*dy) ) / dt; p sorPoisson(p, b, dx, dy, 500, 1e-6); % 投影修正 u(2:end-1,2:end-1) u(2:end-1,2:end-1) - ... dt * (p(3:end,2:end-1) - p(1:end-2,2:end-1))/(2*dx); v(2:end-1,2:end-1) v(2:end-1,2:end-1) - ... dt * (p(2:end-1,3:end) - p(2:end-1,1:end-2))/(2*dy); % 修正后再固定一次边界 u(:, end) utop; u(:, 1) 0; u(1, :) 0; u(end, :) 0; v(:, end) 0; v(:, 1) 0; v(1, :) 0; v(end, :) 0; % 每1000步画一次流场 if mod(n, 1000) 0 [X, Y] meshgrid(x, y); subplot(1,2,1); contourf(X, Y, sqrt(u.^2 v.^2), 30); colorbar; axis equal; title(sprintf(t%.2f, Re%d, n*dt, Re)); subplot(1,2,2); quiver(X(1:2:end,1:2:end), Y(1:2:end,1:2:end), ... u(1:2:end,1:2:end), v(1:2:end,1:2:end)); axis equal; drawnow; end end代码的核心逻辑就是上面说的三步先算不含压力的中间速度再解压力泊松方程最后用压力梯度修正速度。中间两次重设边界条件是为了防止边界被内部计算污染这个习惯建议保留。5.2 快速验收流程拿到代码后不建议直接跑长时程。我的习惯是先跑500步看流场形状确认顶盖下方出现一个大涡、没有NaN然后再跑完整时程。Re100时40×40网格、dt0.002在我的笔记本上跑完整T10大概需要几分钟可以接受。网格加到64×64之后时间步要相应缩小总运行时间会明显上涨。6. 稳定性边界与调试实战跑不动的N种原因6.1 时间步的稳定性上限显式方法最烦人的就是稳定性限制。这个算例里有两道锁一道来自对流项另一道来自扩散项。对流CFL条件要求CFL (|u|max |v|max) * dt / min(dx, dy) 1方腔里速度最大值不超过顶盖速度1所以40×40网格下dxdx0.025dt要小于0.025才能过这一关。扩散项的限制更严格nu * dt / dx² 0.5Re100时nu0.01dx²0.000625算出dt0.03125。两条合起来dt取0.002是非常保守的适合教学想提性能可以试着放大到0.005但要随时盯紧速度场是否出现振荡。网格加密时务必记住空间步长减半显式扩散限制会让时间步缩到原来的四分之一。这就是为什么高Re、细网格下显式方法跑起来让人抓狂想突破就只能上隐式或半隐式。6.2 常见问题排查清单我把自己踩过的坑和身边人常遇到的状况整理成了一份排查顺序遇到问题按照这个顺序查效率最高输出NaN先查CFL条件再看边界条件在投影修正后有没有被重新更新最后检查压力SOR迭代是不是真的收敛了。三者都没毛病但还发散把dt直接除以2再试。压力场出现明显棋盘格教学版非交错网格的常见病。低Re下忍一忍没问题想彻底解决就换交错网格。涡心位置明显偏下或者偏上优先怀疑没跑到稳态方腔Re100需要t20其次怀疑网格太粗。启动阶段速度场剧烈震荡顶盖速度从0阶跃到1导致的用前面给的升速处理就能缓解。SOR迭代长期不收敛打印每次迭代的最大残差看看趋势如果残差在0.01附近抖动下不去检查右端项b的边界是不是也被赋了非零值。6.3 一个小建议跑通之后别着急收工建议把整个主循环包成一个函数输入参数只有Re、Nx、Ny和总时长输出稳态流场。这样后续做Re100、400、1000的参数扫描会非常方便也是我后来做各种数值实验的基础工具。我第一次跑通这个算例的时候前几版都因为dt取太大直接NaN后来养成了习惯任何新网格、新边界条件动手之前先按上面的两个稳定性公式算一遍上限再取一半作为初始dt。这个习惯帮我避开了后面很多无意义的debug时间。二维非定常NS的MATLAB示例其实不难难的是每一步都知道自己为什么这么写把这套流程走通一遍后面的三维、湍流、复杂边界都是在同一个骨架上做加法而已。本文还有配套的精品资源点击获取