企业官网建设流程全解析

简介：这套MATLAB示例面向流体力学初学者与数值计算进阶者，围绕二维非定常Navier-Stokes方程，完整演示从网格定义、边界条件设置、速度场与压力场初始化，到时间推进、压力泊松方程修正及矢量图与流线可视化的求解链路。压缩包共76个文件，其中44个.m源码脚本覆盖组装刚度矩阵、对流项、黏性项、载荷向量、梯度算子与形状函数等核心函数，29张PNG图为不同时刻速度与压力演化结果，另含HTML报告与许可说明，整体仅463KB，目录结构清晰，便于快速定位。已有1932人学习下载。通过可运行代码与配套图形，读者能直观理解不可压缩流场的时间离散、投影方法与边界处理技巧；模块化的函数设计也方便直接移植到课程设计或科研仿真中，作为二次开发的起点。 二维非定常Navier-Stokes方程的MATLAB示例，算是计算流体力学入门里绕不开的一块硬骨头。我当年第一次尝试用MATLAB去解这个方程，对着公式看了半天觉得都懂，结果一运行就发散到飞起，最后才意识到：方程简洁和代码稳定之间隔着的不是“编程能力”，而是对离散格式、边界条件、时间推进稳定性的理解。这篇文章我就用经典的方腔顶盖驱动流（Lid-driven cavity）作为算例，把二维非定常Navier-Stokes从连续方程一路拆到可运行的MATLAB代码，顺便把我在调试过程中踩过、并且估计你也会踩的坑都列出来。适合正在学CFD、或者课程作业里要求用MATLAB自己写NS求解器的同学参考。

1. 这个示例到底解决了什么问题

二维非定常不可压缩Navier-Stokes方程，听起来很唬人，其实核心就两个量：速度场和压力场。但它难在速度的三个分量（二维是两个）不是独立变化的，它们必须满足不可压缩条件，也就是速度场始终无散。这个约束让方程从“能解”变成了“有条件的能解”，也催生出了一大堆数值处理方法。

1.1 方程里的每一项都对应哪段代码

先看最常用的原始变量形式：

其中u为速度矢量，p为压力，Re为雷诺数。左边的时间偏导和非线性对流项，右边是压力梯度项和粘性扩散项。加上不可压条件∇·u=0，这就是完整控制方程。

如果把这个方程和后面的代码对应起来，时间偏导对应时间推进循环，对流项对应速度场与涡量/速度梯度的乘积，扩散项对应二阶中心差分算子，压力梯度项则体现在速度修正和压力泊松方程里。之所以强调这个对应关系，是因为写代码时很容易“只更新了速度，忘了满足无散条件”，这是大多数发散问题的根源。

1.2 为什么选方腔顶盖流而不是圆柱绕流

很多教程喜欢用圆柱绕流展示NS求解器，但圆柱绕流要处理曲线边界、尾流涡街、卡门涡脱落，入门阶段直接上这个大概率劝退。方腔顶盖驱动流则简单很多：矩形计算域、边界都是直线、顶盖以固定速度拖动，其他壁面无滑移。即使网格和格式没那么完美，也能算出一个主涡旋，非常直观。

更重要的是这个算例有充分的文献数据和可视化结果可对比。Re=100时方腔内会出现一个稳定的主涡，涡心大约在坐标(0.62,0.74)附近；Re=1000时角涡变明显，流场会更复杂。所以它既是数值格式的验证平台，又是接触非定常流动的入门案例。

2. 数值方案：为什么绕开投影法

二维不可压NS最常见的数值求解思路是“投影法”，也叫做压力修正法：先不用压力，拿当前速度场显式推进一个中间速度，再求解压力泊松方程，最后用压力梯度修正速度，让速度重新满足无散条件。投影法物理清晰，但压力边界条件和网格排列处理起来比较麻烦，初学时容易在棋盘格振荡上浪费大量时间。

2.1 投影法与涡量流函数法的取舍

投影法的核心步骤可以写成三步：

1. 计算中间速度u*，不考虑压力梯度；
2. 求解压力泊松方程∇²p=∇·u*/Δt；
3. 用压力修正出满足无散条件的速度。

这种思路非常适合三维扩展，因为三维没有流函数这种简单的替代形式。但它也有个隐藏问题：如果压力和速度存储在同一个网格点上（非交错网格），压力泊松方程很容易产生棋盘格振荡，也就是压力场呈现出红黑交替的虚假模式，看起来像棋盘。要解决这个问题，要么使用MAC交错网格，要么用Rhie-Chow动量插值，对新手来说都是额外负担。

所以在这篇示例里，我选择的是涡量-流函数形式，本质上是二维NS方程的一种等价变换，通过取旋度把压力项消掉。方程变成涡量输运方程加流函数泊松方程：

只需要解一个抛物型方程和一个椭圆型方程，不需要处理压力边界条件。代价是无法直接得到压力场，后面如果需要压力，可以再用已知速度场后处理求解压力泊松方程。

2.2 边界涡量怎么给才不是乱写

用涡量流函数法，最容易出错的不是主方程，而是边界涡量。涡量的定义是ω=∂v/∂x−∂u/∂y，在无滑移边界上，涡量其实是由壁面切向速度的法向导数决定的。以方腔为例，顶盖速度为Uwall=1，其他壁面速度为0。

如果只用一阶单边差分，边界涡量的赋值方式是这样的：

• 底壁（y=0）：u=0，ω≈−u(2)/Δy；
• 顶壁（y=1）：u=Uwall，ω≈(u(末)−Uwall)/Δy；
• 左壁（x=0）：v=0，ω≈v(2)/Δx；
• 右壁（x=1）：v=0，ω≈−v(末)/Δx。

注意，这里的“第二行”和“倒数第二行”是指紧邻边界的内网格点。符号错了，流场会出现明显的不对称或直接发散。更精确的做法是用二阶单边差分，但在教学示例里一阶已经能稳定工作。

3. 直接上代码：一个能跑的MATLAB示例

这一段给出一个精简但完整的MATLAB实现思路。完整代码我拆成三个部分：初始化与网格准备、主时间循环、后处理与验证。网格尺寸取N=64，Re=100，初始流场全为零，顶盖突然启动。

3.1 准备网格和泊松求解器

N = 64; % 内部网格点数量 Re = 100; % 雷诺数 Uwall = 1; % 顶盖速度 h = 1 / (N + 1); % 网格间距 u = zeros(N+2, N+2); % 速度分量 u v = zeros(N+2, N+2); % 速度分量 v psi = zeros(N+2, N+2); % 流函数 omega = zeros(N+2, N+2);% 涡量 u(end, :) = Uwall; % 顶盖边界 A = gallery('poisson', N); % 离散拉普拉斯算子，Dirichlet边界

这里的关键是gallery('poisson', N)。它生成的是一个N²×N²的稀疏矩阵，对应正方形区域内部网格点的五点差分格式，边界值固定为0。因为流函数ψ在四壁都等于0，所以只需要对内部网格点求解即可。

需要提醒的是，这个矩阵对应的是“不带h²缩放”的拉普拉斯算子，所以在解泊松方程时右端要乘上h²，否则流函数会差一个网格步长的量级，速度场也会偏小。

3.2 主时间循环的三个关键动作

时间推进用一阶显式欧拉格式就能跑，空间导数全部用二阶中心差分。稳定性限制取对流CFL和扩散稳定性条件的保守值：

dt = min(0.5*h, 0.25*Re*h^2);

Re=100、N=64时，h≈0.0154，Re*h²≈0.0237，所以dt大概取0.005左右。如果雷诺数再大，扩散限制会更严格，时间步长会更小。

主循环如下：

nt = 2000; dt = min(0.5*h, 0.25*Re*h^2); for t = 1:nt % 1. 边界涡量 omega(1, :) = -u(2, :) / h; omega(end, :) = (u(end-1, :) - Uwall) / h; omega(:, 1) = v(:, 2) / h; omega(:, end) = -v(:, end-1) / h; % 2. 求解流函数 omega_in = omega(2:end-1, 2:end-1); psi_vec = A \ (h^2 * omega_in(:)); psi(2:end-1, 2:end-1) = reshape(psi_vec, N, N); % 3. 由流函数恢复速度 u(2:end-1, 2:end-1) = (psi(3:end, 2:end-1) - psi(1:end-2, 2:end-1)) / (2*h); v(2:end-1, 2:end-1) = -(psi(2:end-1, 3:end) - psi(2:end-1, 1:end-2)) / (2*h); % 重新施加边界速度 u(1, :) = 0; u(end, :) = Uwall; u(:, 1) = 0; u(:, end) = 0; v(1, :) = 0; v(end, :) = 0; v(:, 1) = 0; v(:, end) = 0; % 4. 涡量输运方程更新 dwx = (omega(2:end-1, 3:end) - omega(2:end-1, 1:end-2)) / (2*h); dwy = (omega(3:end, 2:end-1) - omega(1:end-2, 2:end-1)) / (2*h); adv = u(2:end-1, 2:end-1) .* dwx + v(2:end-1, 2:end-1) .* dwy; lap = (omega(2:end-1, 3:end) - 2*omega(2:end-1, 2:end-1) + omega(2:end-1, 1:end-2)) / h^2 + ... (omega(3:end, 2:end-1) - 2*omega(2:end-1, 2:end-1) + omega(1:end-2, 2:end-1)) / h^2; omega(2:end-1, 2:end-1) = omega(2:end-1, 2:end-1) + dt * (-adv + lap / Re); end

这个循环里最需要注意的一个细节是：速度u和v在边界上的值必须在每次循环里重新覆盖一次。虽然流函数边界恒为0，理论上内部速度自动满足无滑移，但你在计算内部速度时用到了中心差分，边界点本身并不会被更新，所以必须显式赋值。我见过不少版本在这里忘了加，结果边界条件慢慢失效，流场发散。

3.3 后处理与结果自检

算完之后，最直接的可视化是画涡量云图或者速度矢量图：

[xx, yy] = meshgrid(linspace(0,1,N+2)); figure; contourf(xx, yy, omega, 20); colorbar; title('Vorticity field');

也可以画流函数等值线，方腔顶盖流Re=100时看起来会有一个很清晰的顺时针主涡。如果想定量验证，可以提取竖直中心线上的速度剖面，和已有文献数据对比。Re=100时，沿着x=0.5这条线，u在靠近顶盖处接近1，在腔体中部转为负值，形成回流。这个剖面长什么样，和文献对得上，你的程序基本就是对的。

4. 实战避坑：发散、振荡、性能一个都别放过

这里集中说几个我在调试时踩过，并且很多初学者容易反复踩的问题。每个问题背后都有逻辑，知道了原因，调试起来就不会像无头苍蝇。

4.1 算着算着就发散，先查这三点

第一，检查时间步长是否满足稳定性条件。显式格式下，如果对流CFL数大于1，或扩散项νΔt/h²大于0.5，高频振荡会被逐步放大，然后快速发散。保守一点把dt再缩小一倍，如果流场稳定了，就是dt问题，否则继续查。

第二，检查泊松求解的缩放系数。很多人用gallery('poisson')时忘了乘h²，导致流函数和速度场整体缩小，看起来像是“没有流动”。反过来如果多乘了一个1/h²，速度场又会被放大几倍，也会触发数值失稳。

第三，检查边界涡量符号。特别是左右壁面的v导数方向，很容易因为索引顺序搞反。一个快速测试方法是把顶盖速度设为0，看静止流场是否真的保持不变；如果静止流场出现了涡量，说明边界涡量或速度恢复过程里肯定有错。

4.2 为什么没有压力棋盘格问题，以及怎么处理

在涡量流函数形式里，连续性约束已经通过流函数自动满足了，所以不存在压力棋盘格振荡。如果你用投影法求解原始变量形式，出现棋盘格是常事，这往往是因为速度和压力存储在同一个网格点上。解决办法要么换成MAC交错网格，要么对压力项使用Rhie-Chow插值，要么在压力泊松方程里用足够的五点格式并配合必要的滤波。

如果你只是需要一个能交作业的稳定示例，用涡量流函数法是最省心的。但如果后续课题需要三维和压力信息，那还是要回到投影法，到那个时候再学交错网格也不迟。

4.3 提速经验和下一步扩展

这个示例里，每个时间步都调用了A \ psi_vec，对N=64的网格毫无压力，但如果你把N加到128甚至256，这种“每步直接分解”的方式会明显变慢。更好的做法是在循环外先做一次矩阵分解：

[Lmat, Umat] = lu(A); % 循环内 psi_vec = Umat \ (Lmat \ (h^2 * omega_in(:)));

这样可以省掉每步重新分解矩阵的开销。更彻底的方案是用SOR迭代或者FFT求解泊松方程，前者内存占用小，后者在均匀网格上几乎可以做到O(N²logN)的速度。

如果你想把这个例子扩展下去，有几个方向：改成非定常圆柱绕流，需要引入浸没边界或者贴体网格；把时间推进换成RK3或RK4，在相同时间步长下可以获得更高精度；加入压力后处理，从已知速度场求解压力泊松方程，从而得到完整的NS解。每一步都不难，但都需要先把手上的这个基础示例彻底理解。

我在最初用这个涡量流函数框架时，犯过最无语的一个错误是初始化时忘了给顶盖边界速度赋值，导致算出来整个流场全是零。后来我把代码里“初始条件、边界条件、循环内更新、循环后重新施加边界条件”这四个阶段分开来检查，才慢慢养成调试CFD程序的节奏。如果你也是刚接触，建议先跑通这个示例，再手动改Re和网格数，感受一下非定常流场的变化规律，这比直接追着最新算法要扎实得多。

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