企业官网建设流程全解析

简介：这是一份基于SIMPLE算法、采用交错网格求解非定常不可压缩粘性流体方程的Matlab源码资源，主要面向流体数值计算方向的本硕学生与科研人员，可帮助理解压力速度耦合迭代、交错网格离散以及瞬态时间推进等核心机制。资源包共包含7个文件，其中2个m函数源文件负责核心计算，2个txt说明文件辅助解释，2个png图像展示运行结果，还有1个md文档梳理使用说明，整体大小约470KB，体积轻巧、结构清晰。源码内含主解算器与高阶插值辅助函数，配合运行结果图示，可直接在Matlab中运行并对照分析，便于逐行理解SIMPLE算法的迭代流程、边界条件处理以及交错网格的优势。目前已有104人学习浏览该资源，适合作为本科或硕士阶段的教研学习材料，也可作为课程设计、毕业设计或科研预研的参考框架，为后续开发更复杂的不可压缩流场求解程序打下基础。

1. 为什么说这套 Matlab 源码先用的是 SIMPLE 而不是时间精度

如果只看命名，“非定常 Navier-Stokes 解算器”会让人觉得重点在时间推进格式上，但真正值得拆的其实是 SIMPLE 算法里的压力修正环节。很多第一次写 NS solver 的人都栽在同一处：速度场看着合理，压力场却像棋盘格一样一高一低地交替闪烁，因为常规网格上的中心差分会把奇偶节点之间的信息切断。这套源码通过把 u、v 放在单元面上、把压力放在单元中心，从几何上掐断了这种奇偶失耦的路径，再配合压力修正把速度场投影到无散空间。它适合正在学有限体积法和压力修正算法的研究生，也适合想快速验证自己网格格式的工程师。下面从算法原理、主程序推进、插值函数和参数调优几个层面拆解。

2. SIMPLE 压力修正与交错网格的动态解耦

2.1 非交错网格上的棋盘压力场问题

非交错网格把所有变量都存储在同一套节点上。对一维压力梯度采用中心差分时，(p(i+1)-p(i-1))/(2*dx)完全跳过 i 节点自身的压力值。于是压力场可以呈现1, 0, 1, 0这种锯齿形分布，而梯度计算得出的速度没有被污染，方程组照样收敛。问题是节点 i 处根本没有关于 p(i) 的约束，这个锯齿场只是被矩阵隐式地保存了下来。

交错网格的解法是把速度分量放在控制体的面中心。对于二维顶盖驱动方腔流，u 放在左右面的中点，v 放在上下面的中点，p 放在体心。这样压力梯度dp/dx在 u 方程中用的是相邻两个压力节点的差，dp/dy在 v 方程中用的是相邻压力节点的差，每个方向的梯度都直接作用在同一条网格线上，不再跳过中间点。

2.2 交错网格的变量布置与通量几何

建网格时不能只生成一套坐标。由于变量位置不同，需要同时维护标量坐标和两个速度坐标：

这段网格生成代码演示了坐标数组的组织方式：

% 以简单均匀网格为例，生成交错网格坐标 nx = 64; ny = 64; Lx = 1.0; Ly = 1.0; dx = Lx / nx; dy = Ly / ny; % 标量（压力）单元中心网格 for j = 1:ny for i = 1:nx xc(i,j) = (i - 0.5) * dx; yc(i,j) = (j - 0.5) * dy; end end % u 速度点：位于标量网格的垂直界面（x 方向偏移半个 dx） for j = 1:ny for i = 1:nx+1 xu(i,j) = (i - 1) * dx; yu(i,j) = (j - 0.5) * dy; end end % v 速度点：位于标量网格的水平界面（y 方向偏移半个 dy） for j = 1:ny+1 for i = 1:nx xv(i,j) = (i - 0.5) * dx; yv(i,j) = (j - 1) * dy; end end

代码中nx+1和ny+1是交错网格最容易写错的地方。u 数组在 x 方向比压力数组多一个点，v 数组在 y 方向多一个点，因为边界面上的速度必须存在才能施加无滑移条件。实际求解时，我一般直接把 u 定义为(nx+1) × ny的数组，v 定义为nx × (ny+1)，而不是把所有量都开成同样大小的矩阵，这样索引关系更直观，也避免在循环里反复换算。

2.3 压力修正方程的本质与系数组装

SIMPLE 的预测步先用当前压力场求解动量方程，得到不满足连续性的速度场u*、v*。修正步把修正量p'引入动量方程，得到只与p'相关的速度修正公式。把修正后的速度代入连续性方程，就得到关于p'的五点离散椭圆方程：

% 压力修正方程系数（SIMPLE，二维均匀网格，不可压流） % 系数来自动量离散中的耦合项，A_p 是动量方程对角系数 ap_E = rho * dy / dx; % 左侧邻居贡献系数 ap_W = rho * dy / dx; % 右侧邻居贡献系数 ap_N = rho * dx / dy; % 下方邻居贡献系数 ap_S = rho * dx / dy; % 上方邻居贡献系数 ap_P = ap_E + ap_W + ap_N + ap_S; % 质量残差源项：由 u*、v* 的净通量组成 b = (u_star(i,j) - u_star(i-1,j)) * dy ... + (v_star(i,j) - v_star(i,j-1)) * dx;

这个方程的对角系数永远是周围四个系数之和，因此系数矩阵对称正定，用共轭梯度或直接法都能稳定求解。b项的含义是当前速度场通过控制体积六个面（二维为四边）的净流出量，它不为零时才需要压力修正。当b趋近于机器零时，速度场已经满足连续方程，修正步也就无事可做了。

3. NavierStokesUnsteady.m 的非定常推进流程

3.1 时间项离散与时间步长约束

非定常求解的关键在于时间步长不能随意取。动量预测步如果采用显式离散，稳定条件同时受对流和扩散限制。对于均匀网格，常见做法是同时检查两个条件：对流 CFLu*dt/dx通常要小于 0.5，扩散数nu*dt/dx^2要小于 0.25。实际代码里的步长选择可以写成这样：

% 根据当前最大速度和粘性系数估计时间步长 umax = max(max(abs(u(2:end-1,2:end-1)))); vmax = max(max(abs(v(2:end-1,2:end-1)))); cfl_max = 0.4; dt_conv = cfl_max * min(dx, dy) / max(umax, vmax); dt_diff = 0.2 * min(dx, dy)^2 / nu; dt = min(dt_conv, dt_diff);

dt_conv对高雷诺数流动起主导作用，因为速度越大，跨越一个网格所需时间越短；dt_diff对低雷诺数流动起约束作用。显式推进的优势是每步计算量小，缺点是雷诺数较高时dt被压得很小，导致总步数剧增。如果你要跑高雷诺数算例，可以把动量方程里的扩散项或压力泊松方程做隐式处理，但这份源码的定位是教学演示，显式时间推进配合 SIMPLE 已经足够展示完整流程。

3.2 动量预测步的对流项与粘性项离散

动量预测步在每个速度点上分别计算。以水平速度u为例，离散式包含非稳态项、对流项、扩散项和压力梯度项。对流项需要使用界面上的速度值，这个工作交给 higherorderinterpolation.m 完成：

% NavierStokesUnsteady.m 主循环结构示意 for n = 1:nt % 保存上一时刻速度用于残差统计 u_old = u; v_old = v; % ---- 动量预测步：从 p_n 求 u_star, v_star ---- for j = 2:ny for i = 2:nx ue = higherorderinterpolation(u, i, j, dx, 'quick'); uw = higherorderinterpolation(u, i-1, j, dx, 'quick'); u_conv = (ue^2 - uw^2) / dx; u_diff = nu * (u(i+1,j) - 2*u(i,j) + u(i-1,j)) / dx^2 ... + nu * (u(i,j+1) - 2*u(i,j) + u(i,j-1)) / dy^2; u_star(i,j) = u(i,j) + dt * (-u_conv + u_diff - (p(i+1,j)-p(i,j))/rho/dx); end end % ---- 压力修正方程求解 ---- p_corr = pressure_solver(ap, b, nx, ny); % ---- 速度修正：把无散度修正投影回去 ---- u(2:nx,2:ny) = u_star(2:nx,2:ny) - dt * ... (p_corr(3:nx+1,2:ny) - p_corr(2:nx,2:ny)) / rho / dx; end

这里的higherorderinterpolation传入的是整个速度场和界面索引。注意在交错网格上，u的控制体积中心本身就在界面上，所以计算u^2时还需要把速度插值到u控制体积的左右界面上，这是后面插值函数存在的意义。动量预测步结束后，u_star已经包含压力梯度的作用，但这一压力来自上一时间步，不满足当前时刻的连续条件。

3.3 压力泊松方程的矩阵求解

压力修正方程在二维规则网格上是一个五对角矩阵。Matlab 里最常见的做法是直接构造稀疏矩阵后用反斜杠求解：

% 组装稀疏矩阵与右端项 A = sparse(N, N); bvec = zeros(N, 1); for j = 2:ny-1 for i = 2:nx-1 idx = (j-1)*nx + i; A(idx, idx) = ap_P; A(idx, idx-1) = -ap_W; A(idx, idx+1) = -ap_E; A(idx, idx-nx) = -ap_S; A(idx, idx+nx) = -ap_N; bvec(idx) = bval(i,j); end end % 求解后重新整形为二维数组 p_corr_vec = A \ bvec; p_corr = reshape(p_corr_vec, nx, ny);

sparse在这里是必须的，因为全稠密矩阵在 100×100 网格上就有 1e8 个元素，直接爆内存；而五对角非零元只有约 5e4 个。A \ b对稀疏 SPD 矩阵会自动选择 Cholesky 分解，收敛速度比迭代法更稳定。源包中如果没有单独的压力求解函数，通常会把这段逻辑直接放在NavierStokesUnsteady.m的循环里。

4. higherorderinterpolation.m 的高阶插值与松弛参数调优

4.1 对流界面插值为什么不能只用线性

SIMPLE 算法本身不规定对流项怎么离散。如果界面速度用简单算术平均(u(i,j)+u(i+1,j))/2，在网格 Peclet 数较大时会出现非物理振荡。一阶迎风虽然绝对稳定，但数值耗散大，算出来的方腔流涡心位置会偏移。higherorderinterpolation.m 存在的意义，就是在迎风和中心差分之间取平衡。最典型的选择是 QUICK 格式，它在界面上游方向取三个节点进行二次插值：

% higherorderinterpolation.m 对界面通量的 QUICK 插值示意 % phiC 为当前节点，phiU 为上游邻居，phiD 为下游邻居 % 实际调用时根据流向自动交换 U/D 的角色 function phiI = higherorderinterpolation(phiC, phiU, phiD) % QUICK 插值权重 phiI = 0.5 * (phiC + phiD) - 0.125 * (phiD - 2*phiC + phiU); end

这里的核心是第二项中的二阶差分修正：它有迎风倾向，但比纯迎风少了一个数量级的数值耗散。调用方需要根据界面法向速度的正负判断谁是上游。实现时常见做法是判断ue>0则phiU=phi(i-1,j)，否则phiU=phi(i+2,j)。QUICK 的一个缺点是系数矩阵不再保持对角占优，但在 MATLAB 直接法下这不是问题；如果你使用点迭代求解器，需要在迭代中保证收敛。

4.2 松弛因子与残差收敛判定

SIMPLE 里的压力修正一般不能全量加上，否则修正量过大导致下一迭代步发散。修正速度时也要做欠松弛处理。常用参数范围如下：

实际工程中我一般会把首轮迭代的 α_p 设在 0.8，跑到中期如果压力残差出现周期震荡，就降到 0.5。简单判断方法是打印每一步的质量源项b的绝对值，如果它开始反弹，说明 α_p 偏大。

残差统计不能只观察压力修正量。更可靠的方式是计算动量方程的无量纲残差：

res_u = sum(sum(abs(u - u_old))) / (sum(sum(abs(u))) + 1e-30); res_v = sum(sum(abs(v - v_old))) / (sum(sum(abs(v))) + 1e-30); if max(res_u, res_v) < 1e-6 break; end

注意分母加1e-30是为了防止初始速度场全零时除零，这在低雷诺数启动阶段是常见错误。

4.3 顶盖驱动方腔流的边界条件处理

方腔流是验证 NS solver 的标配算例。边界条件为：三面固壁无滑移，顶盖以恒定速度向右运动。交错网格下设置边界比普通网格更直接，因为速度节点本身就落在边界上。顶盖处v=0，u=1；左右壁面u=0、v=0。压力边界采用零梯度，即边界外层压力等于内侧压力。实际代码里，压力修正方程在边界控制体积上需要关闭垂直于边界的修正系数，否则会出现边界上压力被强行修正的伪效应。

5. 用方腔流基准验证你的解算器

5.1 初级涡心位置对比

对方腔流做验证时，最直观的指标是初级涡心位置。以 Re=100 为例，参照经典文献，涡心应位于约 x≈0.62、y≈0.73 附近。你可以计算出流函数后找极小值点：先对速度场积分得到流函数，再用min找到最小值索引，换算成坐标。偏差在 1% 以内说明网格和插值格式都是对的。如果涡心明显偏向中心，通常是数值耗散过大，把 QUICK 换成了纯迎风导致。

5.2 从残差曲线判断收敛阶段

非定常求解每一时间步内的 SIMPLE 迭代都要看残差曲线。绘制方法是把每个外迭代步的质量源项和动量残差存在数组里：

figure; semilogy(1:niter, res_p(1:niter), 'o-'); xlabel('SIMPLE 迭代步'); ylabel('质量残差（无量纲）'); grid on;

理想曲线是前 20 步快速下降，之后进入平缓下降段。如果曲线在某个水平线上来回振荡，优先检查压力松弛因子是否偏大。如果曲线先降后升，则时间步长过大，导致动量预测步发散。这是判断显式时间推进是否稳定的最快手段，比盯着速度云图直观得多。

5.3 涡量场可视化

最后一步是绘制涡量云图，用它观察顶盖附近的剪切层。涡量定义为 ω = ∂v/∂x - ∂u/∂y，在交错网格上可以直接用二阶中心差分计算：

% 计算涡量场，确保维度一致 omega = zeros(nx, ny); omega(2:nx-1, 2:ny-1) = (v(2:nx-1, 2:ny-1) - v(1:nx-2, 2:ny-1)) / dx ... - (u(2:nx-1, 2:ny-1) - u(2:nx-1, 1:ny-2)) / dy;

注意交错网格上 v 的 x 方向差分需要使用 v 数组在界面上的相邻点，这里v数组比压力数组在 y 方向多一个点，因此索引对应关系必须从 2 开始，避免取到边界外的值。压力修正算法的价值在这个环节体现得最清楚：如果压力场只是近似无散，涡量云图上会出现沿网格线分布的条状伪影；把压力修正做到残差低于 1e-6 之后，这些条状伪影才会从 ω 场中消失。

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