简介这是一套基于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 控制。典型现象是 Re100 时主涡靠近腔体上部Re≥1000 之后主涡中心明显下移右下和左下角开始出现二次涡。SIMPLE 算法在方腔流上的任务就是数值上解出这个压力场和速度场。边界条件需要特别注意顶盖为 u1、v0左右壁和下壁为 uv0。四个角点并不需要显式给定压力压力只是通过连续性方程被整体确定到相差一个常数。整个计算域的边界条件都无法直接提供压力 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(Nx1, Ny1); v zeros(Nx1, Ny1); p zeros(Nx, Ny); % 顶盖切向速度直接写在预留边界行 u(:, Ny1) 1.0; % 其余壁面速度已初始化为 0这段代码里 u 的列编号对应 y 方向u(:, Ny1) 是顶盖上的切向速度。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(i1,j); uS u(i,j-1); uN u(i,j1); % 计算 u(i,j) 两侧的对流通量 Fw rho * 0.5 * (u(i-1,j) u(i,j)) * dy; Fe rho * 0.5 * (u(i,j) u(i1,j)) * dy; Fs rho * 0.5 * (v(i,j-1) v(i1,j-1)) * dx; Fn rho * 0.5 * (v(i,j) v(i1,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(i1,j) aW·p(i-1,j) aN·p(i,j1) 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:Nx1,2:Ny)) * dy ... rho * (v(2:Nx,2:Ny) - v(2:Nx,3:Ny1)) * dx; for iter 1:20 for j 2:Ny-1 for i 2:Nx-1 pprime(i,j) (aE*pprime(i1,j) aW*pprime(i-1,j) ... aN*pprime(i,j1) 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 关键参数和建议值参数建议范围作用说明Re100 ~ 3200方腔流基准测试的常用范围Re 越高越考验网格分辨率Nx, Ny64×64 以上Re100 用 40×40 即可Re3200 至少 128×128alphaU, alphaV0.5 ~ 0.8动量方程欠松弛过大容易引起速度场震荡alphaP0.1 ~ 0.3压力修正的欠松弛和 SOR 超松弛分离残差阈值1e-5 ~ 1e-6基于最大质量通量归一化网格规模不要盲目加大。方腔流的最大可用雷诺数与网格尺度有关经验法则是网格雷诺数 rho·U·dx/mu 不超过 2。也就是说 Re100 用 64 网格时 dx1/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相当于人为削弱了压力对速度的修正量最终表现出来的结果是连续性残差永远压不下去。第三个坑是网格雷诺数过高。Re3200 在用 32×32 网格时哪怕一阶迎风也会在顶盖角落附近出现振荡。解决方式不是简单调低松弛而是加密网格或改用混合格式。可以先从 Re100 开始把主涡位置和中心线速度曲线跑出来再逐步提高 Re。最后建议每次改动参数后画一条方腔中心线 x0.5 上的 u 速度分布并与文献中的基准解对比。峰值位置和曲线形状对网格和收敛状态很敏感比只看残差更实用。本文还有配套的精品资源点击获取