企业官网建设流程全解析

简介：面向流体力学数值模拟学习者的二维非定常纳维斯托克斯方程完整示例，非常适合正在学习计算流体力学或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三个消成ω、ψ两个，连续性约束被流函数自动满足，方程数量少，数值上也稳。但这个方法有个硬伤：边界涡量需要额外推导，而且基本没办法推广到三维。

再说投影法。它的核心思路是把速度更新拆成三步：

1. 先忽略压力梯度，只用对流项和扩散项算一个中间速度u*；
2. 解压力泊松方程，找一个能让中间速度散度归零的压力场；
3. 用压力梯度修正速度，得到满足连续性方程的速度场。

数学上写出来就是：

u* = u^n + dt * ( -(u^n·∇)u^n + (1/Re)∇²u^n ) ∇²p^{n+1} = (1/dt) * ∇·u* u^{n+1} = u* - dt * ∇p^{n+1}

第二步的方程习惯上直接写为∇²p = (∇·u*)/dt，来源是让修正后的速度散度等于零。这是整个算例里最需要理解清楚的一步，后面单独展开。

两种方案放在一起对比，选型逻辑就很清楚了：

我实际用下来的体会是，投影法虽然每一步多一个压力泊松方程要解，但整个求解逻辑与现代CFD的主流框架完全一致，后面想往LES、DNS、或者三维方向走，现在打的地基不会浪费。

3. 离散落地的关键选择：网格布置、差分格式与边界条件

方程选定了，下一步就是把连续方程搬到离散网格上。这里每一步选择都会影响后面的稳定性和收敛性，我按顺序拆开讲。

3.1 为什么要提交错网格

不可压缩CFD里有一个著名的坑叫棋盘压力振荡：如果速度、压力都放在同一个网格节点上，离散后的压力修正方程会出现"相邻节点解耦"，压力场看起来像棋盘一样一高一低交替，实际却是错的。

交错网格是标准解法——u放在x方向的半网格点上，v放在y方向的半网格点上，压力放在网格中心。这样每个压力点周围都有一圈真实的速度分量在"撑场"，信息不会各玩各的。用生活里的类比，这就像两把错开半个齿距的梳子，速度插在压力的齿缝里，谁也躲不开谁。

本文为了代码可读性采用非交错网格的教学版本，低雷诺数下跑方腔没有大问题，压力场会有一点轻微棋盘但不影响速度场主结构。真要上高Re或者复杂边界，请务必换成交错网格，或者引入Rhie-Chow动量插值。

3.2 空间差分和时间推进的选择

对流项和扩散项我都用二阶中心差分。以x方向导数为例：

∂u/∂x ≈ (u(i+1,j) - u(i-1,j)) / (2*dx)

扩散项的拉普拉斯算子就是标准的五点格式。中心差分的优点是二阶精度、实现简单，代价是在Re较高时对流项容易产生数值振荡。如果之后要跑Re>2000的方腔，可以考虑换成迎风格式或者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^{n+1} = u* - dt * ∇p^{n+1}

对两边取散度：

∇·u^{n+1} = ∇·u* - dt * ∇²p^{n+1}

我们希望下一时刻的速度场满足连续性方程，也就是∇·u^{n+1} = 0，于是：

∇²p^{n+1} = (1/dt) * ∇·u*

这个方程里，右端的散度完全由中间速度u*决定，所以每步都要先算出u*、再求散度、再解这个泊松方程，最后才能做速度修正。

4.2 五点离散与SOR迭代

在均匀网格上，二维泊松方程用五点格式离散：

(p(i+1,j) - 2p(i,j) + p(i-1,j))/dx² + (p(i,j+1) - 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(i+1,j)+p(i-1,j))/dx² + (p(i,j+1)+p(i,j-1))/dy² - b(i,j) )

松弛因子ω一般取1.5到1.7，我习惯取1.5，稳。迭代停止条件用最大残差小于1e-6，40×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(i+1,j)+p(i-1,j))/dx2 + ... (p(i,j+1)+p(i,j-1))/dy2 - b(i,j) ); end end if max(abs(p(:) - pOld(:))) < tol break; end end end

4.3 压力边界的处理方式

严格CFD做法是在边界上离散∂p/∂n = 0，工程代码里常用p(1,:)=p(2,:)这类一阶外推。教学版本经常直接让边界压力不更新，也就是代码里SOR只扫内点。这带来的误差在低Re方腔算例里对速度场影响很小，因为压力的绝对值本身不重要，真正起作用的是压力梯度。这个简化可以接受，但心里要清楚它的局限。

如果发现SOR迭代不收敛，优先检查两件事：第一，右端项b是不是只在内部点有值、边界上是不是也被填入了非零散度；第二，迭代内循环里用的是更新后的p还是旧p，SOR公式要求边用新值边覆盖，写错这一行就会从SOR退化成Jacobi，收敛速度立刻掉一个量级。

5. 方腔驱动流验证与MATLAB代码骨架

算例写完之后怎么判断对不对？不能只盯着屏幕看流场漂不漂亮。方腔驱动流的好处是文献数据极其丰富，最常用的是Ghia等人在1982年发表的结果。以主涡涡心位置为例：

跑稳态后，在速度云图里找速度模量最小的区域，坐标误差在3%以内就可以认为代码正确。注意Re=100时大约需要t=20的物理时间才能到近似稳态，别跑几百步就下结论。

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, Nx+1); y = linspace(0, Ly, Ny+1); u = zeros(Nx+1, Ny+1); v = zeros(Nx+1, Ny+1); p = zeros(Nx+1, Ny+1); % ===== 时间推进 ===== 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(Nx+1, Ny+1); 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，然后再跑完整时程。Re=100时，40×40网格、dt=0.002，在我的笔记本上跑完整T=10大概需要几分钟，可以接受。网格加到64×64之后，时间步要相应缩小，总运行时间会明显上涨。

6. 稳定性边界与调试实战：跑不动的N种原因

6.1 时间步的稳定性上限

显式方法最烦人的就是稳定性限制。这个算例里有两道锁：一道来自对流项，另一道来自扩散项。

对流CFL条件要求：

CFL = (|u|max + |v|max) * dt / min(dx, dy) < 1

方腔里速度最大值不超过顶盖速度1，所以40×40网格下dx=dx=0.025，dt要小于0.025才能过这一关。

扩散项的限制更严格：

nu * dt / dx² < 0.5

Re=100时nu=0.01，dx²=0.000625，算出dt<0.03125。两条合起来，dt取0.002是非常保守的，适合教学；想提性能，可以试着放大到0.005，但要随时盯紧速度场是否出现振荡。

网格加密时务必记住：空间步长减半，显式扩散限制会让时间步缩到原来的四分之一。这就是为什么高Re、细网格下显式方法跑起来让人抓狂，想突破就只能上隐式或半隐式。

6.2 常见问题排查清单

我把自己踩过的坑和身边人常遇到的状况整理成了一份排查顺序，遇到问题按照这个顺序查，效率最高：

1. 输出NaN：先查CFL条件，再看边界条件在投影修正后有没有被重新更新，最后检查压力SOR迭代是不是真的收敛了。三者都没毛病但还发散，把dt直接除以2再试。
2. 压力场出现明显棋盘格：教学版非交错网格的常见病。低Re下忍一忍没问题，想彻底解决就换交错网格。
3. 涡心位置明显偏下或者偏上：优先怀疑没跑到稳态，方腔Re=100需要t>=20；其次怀疑网格太粗。
4. 启动阶段速度场剧烈震荡：顶盖速度从0阶跃到1导致的，用前面给的升速处理就能缓解。
5. SOR迭代长期不收敛：打印每次迭代的最大残差看看趋势，如果残差在0.01附近抖动下不去，检查右端项b的边界是不是也被赋了非零值。

6.3 一个小建议

跑通之后别着急收工，建议把整个主循环包成一个函数，输入参数只有Re、Nx、Ny和总时长，输出稳态流场。这样后续做Re=100、400、1000的参数扫描会非常方便，也是我后来做各种数值实验的基础工具。

我第一次跑通这个算例的时候，前几版都因为dt取太大直接NaN，后来养成了习惯：任何新网格、新边界条件动手之前，先按上面的两个稳定性公式算一遍上限，再取一半作为初始dt。这个习惯帮我避开了后面很多无意义的debug时间。二维非定常NS的MATLAB示例其实不难，难的是每一步都知道自己为什么这么写；把这套流程走通一遍，后面的三维、湍流、复杂边界，都是在同一个骨架上做加法而已。

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