企业官网建设流程全解析

简介：这是一套以MATLAB为平台、面向二维非定常Navier-Stokes方程数值求解的完整示例代码包，适用于流体力学相关专业的学生、科研工作者以及希望掌握CFD基础编程的工程师，能直观展示不可压缩流体在二维空间中的速度场与压力场随时间演化过程。压缩包共76个文件，包含44个.m格式的脚本源码，覆盖计算网格生成、基函数与矩阵组装、泊松方程压力修正、边界条件约束、时间推进等完整流程；另有29张PNG结果图便于查看流场模拟效果，以及HTML说明文档和license.txt许可说明，整体体积仅463KB，轻量而系统。目前已有1935人学习下载，可作为非定常流动数值模拟的入门范例或二次开发基础。代码中大量采用有限元/有限差分混合思路，并包含对流矩阵、粘性矩阵、载荷向量等核心子程序，结合quiver和contourf等可视化命令，能帮助使用者快速掌握从方程离散到程序实现的完整链路。通过修改初始条件、边界类型和网格参数，可扩展至不同非定常流动场景，适合教学演示、课程设计或科研预研使用。

1. 二维非定常N-S的MATLAB示例：为什么一上来就选涡量-流函数法

这份二维 Navier-Stokes 非定常求解的 MATLAB 示例，是我拆过的下载资源里最适合新手快速跑通不可压缩方腔流的一类。它没有直接解速度和压力，而是把方程组改写成涡量和流函数两个标量场，用显式时间推进加一个泊松迭代，二百行内就能出图。看起来像绕远了，实际上绕开了最麻烦的压力-速度耦合和棋盘格振荡，也把“非定常”这个关键词落在了一个能动态演示的算法上。

如果你是为了课程设计、仿真实验或者刚接触 CFD，想验证自己对差分格式的理解，这份资源很合适。它会回答“公式到底怎么变成矩阵循环、边界条件到底在哪一步写、dt 取多大才不会爆”。但如果你是做工业级湍流或者复杂几何，那这套均匀网格方腔示例只能当起点，不能当终点。

下面我从方程离散开始，一步步拆出能抄进 MATLAB 的代码，再告诉你哪些地方会翻车。

2. 方程与离散：把偏微分方程变成可循环的差分格式

2.1 控制方程与无量纲化

二维不可压缩 N-S 的原始变量形式是连续方程加两个动量方程，里面速度和压力耦合在一起，求解时要处理压力投影或者人工压缩性。这份资源采用涡量-流函数法，先把涡量定义为：

ω = ∂v/∂x - ∂u/∂y

再把流函数插进去，u = ∂ψ/∂y，v = -∂ψ/∂x。代入动量方程并取旋度，压力和梯度项会消掉，得到两个主方程：

∂ω/∂t + u ∂ω/∂x + v ∂ω/∂y = (1/Re)(∂²ω/∂x² + ∂²ω/∂y²)

∂²ψ/∂x² + ∂²ψ/∂y² = -ω

这两个方程对应两种完全不同的数学角色：涡量输运方程是抛物型/双曲型的混合，负责把流动的旋转“带”到下游；流函数泊松方程是椭圆型的，每一时刻都要解一次，用来从涡量恢复出速度场。因为不可压缩条件被流函数自动满足，所以数值上不需要再额外守一个连续方程残差。

资源里所有计算都是无量纲的。方腔边长取 1，顶盖速度取 1，特征时间就是 L/U0。这样雷诺数直接作为单参数出现，Re=100 时粘性效应比较明显，Re=1000 时可以看到角涡和薄一点的边界层。

写代码时不用刻意做单位换算，直接把 U0=1、Lx=1、Ly=1 放进去就行。

2.2 空间差分格式与边界条件

网格是均匀的矩形网格，nx 列，ny 行。对应关系在 MATLAB 里容易混，我统一约定：矩阵第一维是 y 方向，第二维是 x 方向。也就是说psi(i,j)表示第 i 个 y 节点、第 j 个 x 节点。

空间导数全部用二阶中心差分：

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

∂²ω/∂x² ≈ (ω(i,j+1) - 2ω(i,j) + ω(i,j-1)) / dx²

速度项从流函数差分得到：

u(i,j) ≈ (ψ(i+1,j) - ψ(i-1,j)) / (2dy)

v(i,j) ≈ -(ψ(i,j+1) - ψ(i,j-1)) / (2dx)

边界条件是方腔顶盖驱动流，四周流函数都固定为 0，速度无滑移。顶盖给 u=1，其余壁面速度为零。最麻烦的是涡量边界条件，它不是物理边界值，而是需要从无滑移约束推出来的数值边界条件。示例中用的是虚拟点导出的二阶格式：

ω(1,:) = -2(ψ(2,:) - ψ(1,:)) / dy² （底壁）

ω(ny,:) = -2(ψ(ny-1,:) - ψ(ny,:)) / dy² - 2U0/dy （顶盖）

ω(:,1) = -2(ψ(:,2) - ψ(:,1)) / dx² （左壁）

ω(:,nx) = -2(ψ(:,nx-1) - ψ(:,nx)) / dx² （右壁）

这个顶盖公式里的“-2U0/dy”就是移动壁面的贡献项，符号一旦写反，流函数中心涡旋方向会变成逆时针，后处理时“一眼就能看出不对劲”。我当年第一次跑的时候就是在这里翻的车。

提示：如果你的数组是zeros(ny,nx)，顶盖在最后一行，也就是omega(ny,:)，别写成omega(:,ny)，否则 MATLAB 不会报错，但结果完全错乱。

3. 时间推进与泊松求解：显式 Euler 到底行不行

3.1 涡量输运的时间步进

这份资源先用显式 Euler 做非定常推进，公式很简单：

omega_new = omega_old + dt * RHS

其中 RHS 是对流项和粘性项的组合。对应 MATLAB 代码，核心片段如下：

% 涡量输运右端项：对流 + 扩散 dwx(2:ny-1, 2:nx-1) = (omega(2:ny-1, 3:nx) - omega(2:ny-1, 1:nx-2)) / (2*dx); dwy(2:ny-1, 2:nx-1) = (omega(3:ny, 2:nx-1) - omega(1:ny-2, 2:nx-1)) / (2*dy); d2wx(2:ny-1, 2:nx-1) = (omega(2:ny-1, 3:nx) - 2*omega(2:ny-1, 2:nx-1) + omega(2:ny-1, 1:nx-2)) / dx^2; d2wy(2:ny-1, 2:nx-1) = (omega(3:ny, 2:nx-1) - 2*omega(2:ny-1, 2:nx-1) + omega(1:ny-2, 2:nx-1)) / dy^2; rhs = -u(2:ny-1, 2:nx-1).*dwx(2:ny-1, 2:nx-1) ... - v(2:ny-1, 2:nx-1).*dwy(2:ny-1, 2:nx-1) ... + (d2wx(2:ny-1, 2:nx-1) + d2wy(2:ny-1, 2:nx-1)) / Re; omega(2:ny-1, 2:nx-1) = omega(2:ny-1, 2:nx-1) + dt * rhs;

这里的dwx是涡量对 x 的导数，dwy是涡量对 y 的导数，d2wx、d2wy是二阶导数。代码只更新内点，边界点的涡量在另外的边界条件语句里写。注意u和v的尺寸是ny×nx，和omega一致，所以点乘时不需要转置。

显式 Euler 最大的问题不是精度，而是稳定性。粘性扩散项有步长限制，二维情况下近似为：

dt ≤ Re * dx² / 4

对流项限制近似为：

dt ≤ dx / max(u, v)

以 65×65 网格为例，dx≈0.015625，Re=100，粘性限制算出来是 dt ≤ 100×0.015625²/4 ≈ 0.0061。对流限制是 0.015625/1 ≈ 0.0156。所以 dt 取 0.002 是安全的，既不会爆，也不会慢到没耐心等。Re 提高后，粘性限制会成比例变松，但因为薄边界层需要对流限制更严格，建议至少乘个 0.5 的安全系数。

3.2 流函数泊松方程的逐点松弛

每一时间步更新完涡量，都要解一次泊松方程：

∂²ψ/∂x² + ∂²ψ/∂y² = -ω

这份资源用逐点 SOR 迭代求解，松弛因子 w 在 1 到 2 之间。MATLAB 函数如下：

function psi = poisson_sor(psi, omega, dx, dy, w, itmax, tol) % 逐点SOR求解流函数泊松方程 % psi: 上一时刻的流函数，作为迭代初值 % omega: 当前涡量 % w: 松弛因子，推荐1.5 % itmax: 最大迭代次数 % tol: 收敛残差 [ny, nx] = size(psi); dx2 = dx * dx; dy2 = dy * dy; beta2 = (dx/dy)^2; denom = 2*(1 + beta2); for k = 1:itmax psi_old = psi; % 内点更新 psi(2:end-1, 2:end-1) = ... (1-w) * psi(2:end-1, 2:end-1) + ... (w/denom) * ( psi(3:end, 2:end-1) + psi(1:end-2, 2:end-1) + ... beta2 * (psi(2:end-1, 3:end) + psi(2:end-1, 1:end-2)) + ... dx2 * omega(2:end-1, 2:end-1) ); % 边界流函数固定为0 psi(1,:) = 0; psi(end,:) = 0; psi(:,1) = 0; psi(:,end) = 0; err = max(abs(psi(:) - psi_old(:)), [], 'all'); if err < tol break; end end end

这个函数的关键是把上一时刻的psi作为迭代初值传进来。因为时间步长很小，相邻两步的流函数差别不大，所以实际运行中 SOR 往往 30 到 80 轮就收敛到 1e-6，不需要每次都从零开始算。如果设置成每次从零开始，计算量会大好几倍，但初学者很容易不知不觉这么写。

4. MATLAB 代码：从初始化到出图的六步实现

4.1 初始化与参数设置

下载资源里是一整套完整的ns2d_cavity.m，下面挑最核心的初始化部分讲。

clear; clc; % 物理与网格参数 Re = 100; % 雷诺数 U0 = 1; % 顶盖无量纲速度 Lx = 1; Ly = 1; % 方腔尺寸 nx = 65; ny = 65; % 网格点数 dx = Lx/(nx-1); dy = Ly/(ny-1); x = linspace(0, Lx, nx); y = linspace(0, Ly, ny); % 时间参数 dt = 0.002; % 时间步长 T = 5; % 总的无量纲时间 nt = ceil(T/dt); plot_step = round(0.1/dt); % 每0.1个时间单位画一帧 % 场变量初始化 psi = zeros(ny, nx); omega = zeros(ny, nx); u = zeros(ny, nx); v = zeros(ny, nx);

psi、omega、u、v四个变量全部是ny×nx的矩阵。nt在 dt=0.002、T=5 时是 2500 步，65×65 网格在普通笔记本上跑大约一两分钟。如果你想快速验证流程，可以先把网格改成 33×33，dt 相应增大到 0.005，几分钟推完，再用细网格出正式图。

4.2 主循环代码骨架

主循环的顺序是固定五步：从流函数算速度、处理速度边界、推进涡量、更新涡量边界、解泊松方程。顺序不能乱。

for n = 1:nt % 1. 由流函数计算内部点速度 u(2:end-1, 2:end-1) = (psi(3:end, 2:end-1) - psi(1:end-2, 2:end-1)) / (2*dy); v(2:end-1, 2:end-1) = -(psi(2:end-1, 3:end) - psi(2:end-1, 1:end-2)) / (2*dx); % 2. 速度边界：四壁无滑移，顶盖向右 u(end,:) = U0; u(1,:) = 0; v(1,:) = 0; u(:,1) = 0; v(:,1) = 0; u(:,end) = 0; v(:,end) = 0; v(end,:) = 0; % 3. 涡量输运右端项（见3.1代码，封装成vorticity_rhs） rhs = vorticity_rhs(u, v, omega, Re, dx, dy); omega(2:end-1, 2:end-1) = omega(2:end-1, 2:end-1) + dt * rhs; % 4. 涡量边界条件 omega(1,:) = -2*(psi(2,:) - psi(1,:)) / dy^2; omega(end,:) = -2*(psi(end-1,:) - psi(end,:)) / dy^2 - 2*U0/dy; omega(:,1) = -2*(psi(:,2) - psi(:,1)) / dx^2; omega(:,end) = -2*(psi(:,end-1) - psi(:,end)) / dx^2; % 5. 解流函数泊松方程 psi = poisson_sor(psi, omega, dx, dy, 1.5, 200, 1e-6); % 每plot_step步输出一帧 if mod(n, plot_step) == 0 % 后处理代码，见4.3 end end

第 4 步的涡量边界一定要在用完psi之后再更新。如果你把边界涡量放在psi更新之后去算，那相当于用旧涡量配合新流函数，时间上混了一拍，瞬态结果会失真。这个顺序对初学者来说非常容易搞反。

4.3 后处理与动画输出

资源包里还带了三个后处理片段，分别画流函数等值线、速度矢量和动画。

% 流函数等值线 figure; contourf(x, y, psi, 30); axis equal; colorbar; title('Streamfunction'); % 速度矢量 figure; quiver(x, y, u, v, 2); % 动画：每输出一帧刷新一次 figure; for k = 1:num_frames contourf(x, y, psi_stack(:,:,k), 30); axis equal; drawnow; end

注意contourf的第一、二个参数是 x 和 y 向量，第三个参数是ny×nx的矩阵，顺序反了会直接把坐标轴搞混。速度矢量的quiver也一样，u和v的尺寸要与 x、y 匹配。画动画时如果发现流速越来越慢，大概率是drawnow卡了预分配，后续建议先把所有帧存成三维矩阵，最后统一重绘。

5. 避坑指南：非定常 N-S 求解最常见的六个翻车现场

5.1 现象：NaN 满天飞，程序跑到一半戛然而止

原因：显式 Euler 的 dt 超过了粘性稳定极限。特别是 Re 不太高、网格比较密的时候，粘性项对 dt 的限制是主导。Re=100、dx=0.0156 时，理论限制是 0.0061，但如果你用 129×129 网格，dx 变成 0.0078，限制直接缩到 0.0015，原来的 dt 一定爆。

解决：把 dt 先缩到理论值的 0.5 倍，比如 dx=0.0078 时取 dt=0.0007。如果问题照旧，再检查是不是初始涡量边界给了一个巨大的尖峰，这时把顶盖速度改成前几步平滑加载，比如 U0_eff = U0 * min(t/0.1, 1)，能明显减轻瞬时冲击。

5.2 现象：流函数画出来像“梳子”，锯齿状

原因：边界涡量符号写反，或者顶盖移动方向与流函数梯度的耦合对不上。出现锯齿的另一个常见原因是速度差分时用了向前差分，而不是中心差分。

解决：先确认顶盖涡量是-2*(psi(end-1,:)-psi(end,:))/dy^2 - 2*U0/dy，不是加号。再用短时间步跑 200 步，看quiver里靠近顶盖的速度是不是指向右。如果指向左，就把 U0 的符号反过来，或者检查u和v差分公式里的负号。这里没有玄学，方向错了就是错了。

5.3 现象：泊松迭代一千次还不收敛

原因：大部分情况是松弛因子取 1，或者边界流函数没有固定成 0。逐点更新时如果psi(1,:)没有强制赋零，边界值会被内点迭代慢慢“拖走”，残差永远下不去。

解决：把 SOR 的 w 取 1.5，收敛残差放宽到 1e-6，迭代上限 200。每次传入上一时刻的 psi 作为初值，残差一般能在几十轮内降下来。如果 200 轮还不行，优先检查边界赋值语句，再检查网格步长dx、dy是否在分母上写反。

5.4 现象：速度场左右不对称，左边明显比右边大

原因：MATLAB 的矩阵索引和 x/y 坐标关系没理顺。矩阵第一维是行、第二维是列，行对应 y，列对应 x。如果你在差分代码里把psi(2:end-1, 3:end)和psi(1:end-2, 3:end)混用，x 方向导数的截断位置就偏移了，左右自然不对称。

解决：写代码之前先打印size(psi)，并在注释里写清“第 1 维是 y”。可以做一个快速自检：把顶盖速度设成 0，任其衰减，理论上每个点速度都应该是 0；如果出现非零速度，说明差分模板的索引偏移了。

5.5 现象：残差降不下去，等值线一直飘

原因：显式 Euler 的数值耗散太大，非定常推进的瞬态响应还没有充分衰减；或者 T 设置得太短，中心涡量还没稳定下来。

解决：资源包附带 RK2 版本，用中点格式可以让 dt 适当放大。改用omega_mid = omega + dt*rhs1，再用omega = omega + 0.5*dt*(rhs1+rhs2)做一次校正。每一轮中点计算也要调用一次泊松和边界更新，代码会变长，但瞬态曲线的抖动会明显小很多。

5.6 现象：Re=1000 之后结果抖动，出现非物理振荡

原因：网格分辨率不足。方腔右上角和左下角的旋涡尺度随 Re 变大而变小，65×65 网格在 Re=1000 时根本无法分辨边界层内的速度反转。这不是数值格式的问题，是分辨率不够。

解决：粗略估计网格数要跟随 Re 的平方根增长，Re=1000 用 129×129，Re=5000 至少 257×257。若资源包里的默认网格是 65×65，请按这个规则手动改。

6. 验证与进阶：用网格收敛性检验你的求解器

6.1 网格收敛性：在同一 Re 下比较速度剖面

只看到一个漂亮的涡流图不代表求解器正确，我习惯用中线速度剖面做验证。取 x=0.5 处的 u(y) 曲线，用 65×65 和 129×129 两套网格分别跑，比较两条曲线的偏差。收敛的二阶格式应该随网格加密，误差大约按 1/4 的比例下降。

mid = (nx+1)/2; % 奇数网格，x=0.5 正好在节点上 u_mid = u(:, mid); % 纵向速度剖面 plot(u_mid, y, 'linewidth', 1.5);

判断标准是：主涡附近的速度峰值位置一致，靠近顶盖处速度从 1 迅速衰减到负值，下壁附近回流速度幅度合理。如果剖面曲线像波浪一样，多半是涡量边界条件没收敛。

6.2 进阶扩展：RK2 与三维扩展

在资源包的进阶版里，显式 Euler 被换成了二阶 Runge-Kutta。关键改动是把一次推进变成两次半步，并且每一步都要更新涡量边界和流函数：

rhs1 = vorticity_rhs(u, v, omega, Re, dx, dy); omega_mid = omega + dt * rhs1; % 用omega_mid重新调一次泊松，得到psi_mid，再算速度 rhs2 = vorticity_rhs(u_mid, v_mid, omega_mid, Re, dx, dy); omega = omega + 0.5 * dt * (rhs1 + rhs2);

这套思路放大到三维并不难，但第三步是通用的：涡量输运方程改成三维形式，泊松方程也要改成三维拉普拉斯，网格数会从平方涨到立方。所以我一般建议先把二维验证做透，再往三维推。从那以后，我每次换网格或者换雷诺数，都会先跑一段短时程中心涡量时间序列，等它平滑趋于定值后再开正式的长时程计算。这个习惯帮我挡掉了不知道多少次“看起来还行、其实已经发散”的假象，希望也能帮到你。

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