MATLAB实现SIMPLE算法求解方腔驱动流：从方程到收敛

简介：这是一套基于MATLAB开发的二维方腔驱动流模拟程序，主要面向计算流体力学初学者、工程技术人员以及需要快速实现SIMPLE算法代码的开发者。资源以SIMPLE压力修正框架为内核，将方腔流动的网格离散、边界条件施加、代数方程求解与速度压力修正整合为清晰流程，能够直观展示不可压缩N-S方程在标准网格下的数值求解思路。压缩包共11个文件，其中10个为.m格式的MATLAB脚本，1个为txt说明文件，整体大小约8KB。脚本按功能划分明确，涉及系数矩阵生成、水平与垂直动量方程更新、五对角方程组求解、压力修正、速度修正以及无散度检查等关键环节，读者可对照源码逐模块学习SIMPLE算法的迭代逻辑，并理解方腔流算例中无滑移壁面条件的处理方式。目前已有284人学习下载，适合用于课程实验、毕业设计或作为进一步研究对流、湍流等复杂流动问题的起点。

1. 为什么说方腔驱动是MATLAB里验证SIMPLE算法的最佳起手式

如果你刚拿到一个叫 cavityFlow2D 的 MATLAB 包，第一反应多半是“顶盖一拖，半小时之内把涡画出来”。方腔驱动流确实没有复杂几何，但它在验证 simple 算法上的作用并不轻：顶盖以恒定速度拖动封闭方腔内的流体，流动在腔内形成主涡和两个角涡，整个现象只看 Re 一个参数。SIMPLE 算法最难的环节是压力与速度双向耦合，而这个案例恰好把不可压缩 Navier-Stokes 方程从离散到收敛的全部要素都覆盖了：交错网格、通量平衡、压力修正矩阵、欠松弛迭代。除了 CFD 入门的在校生，这类代码对做工业仿真的工程师也有参考价值，因为它很适合作自己求解器的最小验证床。

下面按“先立理论、再写代码、然后调参数、最后看残差排错”的顺序展开，内容即使没有那套 zip 包，也能照着在 MATLAB 里重建一遍。

2. SIMPLE算法解方腔流动：先建立方程和离散模型

2.1 方腔流控制方程和边界条件的无量纲形式

方腔流动的物理域很简单：一个边长为 1 的方腔，顶盖以速度 U0 向右移动，其余三面固定。选择 U0 和方腔边长 L 作为基准量，压力用 ρU0² 归一化，定常二维不可压流动的 Navier-Stokes 方程可以写成

∂u/∂x + ∂v/∂y = 0

∂(uu)/∂x + ∂(uv)/∂y = -∂p/∂x + (1/Re)(∂²u/∂x² + ∂²u/∂y²)

∂(uv)/∂x + ∂(vv)/∂y = -∂p/∂y + (1/Re)(∂²v/∂x² + ∂²v/∂y²)

其中 Re = U0L/ν。无量纲化之后，流场形态只受 Re 控制。典型现象是 Re=100 时主涡靠近腔体上部，Re≥1000 之后主涡中心明显下移，右下和左下角开始出现二次涡。SIMPLE 算法在方腔流上的任务就是数值上解出这个压力场和速度场。

边界条件需要特别注意：顶盖为 u=1、v=0，左右壁和下壁为 u=v=0。四个角点并不需要显式给定压力，压力只是通过连续性方程被整体确定到相差一个常数。整个计算域的边界条件都无法直接提供压力 Dirichlet 条件，这对后面压力修正方程的结构有直接影响。

2.2 SIMPLE算法的预测-修正循环

SIMPLE 即 Semi-Implicit Method for Pressure Linked Equations。它不直接解原始联立方程，而是把每个迭代周期拆成“先猜压力解动量、再用连续性方程修正压力”的两步：

1. 给定初场 u, v, p 2. 用当前压力 p* 求解动量方程，得到中间速度 u*, v* 3. 构造速度修正 u' = d_e(p'_P - p'_E)，v' = d_n(p'_P - p'_N) 4. 把修正速度代入连续性方程，得到压力修正方程 5. 解出 p'，计算 p = p* + αp p'，再用 p' 更新 u、v 6. 回到第 2 步，直到质量通量残差低于阈值

这里的 d_e 和 d_n 来自动量方程系数，物理意义是“单位压力差能产生的速度修正量”。第 4 步得到的不是原始泊松方程，而是把速度修正消去后得到的只含 p' 的五点方程。整个算法的半隐式特征就体现在这里：新求出的压力 p' 被用于更新速度，但速度修正没有反过去影响动量方程系数，因此每个迭代周期内都不需要重新组装动量矩阵。

实际写 MATLAB 代码时，不要把这个过程组织成一个大函数，最好拆成“解动量”“算残差”“解压力修正”“回代修正”四个独立函数。这样调试时能单独看每一步的残差。

2.3 为什么交错网格是SIMPLE的默认选择

如果压力和速度都定义在同一个网格中心，SIMPLE 算法会出现棋盘式压力分布。原因是同位网格上压力梯度由相隔两个网格的节点差分得到，常数倍的正负交替压力场不会被压力梯度感知，于是压力场在迭代中很难被消除振荡。

交错网格的做法是把速度放在网格面上：u 放在 x 方向的垂直面上，v 放在 y 方向的水平面上，压力放在单元中心。这样压力梯度直接用相邻两个压力值相减得到，间隔只有 Δx 或 Δy，棋盘模态立刻会被压力修正方程压制。付出的代价是变量数组尺寸不一致、边界处理变繁琐，但这在方腔流这种结构化网格上完全可控。

SIMPLE 在 MATLAB 里最常见的实现方式也是交错网格。拿到 cavityFlow2D 这类代码时，第一件事应该是确认它的压力、u、v 三个数组的尺寸关系，而不是急着运行。

3. 用MATLAB从零搭一个可复现的SIMPLE求解核心

3.1 按交错网格定义数组和速度边界

方腔流适合用均匀网格。设压力单元数为 Nx×Ny，为了方便施加顶盖速度边界，可以在 u、v 两个方向都预留边界节点：

Nx = 40; Ny = 40; dx = 1/Nx; dy = 1/Ny; Re = 100; rho = 1; mu = 1/Re; % u,v 均保留边界节点；p 在单元中心 u = zeros(Nx+1, Ny+1); v = zeros(Nx+1, Ny+1); p = zeros(Nx, Ny); % 顶盖切向速度直接写在预留边界行 u(:, Ny+1) = 1.0; % 其余壁面速度已初始化为 0

这段代码里 u 的列编号对应 y 方向，u(:, Ny+1) 是顶盖上的切向速度。v 在顶盖上是法向速度，所以仍然保持 0。角点速度虽然没有直接出现在离散方程里，但在输出流场时可能要用到，通常用相邻节点的平均值代替，不要在角点上随意赋非零值。

这里的数组尺寸是一种扩展交错网格写法，它与典型 CFD 教材里的 p 中心、u 面、v 面定义等价，但保留边界节点后索引更直观。实际写动量方程时，只需要把内部速度节点当作未知量，外部边界节点则作为已知值参与系数计算。

3.2 动量方程的系数组装与亚松弛

对二维均匀网格上的速度节点 u(i,j)，离散后的动量方程可以写成标准形式

aP·u(i,j) = aW·uW + aE·uE + aS·uS + aN·uN + (p(i-1,j)-p(i,j))·dy

其中扩散系数是 mu·dy/dx，对流系数需要处理迎风方向。下面是一段以 x 方向速度为例的 MATLAB 参考片段，采用一阶迎风混合项：

for j = 2:Ny for i = 2:Nx uW = u(i-1,j); uE = u(i+1,j); uS = u(i,j-1); uN = u(i,j+1); % 计算 u(i,j) 两侧的对流通量 Fw = rho * 0.5 * (u(i-1,j) + u(i,j)) * dy; Fe = rho * 0.5 * (u(i,j) + u(i+1,j)) * dy; Fs = rho * 0.5 * (v(i,j-1) + v(i+1,j-1)) * dx; Fn = rho * 0.5 * (v(i,j) + v(i+1,j)) * dx; Dw = mu * dy / dx; De = mu * dy / dx; Ds = mu * dx / dy; Dn = mu * dx / dy; % 迎风贡献：只加入上风方向系数 aW = Dw + max(Fw, 0); aE = De + max(-Fe, 0); aS = Ds + max(Fs, 0); aN = Dn + max(-Fn, 0); aP = aW + aE + aS + aN + (Fe - Fw + Fn - Fs); u(i,j) = (aW*uW + aE*uE + aS*uS + aN*uN ... + (p(i-1,j) - p(i,j)) * dy) / aP; end end

代码中 max(Fw,0) 表示当流动方向为正时，西侧邻居对当前节点的影响更强；max(-Fe,0) 同理。aP 里加上的净通量项保证了离散方程在均匀流场中能精确成立。上步解出的 u 还需要做亚松弛，常见做法是

u_new = alphaU * u_calc + (1 - alphaU) * u_old

不要把这个松弛直接加到 aP 上，否则残差和收敛曲线会变得很难读。alphaU 通常取 0.5 到 0.8。

3.3 压力修正方程与SOR求解

动量方程解完后，当前速度场一般不满足连续性。每个单元的质量不平衡量 b(i,j) 是压力修正方程的源项。压力修正方程的标准形式为

aP·p'(i,j) = aE·p'(i+1,j) + aW·p'(i-1,j) + aN·p'(i,j+1) + aS·p'(i,j-1) + b(i,j)

在 MATLAB 里不必组装全局稀疏矩阵，直接做几轮 SOR 扫描即可：

b = zeros(Nx, Ny); % 由交错网格速度计算单元净质量流出量 b(2:Nx, 2:Ny) = rho * (u(2:Nx,2:Ny) - u(3:Nx+1,2:Ny)) * dy ... + rho * (v(2:Nx,2:Ny) - v(2:Nx,3:Ny+1)) * dx; for iter = 1:20 for j = 2:Ny-1 for i = 2:Nx-1 pprime(i,j) = (aE*pprime(i+1,j) + aW*pprime(i-1,j) + ... aN*pprime(i,j+1) + aS*pprime(i,j-1) + b(i,j)) / aP; pprime(i,j) = omega * pprime(i,j) + (1-omega) * pprime(i,j); end end end % 速度修正，d 来自动量系数 u(2:Nx,2:Ny) = u(2:Nx,2:Ny) + d_u .* (pprime(1:Nx-1,2:Ny) - pprime(2:Nx,2:Ny)); v(2:Nx,2:Ny) = v(2:Nx,2:Ny) + d_v .* (pprime(2:Nx,1:Ny-1) - pprime(2:Nx,2:Ny)); % 压力修正：只在这里使用 alphaP p = p + alphaP * pprime;

这里的 aE、aW、aN、aS 由质量流量和 d_u、d_v 组成，aP 是它们的和。注意 pprime 是一个 Nx×Ny 数组，速度修正时索引要取相邻压力差。SOR 的超松弛 omega 一般取 1.2 到 1.8，但它是用于内迭代收敛的，跟 SIMPLE 外循环的 alphaP 不是一回事。alphaP 过大时压力修正会过冲，流场容易出现周期性振荡。

4. 在MATLAB里把方腔流跑起来：参数设置、迭代流程与后处理

4.1 主迭代循环

把上一章的模块串成主循环，推荐写成如下形式：

for iter = 1:maxIter u_old = u; v_old = v; % 每个外迭代周期内解 2~3 次动量方程 for inner = 1:3 u = solve_momentum_u(u, v, p, rho, mu, dx, dy, alphaU); v = solve_momentum_v(u, v, p, rho, mu, dx, dy, alphaV); end b = compute_mass_residual(u, v, rho, dx, dy); pprime = solve_pressure_correction(b, u, v, rho, dx, dy); [u, v, p] = correct_velocity_pressure(u, v, p, pprime, alphaP); res_u = norm(u(:) - u_old(:)) / (norm(u(:)) + eps); if res_u < 1e-6 break; end end

每个函数都接收同组网格参数，避免在多个脚本之间复制全局变量。动量方程内部迭代次数并不需要太多，因为 pressure-velocity 耦合由外层的压力修正来收敛。如果发现残差曲线平走，先检查单次动量迭代是否已经收敛。

4.2 关键参数和建议值

网格规模不要盲目加大。方腔流的最大可用雷诺数与网格尺度有关，经验法则是网格雷诺数 rho·U·dx/mu 不超过 2。也就是说 Re=100 用 64 网格时 dx=1/64，网格雷诺数只有 1.56，小于 2，比较安全。Re 提高到 1000 时，同样的 64 网格网格雷诺数超过 15，一阶迎风还能勉强跑，二阶中心差分会直接发散。

4.3 用流函数和速度矢量验证流场结构

等值线比箭头图更能说明方腔流的涡结构。交错网格的 u、v 不在同一组节点上，做可视化前需要先插值到同一组坐标：

[X, Y] = meshgrid(linspace(dx/2, 1-dx/2, Nx), ... linspace(dy/2, 1-dy/2, Ny)); xu = (0:Nx) * dx; % u 节点 x 坐标 yu = (0:Ny) * dy; % u 节点 y 坐标 [Ux, Uy] = meshgrid(xu, yu); uq = interp2(Ux', Uy', u', X', Y'); [Xv, Yv] = meshgrid(xu, yu); vq = interp2(Xv', Yv', v', X', Y'); streamslice(X', Y', uq, vq, 3); axis equal; xlim([0 1]); ylim([0 1]);

在 MATLAB 里交错网格的数组通常会有一维尺寸比其他方向多 1，直接用 streamslice 前必须插值。更简单的验证方式是直接画速度矢量图，然后手动找零速度点；但从残差和涡心位置判断版本一致性，还是流线图更直观。

5. 收敛判断与三个容易踩的坑

5.1 用质量通量残差而不是速度变化判断收敛

很多初学实现会把残差定义为两次迭代的速度差，但速度差减小不代表连续性方程被满足。更可靠的指标是每个单元的质量净流出量：

Rmean = mean(abs(b(:))) / rho; Rmax = max(abs(b(:))) / rho;

Rmean 反映总体质量守恒水平，Rmax 反映局部最大不平衡。正常的 SIMPLE 迭代曲线应该是前几十步快速下降，之后进入缓慢下降段。如果 Rmean 在某个值附近震荡，通常是 alphaP 太大；如果阶梯状上升，多半是压力修正方程中的 SOR 内迭代没有收敛。

5.2 三个容易踩的坑

第一个坑是压力修正方程没有固定参考点。压力只有 Neumann 边界条件，压力修正矩阵是奇异的。SOR 或者 MATLAB 自带的迭代求解器不会直接报错，但 pprime 会整体飘移。常见做法是每轮 SOR 后对 pprime 减去平均值，或者在求解前把 p(1,1) 固定为 0。

第二个坑是把 alphaP 同时用于速度修正。SIMPLE 推导过程中，速度修正和压力修正之间是强耦合的，速度修正是由 p' 直接算出来的。如果在速度修正上再乘一个 alphaP，相当于人为削弱了压力对速度的修正量，最终表现出来的结果是连续性残差永远压不下去。

第三个坑是网格雷诺数过高。Re=3200 在用 32×32 网格时，哪怕一阶迎风也会在顶盖角落附近出现振荡。解决方式不是简单调低松弛，而是加密网格或改用混合格式。可以先从 Re=100 开始，把主涡位置和中心线速度曲线跑出来，再逐步提高 Re。

最后建议每次改动参数后，画一条方腔中心线 x=0.5 上的 u 速度分布，并与文献中的基准解对比。峰值位置和曲线形状对网格和收敛状态很敏感，比只看残差更实用。

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