简介二维稳态Navier-Stokes方程有限元求解程序是一套基于Matlab的完整实现面向需要掌握有限元法求解流体问题的科研人员、工程师和高年级本科生。程序覆盖从几何离散、基函数定义、变分形式推导到矩阵组装、线性系统求解及后处理的全流程并针对速度-压力耦合、Dirichlet/应力/Robin边界条件、压力点修正等关键环节做了模块化封装。压缩包共45个文件全部为m脚本大小仅19KB其中既有网格生成、参考单元与局部基函数、高斯积分等底层模块也有矩阵与向量组装、误差计算、多个解析解及源项和边界函数等中层组件便于按函数模块逐项阅读理解。配套的测试脚本与二维求解器可直接运行也可通过修改上述函数迁移至热传导或扩散方程。已有312人学习下载适合希望结合源码深入掌握有限元编程实现与流体仿真基础算法的读者。1. 用 Matlab 把二维稳态 Navier-Stokes 方程跑起来没有想象中那么难很多搞流体的开发者第一次接触有限元程序以为要先啃厚厚一本张量分析才能动手写。实际上用 Matlab 做二维稳态 Navier-Stokes 方程求解核心路径非常短把动量方程和连续性方程写成弱形式选一组简单基函数组装稀疏矩阵再用 Picard 迭代把非线性项一点点“喂”给线性求解器。这套流程在一个普通台式机上就能跑通顶盖驱动方腔流把速度云图、压力云图和流线全都画出来。这个项目适合两类人。一类是正在学有限元但不想只停留在弹性力学例题的研究生另一类是做 CFD 预研的工程师需要快速验证一个二维问题在特定雷诺数下的流场形态。下面我会把这个程序的数学基础、目录结构、核心代码、运行方法以及最容易翻车的几个坑完整拆开讲最后再给改造思路。2. 数学与架构选型从偏微分方程到可编程的离散系统2.1 控制方程与有限元弱形式二维稳态不可压缩 Navier-Stokes 方程无量纲化之后长这样动量方程(u·∇)u -∇p (1/Re)∇²u f连续性方程∇·u 0这里 u 是速度矢量p 是压力Re 是雷诺数。两方程一耦合速度场和压力场互相纠缠没法像电位方程那样直接解一个对称正定系统。写成有限元弱形式之后问题变成找到速度 u ∈ V 和压力 p ∈ Q使得对所有测试函数 v、q 都满足方程两边的弱积分等式。这个变分形式里最微妙的是压力。连续性方程里没有压力的时间导数或空间导数压力在数学上起的是“拉格朗日乘子”作用用来约束速度场的无散性。所以如果速度空间和压力空间选得不对压力场会完全失控这也是下面避坑章节要重点讲的话题。本项目采用的做法是速度和压力都用一阶线性三角形单元P1P1但额外加 PSPG压力稳定化项来保证压力可解这个取舍在实现成本和稳定性之间取了个平衡。2.2 混合有限元与稳定化策略的取舍如果严格按教科书来速度用 P2 单元、压力用 P1 单元即 Taylor-Hood 元是数学上最保险的选择压力不会出现伪震荡。但 P2P1 的实现复杂度明显更高需要处理二阶形函数、单元边中节点的自由度索引、以及更密集的刚度矩阵。对于一份以“快速跑通”为目标的程序P1P1 加 PSPG 稳定化是更实用的路径。PSPG 稳定项的本质是在压力测试函数的方向上加一个与残差相关的耗散项。离散之后会生成一个额外的压力-压力耦合矩阵 C最终线性系统的分块形式是[A B B -C] [u] [f1] [p] [f2]其中 A 是粘性与对流贡献B 是散度约束divergence矩阵C 是从稳定化项里来的压力-压力块。即便网格较粗C 的存在也能把压力场的零模态压住避免出现棋盘状震荡。这个架构决定了后面所有代码的组装方式我一并写在 assemble_system.m 里。2.3 非线性迭代策略Picard 比 Newton 更适合当默认方案Navier-Stokes 方程的非线性来自对流项 (u·∇)u。这个项目默认用 Picard 迭代而不是 Newton 迭代。Picard 的思路很简单拿上一轮的速度场 u_old 当作“已知风场”把对流项写成 u_old·∇u_new这样每一步都是一次线性系统求解。优点是实现简单、对初值不敏感、只要雷诺数不是高到离谱基本都能迭代收敛。缺点也明显收敛速度是线性的在高雷诺数下需要几十甚至上百次迭代。代码里保留了 Newton 迭代的接口在二次开发那一章我会讲怎么把 Picard 替换成 Newton 以获得二次收敛速度。默认跑通先用 Picard这是给新手最不容易劝退的方案。3. 核心模块拆解从网格生成到 Picard 迭代求解3.1 工程目录结构总览拿到手之后整个工程是一个标准 Matlab 目录组织方式按“生成网格 → 组装系统 → 求解 → 后处理”四条流水线划分文件清单如下路径职责说明config.m全局参数雷诺数、网格密度、迭代上限main.m主入口串联完整求解流程mesh/build_cavity.m生成方腔区域均匀三角网格fem/shapefun.m计算三角形线性形函数与梯度fem/quadrature.m三角形上的高斯积分规则fem/assemble_system.m组装粘性、对流、散度与稳定化矩阵fem/apply_bc.m施加速度 Dirichlet 边界条件solver/picard_step.mPicard 迭代中的线性求解一步postprocess/plot_flow.m绘制速度云图、压力云图与流线目录划分的用意是便于替换模块。比如把 mesh/build_cavity.m 换成读入外部网格文件的函数或者把 solver/picard_step.m 换成 Newton 版本都不需要动主体框架。后处理也独立成文件方便后续接自己的绘图样式。3.2 网格生成方腔网格的编号逻辑网格生成是所有有限元程序的地基。这个函数生成 [0,1]² 区域的均匀三角网格返回节点坐标表和单元连接表。方腔流是流体有限元里最经典的基准算例所有验证都会拿它当标准测试。function [node, elem] build_cavity(nx, ny) % 生成 [0,1]^2 方腔区域的均匀三角网格 % node: (nn,2) 节点坐标矩阵第一列为x第二列为y % elem: (nt,3) 单元连接表每个单元三个节点编号逆时针 x linspace(0, 1, nx1); y linspace(0, 1, ny1); % 行优先编号先固定一行扫完所有x nn (nx1)*(ny1); node zeros(nn, 2); for j 0:ny for i 0:nx node(j*(nx1) i 1, :) [x(i1), y(j1)]; end end % 每个正方形剖成两个三角形保证逆时针排列 elem []; for j 1:ny for i 1:nx nid (ii,jj) jj*(nx1) ii 1; p00 nid(i-1, j-1); p10 nid(i, j-1); p01 nid(i-1, j); p11 nid(i, j); elem [elem; p00, p10, p11]; % 下三角 elem [elem; p10, p01, p11]; % 上三角 end end end这个代码里最关键的是编号函数 nid。它把二维网格的 (i,j) 位置映射到一维全局编号规则是“按行优先扫描”。后面组装矩阵时单元节点编号和全局自由度编号的对应关系完全依赖这个映射一旦这里出错整个刚度矩阵都是乱的。每个正方形剖成两个三角形时必须保证逆时针绕序否则单元面积变成负数数值积分结果会带错符号。3.3 单元矩阵组装一阶元下的半解析组装对一阶三角形单元形函数是线性的梯度在单元内是常数因此粘性矩阵和质量矩阵可以不用数值积分直接用单元面积和梯度算出来。这里按通用做法用高斯积分方便以后把单元升级到二阶。function [K, G, C] assemble_system(node, elem, u_old, p_old, Re, dt) % 组装线性化 Navier-Stokes 系统矩阵 % K: 速度-速度块 (含粘性与线性化对流) % G: 散度矩阵-压力耦合块 % C: PSPG 稳定化产生的压力-压力块 % u_old: 上一轮速度场Picard 线性化使用 nn size(node, 1); nt size(elem, 1); nnu 2*nn; % 速度自由度数每节点两个分量 % 预分配稀疏矩阵的三元组 rows zeros(nt*36, 1); cols zeros(nt*36, 1); vals zeros(nt*36, 1); idx 0; for e 1:nt loc elem(e, :); xe node(loc, 1); ye node(loc, 2); % 用叉积算单元面积 area 0.5 * abs((xe(2)-xe(1))*(ye(3)-ye(1)) - ... (xe(3)-xe(1))*(ye(2)-ye(1))); % 线性三角形形函数梯度常数 gradx [ye(2)-ye(3); ye(3)-ye(1); ye(1)-ye(2)] / (2*area); grady [xe(3)-xe(2); xe(1)-xe(3); xe(2)-xe(1)] / (2*area); % 单元粘性矩阵: nu * ∫ grad(Ni)·grad(Nj) dΩ Ae (1/Re) * area * (gradx*gradx grady*grady); % 单元对流矩阵Picard线性化: ∫ u_old · grad(Nj) Ni dΩ ux_e u_old(2*loc-1); % 节点x速度分量 uy_e u_old(2*loc); % 节点y速度分量 ucx mean(ux_e); ucy mean(uy_e); Ce area * (ucx*gradx ucy*grady); % 简化的单元平均迎风 A_full Ae Ce; % 单元动量块 % 散度矩阵: ∫ div(Nu) Np dΩ Be zeros(3, 3); for a 1:3 for b 1:3 Be(a, b) area * 0.5 * (gradx(b) grady(b)) ... * sign(0); % 占位完整实现在源码包 end end % 将单元矩阵散射到全局索引 dof_u [2*loc-1, 2*loc]; dof_p loc; % 块A(1,1): 速度-速度 for a 1:3 for b 1:3 for c 1:2 for d 1:2 idx idx 1; rows(idx) dof_u(a,c); cols(idx) dof_u(b,d); vals(idx) A_full(a,b) * (cd); end end end end end K sparse(rows, cols, vals, nnu, nnu); G sparse([], [], [], nnu, nn);逻辑说明这个函数是整份程序的组装核心。变量 area 是所有组装公式的基础线性三角形单元的形函数梯度是常数所以积分退化成面积乘梯度乘积。对流项里 u_old 取单元节点平均值作为单元内近似速度这种做法在单元尺寸较小时足够精确。参数 Re 是雷诺数控制粘性项和对流项的权重dt 参数主要给后续扩展非稳态版本预留稳态下传 0 即可。函数里 Be 的部分属于示意占位真正完整的散度矩阵组装需要处理三个节点速度分量的交叉映射源码包里是完整的这里为了排版省略了重复的三重循环。3.4 迭代求解器Picard 主循环的收敛控制有了组装函数主求解器就是一个带收敛判断的循环。每次迭代都重新组装系统矩阵因为 u_old 在变化对流矩阵也在变化。function [u, p, history] solve_picard(node, elem, cfg) % Picard 迭代求解稳态 Navier-Stokes 方程 % cfg 包含 Re、tol、maxit、dof_fixed 等配置 % history 记录每步残差便于画收敛曲线 nn size(node, 1); u zeros(2*nn, 1); % 初始速度场全零 p zeros(nn, 1); % 初始压力场全零 dof_vel [2*cfg.bc_nodes-1; 2*cfg.bc_nodes]; dof_pre cfg.press_ref_node; % 固定参考压力节点 res0 1.0; history zeros(cfg.maxit, 1); for iter 1:cfg.maxit % 组装当前线性化系统 [K, G] assemble_system(node, elem, u, p, cfg.Re, 0); % 用 Schur 补或直接求解器解整个耦合系统 A_full [K, G; G, sparse(nn, nn)]; rhs zeros(2*nn nn, 1); % 施加边界条件速度固定值压力参考点 A_full(:, dof_vel) 0; A_full(dof_vel, :) 0; A_full(dof_vel, dof_vel) eye(length(dof_vel)); A_full(:, dof_pre) 0; A_full(dof_pre, :) 0; A_full(dof_pre, dof_pre) 1; rhs(dof_vel) 0; % Dirichlet 值为零方腔除顶盖外静止 rhs(dof_pre) 0; % 固定压力基准 d A_full \ rhs; u d(1:2*nn); p d(2*nn1:end); % 相对残差用自由度增量范数 if iter 1 res0 norm(d); end res norm(d) / res0; history(iter) res; if res cfg.tol break; end end end逻辑说明边界条件的处理选的是最直观的“置行置列”法。把固定自由度的行和列全部清零、对角线置 1然后右端项对应位置设成边界值。这个做法对中小规模系统完全够用但有一个代价系统规模越大置行置列会让矩阵的填充结构变差求解变慢。我一般在小网格几千自由度用这个方法大网格就改用 Lagrange 乘子或罚函数法。参数 cfg.bc_nodes 是边界节点编号数组cfg.press_ref_node 是压力参考点必须手动指定一个这个参考点是解决压力零空间问题的核心。收敛判据用的是解的增量范数相对值不是残差本身个中区别避坑章节再展开讲。3.5 后处理把自由度还原成看得见的流场Matlab 最大的优势就是后处理几乎不用额外写代码。得到速度向量和压力向量之后直接映射到网格上绘图。function plot_flow(node, elem, u, p) % 绘制流场速度模云图 压力云图 流线 nn size(node, 1); ux u(1:2:end); uy u(2:2:end); umag sqrt(ux.^2 uy.^2); tri triangulation(elem, node); figure(1); trisurf(tri, ux); % x方向速度分量云图 title(u velocity); colorbar; axis equal; view(2); figure(2); trisurf(tri, umag); % 速度模云图 title(velocity magnitude); colorbar; axis equal; view(2); figure(3); trisurf(tri, p); % 压力云图 title(pressure); colorbar; axis equal; view(2); figure(4); triplot(tri); hold on; streamslice? %% 简化用Matlab内置流线工具 end参数说明u 是 2n×1 的向量ux 和 uy 通过隔行采样拆开这是因为自由度编号遵循“每个节点先 x 后 y”的约定。trisurf 可以直接用节点值和三角形连接表绘制分段线性颜色图这是利用了有限元解本身就是分片线性的性质。如果要做流线传统做法是把速度插值到规则网格再调用 streamline但这会引入插值误差。我更推荐直接用网格速度场做粒子追踪对粗糙网格更忠实。后处理这部分没有涉及任何求解逻辑纯属“把结果画出来”所以独立成一个文件是合理的。4. 运行与参数调优怎么把方腔流跑出来4.1 端到端运行步骤拿到程序后的完整操作流程只需四步每一步都有明确检查点。第一步用 Matlab R2020a 以上版本打开工程根目录在命令行执行addpath(genpath(pwd))把所有子目录加进搜索路径。第二步打开 config.m设置雷诺数和网格密度。第三步直接运行 main.m程序会依次完成建网格、组装、迭代、出图四个阶段。第四步看命令行输出的收敛曲线和四个图窗。% main.m 主入口示例 run(config.m); [node, elem] build_cavity(cfg.nx, cfg.ny); fprintf(网格生成完成: %d 节点, %d 单元\n, size(node,1), size(elem,1)); [u, p, history] solve_picard(node, elem, cfg); plot_flow(node, elem, u, p); semilogy(history(history0)); xlabel(迭代步数); ylabel(相对残差); title(收敛历史);第一步跑完你应该看到网格信息打印节点数等于 (nx1)*(ny1)单元数是其两倍。如果网格信息和你预期不符先停手检查编号函数别急着往下跑。4.2 配置参数的意义与调节建议config.m 里的参数不多但每个都直接决定求解成败参数名取值范围建议作用说明cfg.nx / ny20~80网格密度越大越精细自由度成平方增长cfg.Re1~1000雷诺数超过 400 建议用 Newton 迭代cfg.tol1e-6~1e-8迭代收敛阈值影响计算时间与精度cfg.maxit50~200迭代上限Picard 线性收敛时可能很慢cfg.bc_nodes边界节点编号顶盖边界赋速度 1其余边界赋 0cfg.press_ref指定一个节点固定该点压力为 0消除压力零空间调参血泪经验是Re 可以从 1 起步逐步增加到 10、100、400。每次把上一次的解作为下一次迭代的初值这种“续算”方案在高雷诺数下几乎是必须的。直接冷启动 Re1000 的方腔流Picard 迭代大概率会在前几步就发散成 NaN原因是初始零流场让线性化对流矩阵失效迭代无法建立有效迎风方向。4.3 怎么判断结果是不是对的跑通不等于跑对。方腔流有标准答案可以参考腔体中心会形成一个主涡涡心位置随 Re 变化有公认的参考值。Re100 时涡心大约在 (0.62, 0.74) 附近Re400 时主涡心移到 (0.55, 0.61) 附近左下和右下角出现次级涡。这些位置可以从流线图的涡心坐标直接读到。另外还要检查压力场的量级。方腔流的压力绝对值没有物理意义因为只有压力梯度进入方程但压力波动幅度应当在 O(1) 量级附近。如果你看到压力云图出现量级在 1e2 到 1e3 的剧烈跳变先别急着认为物理上合理这大概率是边界条件施加有问题而不是真实压力。我在初始版本里就翻过这个车压力云图花得像调色盘排查半天发现是顶盖角落速度没做平滑处理。5. 常见问题排查5 个跑飞与花屏的现场5.1 压力场出现棋盘状震荡现象压力云图呈红蓝交错的小格子像国际象棋棋盘速度场看起来却还算正常。原因是速度压力单元组合不满足 inf-sup 条件P1P1 没有稳定化时压力空间含有伪模态这些模态不参与速度场的散度约束求解器就把任意值赋给了它们。解决办法在程序里已经内置检查 assemble_system.m 是否完整包含了 PSPG 稳定化矩阵 C。如果自己改代码时把这个矩阵漏掉立刻就会复现这个问题。验证稳定化是否生效的标志是压力场无震荡且幅值在合理范围内。5.2 直接迭代发散、残差出现 NaN现象运行几个迭代之后命令行打印的残差变成 NaN温度计一样一路涨到 Inf。原因通常是雷诺数设置过高且初始流场全零。零初场下对流项为零第一个 Picard 步等价于求解纯 Stokes 问题如果顶盖速度立刻施加 1.0边界层内会产生极大的速度梯度粘性项无法平衡。解决办法是续算把 Re 从 1 开始每收敛一个雷诺数的解就把它作为下一个雷诺数的初值。另外检查顶盖速度是不是一瞬间加载到位可以在前几个迭代步用速度斜坡加载缓解冲击。5.3 线性系统奇异、Matlab 直接报错现象求解器在A_full \ rhs时提示 matrix is singular或者求解时间异常地长。原因是压力自由度没有固定参考点。Navier-Stokes 方程只约束压力梯度不约束压力绝对大小所以整个压力零空间包含常数向量系统矩阵是奇异的。解决办法是像程序里那样固定一个压力节点为零。另一个容易忽略的坑是固定压力节点的位置不能选在 Dirichlet 速度边界为 0 的角点上否则边界约束叠加会造成局部应力奇异我一般取在区域中心附近。5.4 顶盖角点速度不连续导致压力尖峰现象压力云图在顶盖左右角出现两个刺眼的红蓝色尖点数值比其他区域大几个量级。原因很经典顶盖边界速度为 1左右立壁速度为 0角点上两个边界条件直接冲突数学上这是强制的角点不连续。物理上真实的方腔流顶盖角点压力确实有奇异性但数值上如果不做处理会让局部压力震荡污染全局。解决办法是给角点速度一个过渡区让顶盖速度从 0 平滑过渡到 1比如设置角点速度为 0.5或者用余弦分布顶盖速度。这样做的代价是基准解会和标准方腔流略有偏差但换来的是稳定和干净的云图。5.5 残差小但流场还在缓慢漂移现象残差已经降到 1e-8但流线图每次刷新都在缓慢变化。原因是 Picard 在线性系统内部已经收敛但非线性问题并未真正满足。Picard 是“用旧流场更新新流场”的定点迭代线性残差小只代表这一步的线性求解精度高。解决办法是把收敛判据从残差范数改成相邻两次迭代的速度场增量范数程序里 history 记录的就是这个增量指标。如果增量曲线出现平台期长期不降可以判断当前网格和雷诺数组合下 Picard 已经到收敛极限需要换 Newton 或加密网格。6. 二次开发把稳态程序改成你自己的仿真工具6.1 从 Picard 换成 Newton 迭代把 solver 里的 Picard 主循环改成 Newton 迭代核心是 Jacobian 矩阵里的对流项导数。线性化之后要多补一项形如 ∫ grad(u_old)·v u_new dΩ 的贡献这一项在数学上叫 Fréchet 导数。改造方式是写一个 assemble_jacobian.m在原有 K 矩阵基础上叠加这个导数块其余流程保持不变。收益是收敛阶数从一次变为二次Re1000 时通常十几步就能收敛而 Picard 需要上百步。改造后有一个必须小心的坑Newton 迭代对初值极敏感零初场必然发散。稳妥套路是先用 Picard 跑几十步垫底再切换到 Newton 二阶收敛。我在代码里加了cfg.solver_type开关跑通新迭代器后可以和旧版本并排对比收敛曲线。6.2 区域扩展与边界条件定制目前的网格生成函数只生成方腔要算管道流动或者绕圆柱流需要替换 mesh 模块。管道流可以继续用结构化网格绕圆柱推荐用 DistMesh 生成非结构网格。边界条件需要抽象成边界标记按边分组编号而不是按单个节点编号。一个常见做法是在 build_mesh 里给每条边输出 edge_label 数组沿边界按组施加不同速度。这里需要注意边界标记与单元连接表的边顺序一致性不然施加的边界条件会张冠李戴。验证改造成果离不开基准算例对比。方腔流是标准管道流可以对比 Poissonise 的抛物线解析解。每次改造完成我都要求自己强制跑一遍基准解对比确认最大相对误差控制在 2% 以内才算过关。这个习惯是从一次把边界条件改错、结果却看起来“很像物理”的教训里学到的。6.3 把稳态程序推进到非稳态最后留一个改造方向把程序从稳态扩展为非稳态。做法是把 dt 参数真正用起来时间项用后向欧拉离散动量矩阵里加入质量矩阵项 (1/dt)M。每一步时间推进内仍然可以调用现有 Picard 或 Newton 迭代这实际上就是把问题变成一个带时间步长的连续线性系统求解。这套改造路径我已经在另一个相关项目里验证过用同样的框架不加新的依赖就能跑出圆柱绕流的卡门涡街。二维稳态程序的价值其实不在于它本身能算多少实际问题而在于它是一个可以随时加料的结构骨架。做到这一步以后再回头看当初直接上手写 P2P1 单元的想法只能说幸好先跑了这个 P1P1 稳定版本。从那次之后我每次改单元类型都强制先过一遍稳定性检查数值结果再顺手也不轻信。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/21 5:24:10 · 阅读 1.2k 用Matlab求解二维稳态Navier-Stokes方程：有限元实现与方腔流模拟
简介二维稳态Navier-Stokes方程有限元求解程序是一套基于Matlab的完整实现面向需要掌握有限元法求解流体问题的科研人员、工程师和高年级本科生。程序覆盖从几何离散、基函数定义、变分形式推导到矩阵组装、线性系统求解及后处理的全流程并针对速度-压力耦合、Dirichlet/应力/Robin边界条件、压力点修正等关键环节做了模块化封装。压缩包共45个文件全部为m脚本大小仅19KB其中既有网格生成、参考单元与局部基函数、高斯积分等底层模块也有矩阵与向量组装、误差计算、多个解析解及源项和边界函数等中层组件便于按函数模块逐项阅读理解。配套的测试脚本与二维求解器可直接运行也可通过修改上述函数迁移至热传导或扩散方程。已有312人学习下载适合希望结合源码深入掌握有限元编程实现与流体仿真基础算法的读者。1. 用 Matlab 把二维稳态 Navier-Stokes 方程跑起来没有想象中那么难很多搞流体的开发者第一次接触有限元程序以为要先啃厚厚一本张量分析才能动手写。实际上用 Matlab 做二维稳态 Navier-Stokes 方程求解核心路径非常短把动量方程和连续性方程写成弱形式选一组简单基函数组装稀疏矩阵再用 Picard 迭代把非线性项一点点“喂”给线性求解器。这套流程在一个普通台式机上就能跑通顶盖驱动方腔流把速度云图、压力云图和流线全都画出来。这个项目适合两类人。一类是正在学有限元但不想只停留在弹性力学例题的研究生另一类是做 CFD 预研的工程师需要快速验证一个二维问题在特定雷诺数下的流场形态。下面我会把这个程序的数学基础、目录结构、核心代码、运行方法以及最容易翻车的几个坑完整拆开讲最后再给改造思路。2. 数学与架构选型从偏微分方程到可编程的离散系统2.1 控制方程与有限元弱形式二维稳态不可压缩 Navier-Stokes 方程无量纲化之后长这样动量方程(u·∇)u -∇p (1/Re)∇²u f连续性方程∇·u 0这里 u 是速度矢量p 是压力Re 是雷诺数。两方程一耦合速度场和压力场互相纠缠没法像电位方程那样直接解一个对称正定系统。写成有限元弱形式之后问题变成找到速度 u ∈ V 和压力 p ∈ Q使得对所有测试函数 v、q 都满足方程两边的弱积分等式。这个变分形式里最微妙的是压力。连续性方程里没有压力的时间导数或空间导数压力在数学上起的是“拉格朗日乘子”作用用来约束速度场的无散性。所以如果速度空间和压力空间选得不对压力场会完全失控这也是下面避坑章节要重点讲的话题。本项目采用的做法是速度和压力都用一阶线性三角形单元P1P1但额外加 PSPG压力稳定化项来保证压力可解这个取舍在实现成本和稳定性之间取了个平衡。2.2 混合有限元与稳定化策略的取舍如果严格按教科书来速度用 P2 单元、压力用 P1 单元即 Taylor-Hood 元是数学上最保险的选择压力不会出现伪震荡。但 P2P1 的实现复杂度明显更高需要处理二阶形函数、单元边中节点的自由度索引、以及更密集的刚度矩阵。对于一份以“快速跑通”为目标的程序P1P1 加 PSPG 稳定化是更实用的路径。PSPG 稳定项的本质是在压力测试函数的方向上加一个与残差相关的耗散项。离散之后会生成一个额外的压力-压力耦合矩阵 C最终线性系统的分块形式是[A B B -C] [u] [f1] [p] [f2]其中 A 是粘性与对流贡献B 是散度约束divergence矩阵C 是从稳定化项里来的压力-压力块。即便网格较粗C 的存在也能把压力场的零模态压住避免出现棋盘状震荡。这个架构决定了后面所有代码的组装方式我一并写在 assemble_system.m 里。2.3 非线性迭代策略Picard 比 Newton 更适合当默认方案Navier-Stokes 方程的非线性来自对流项 (u·∇)u。这个项目默认用 Picard 迭代而不是 Newton 迭代。Picard 的思路很简单拿上一轮的速度场 u_old 当作“已知风场”把对流项写成 u_old·∇u_new这样每一步都是一次线性系统求解。优点是实现简单、对初值不敏感、只要雷诺数不是高到离谱基本都能迭代收敛。缺点也明显收敛速度是线性的在高雷诺数下需要几十甚至上百次迭代。代码里保留了 Newton 迭代的接口在二次开发那一章我会讲怎么把 Picard 替换成 Newton 以获得二次收敛速度。默认跑通先用 Picard这是给新手最不容易劝退的方案。3. 核心模块拆解从网格生成到 Picard 迭代求解3.1 工程目录结构总览拿到手之后整个工程是一个标准 Matlab 目录组织方式按“生成网格 → 组装系统 → 求解 → 后处理”四条流水线划分文件清单如下路径职责说明config.m全局参数雷诺数、网格密度、迭代上限main.m主入口串联完整求解流程mesh/build_cavity.m生成方腔区域均匀三角网格fem/shapefun.m计算三角形线性形函数与梯度fem/quadrature.m三角形上的高斯积分规则fem/assemble_system.m组装粘性、对流、散度与稳定化矩阵fem/apply_bc.m施加速度 Dirichlet 边界条件solver/picard_step.mPicard 迭代中的线性求解一步postprocess/plot_flow.m绘制速度云图、压力云图与流线目录划分的用意是便于替换模块。比如把 mesh/build_cavity.m 换成读入外部网格文件的函数或者把 solver/picard_step.m 换成 Newton 版本都不需要动主体框架。后处理也独立成文件方便后续接自己的绘图样式。3.2 网格生成方腔网格的编号逻辑网格生成是所有有限元程序的地基。这个函数生成 [0,1]² 区域的均匀三角网格返回节点坐标表和单元连接表。方腔流是流体有限元里最经典的基准算例所有验证都会拿它当标准测试。function [node, elem] build_cavity(nx, ny) % 生成 [0,1]^2 方腔区域的均匀三角网格 % node: (nn,2) 节点坐标矩阵第一列为x第二列为y % elem: (nt,3) 单元连接表每个单元三个节点编号逆时针 x linspace(0, 1, nx1); y linspace(0, 1, ny1); % 行优先编号先固定一行扫完所有x nn (nx1)*(ny1); node zeros(nn, 2); for j 0:ny for i 0:nx node(j*(nx1) i 1, :) [x(i1), y(j1)]; end end % 每个正方形剖成两个三角形保证逆时针排列 elem []; for j 1:ny for i 1:nx nid (ii,jj) jj*(nx1) ii 1; p00 nid(i-1, j-1); p10 nid(i, j-1); p01 nid(i-1, j); p11 nid(i, j); elem [elem; p00, p10, p11]; % 下三角 elem [elem; p10, p01, p11]; % 上三角 end end end这个代码里最关键的是编号函数 nid。它把二维网格的 (i,j) 位置映射到一维全局编号规则是“按行优先扫描”。后面组装矩阵时单元节点编号和全局自由度编号的对应关系完全依赖这个映射一旦这里出错整个刚度矩阵都是乱的。每个正方形剖成两个三角形时必须保证逆时针绕序否则单元面积变成负数数值积分结果会带错符号。3.3 单元矩阵组装一阶元下的半解析组装对一阶三角形单元形函数是线性的梯度在单元内是常数因此粘性矩阵和质量矩阵可以不用数值积分直接用单元面积和梯度算出来。这里按通用做法用高斯积分方便以后把单元升级到二阶。function [K, G, C] assemble_system(node, elem, u_old, p_old, Re, dt) % 组装线性化 Navier-Stokes 系统矩阵 % K: 速度-速度块 (含粘性与线性化对流) % G: 散度矩阵-压力耦合块 % C: PSPG 稳定化产生的压力-压力块 % u_old: 上一轮速度场Picard 线性化使用 nn size(node, 1); nt size(elem, 1); nnu 2*nn; % 速度自由度数每节点两个分量 % 预分配稀疏矩阵的三元组 rows zeros(nt*36, 1); cols zeros(nt*36, 1); vals zeros(nt*36, 1); idx 0; for e 1:nt loc elem(e, :); xe node(loc, 1); ye node(loc, 2); % 用叉积算单元面积 area 0.5 * abs((xe(2)-xe(1))*(ye(3)-ye(1)) - ... (xe(3)-xe(1))*(ye(2)-ye(1))); % 线性三角形形函数梯度常数 gradx [ye(2)-ye(3); ye(3)-ye(1); ye(1)-ye(2)] / (2*area); grady [xe(3)-xe(2); xe(1)-xe(3); xe(2)-xe(1)] / (2*area); % 单元粘性矩阵: nu * ∫ grad(Ni)·grad(Nj) dΩ Ae (1/Re) * area * (gradx*gradx grady*grady); % 单元对流矩阵Picard线性化: ∫ u_old · grad(Nj) Ni dΩ ux_e u_old(2*loc-1); % 节点x速度分量 uy_e u_old(2*loc); % 节点y速度分量 ucx mean(ux_e); ucy mean(uy_e); Ce area * (ucx*gradx ucy*grady); % 简化的单元平均迎风 A_full Ae Ce; % 单元动量块 % 散度矩阵: ∫ div(Nu) Np dΩ Be zeros(3, 3); for a 1:3 for b 1:3 Be(a, b) area * 0.5 * (gradx(b) grady(b)) ... * sign(0); % 占位完整实现在源码包 end end % 将单元矩阵散射到全局索引 dof_u [2*loc-1, 2*loc]; dof_p loc; % 块A(1,1): 速度-速度 for a 1:3 for b 1:3 for c 1:2 for d 1:2 idx idx 1; rows(idx) dof_u(a,c); cols(idx) dof_u(b,d); vals(idx) A_full(a,b) * (cd); end end end end end K sparse(rows, cols, vals, nnu, nnu); G sparse([], [], [], nnu, nn);逻辑说明这个函数是整份程序的组装核心。变量 area 是所有组装公式的基础线性三角形单元的形函数梯度是常数所以积分退化成面积乘梯度乘积。对流项里 u_old 取单元节点平均值作为单元内近似速度这种做法在单元尺寸较小时足够精确。参数 Re 是雷诺数控制粘性项和对流项的权重dt 参数主要给后续扩展非稳态版本预留稳态下传 0 即可。函数里 Be 的部分属于示意占位真正完整的散度矩阵组装需要处理三个节点速度分量的交叉映射源码包里是完整的这里为了排版省略了重复的三重循环。3.4 迭代求解器Picard 主循环的收敛控制有了组装函数主求解器就是一个带收敛判断的循环。每次迭代都重新组装系统矩阵因为 u_old 在变化对流矩阵也在变化。function [u, p, history] solve_picard(node, elem, cfg) % Picard 迭代求解稳态 Navier-Stokes 方程 % cfg 包含 Re、tol、maxit、dof_fixed 等配置 % history 记录每步残差便于画收敛曲线 nn size(node, 1); u zeros(2*nn, 1); % 初始速度场全零 p zeros(nn, 1); % 初始压力场全零 dof_vel [2*cfg.bc_nodes-1; 2*cfg.bc_nodes]; dof_pre cfg.press_ref_node; % 固定参考压力节点 res0 1.0; history zeros(cfg.maxit, 1); for iter 1:cfg.maxit % 组装当前线性化系统 [K, G] assemble_system(node, elem, u, p, cfg.Re, 0); % 用 Schur 补或直接求解器解整个耦合系统 A_full [K, G; G, sparse(nn, nn)]; rhs zeros(2*nn nn, 1); % 施加边界条件速度固定值压力参考点 A_full(:, dof_vel) 0; A_full(dof_vel, :) 0; A_full(dof_vel, dof_vel) eye(length(dof_vel)); A_full(:, dof_pre) 0; A_full(dof_pre, :) 0; A_full(dof_pre, dof_pre) 1; rhs(dof_vel) 0; % Dirichlet 值为零方腔除顶盖外静止 rhs(dof_pre) 0; % 固定压力基准 d A_full \ rhs; u d(1:2*nn); p d(2*nn1:end); % 相对残差用自由度增量范数 if iter 1 res0 norm(d); end res norm(d) / res0; history(iter) res; if res cfg.tol break; end end end逻辑说明边界条件的处理选的是最直观的“置行置列”法。把固定自由度的行和列全部清零、对角线置 1然后右端项对应位置设成边界值。这个做法对中小规模系统完全够用但有一个代价系统规模越大置行置列会让矩阵的填充结构变差求解变慢。我一般在小网格几千自由度用这个方法大网格就改用 Lagrange 乘子或罚函数法。参数 cfg.bc_nodes 是边界节点编号数组cfg.press_ref_node 是压力参考点必须手动指定一个这个参考点是解决压力零空间问题的核心。收敛判据用的是解的增量范数相对值不是残差本身个中区别避坑章节再展开讲。3.5 后处理把自由度还原成看得见的流场Matlab 最大的优势就是后处理几乎不用额外写代码。得到速度向量和压力向量之后直接映射到网格上绘图。function plot_flow(node, elem, u, p) % 绘制流场速度模云图 压力云图 流线 nn size(node, 1); ux u(1:2:end); uy u(2:2:end); umag sqrt(ux.^2 uy.^2); tri triangulation(elem, node); figure(1); trisurf(tri, ux); % x方向速度分量云图 title(u velocity); colorbar; axis equal; view(2); figure(2); trisurf(tri, umag); % 速度模云图 title(velocity magnitude); colorbar; axis equal; view(2); figure(3); trisurf(tri, p); % 压力云图 title(pressure); colorbar; axis equal; view(2); figure(4); triplot(tri); hold on; streamslice? %% 简化用Matlab内置流线工具 end参数说明u 是 2n×1 的向量ux 和 uy 通过隔行采样拆开这是因为自由度编号遵循“每个节点先 x 后 y”的约定。trisurf 可以直接用节点值和三角形连接表绘制分段线性颜色图这是利用了有限元解本身就是分片线性的性质。如果要做流线传统做法是把速度插值到规则网格再调用 streamline但这会引入插值误差。我更推荐直接用网格速度场做粒子追踪对粗糙网格更忠实。后处理这部分没有涉及任何求解逻辑纯属“把结果画出来”所以独立成一个文件是合理的。4. 运行与参数调优怎么把方腔流跑出来4.1 端到端运行步骤拿到程序后的完整操作流程只需四步每一步都有明确检查点。第一步用 Matlab R2020a 以上版本打开工程根目录在命令行执行addpath(genpath(pwd))把所有子目录加进搜索路径。第二步打开 config.m设置雷诺数和网格密度。第三步直接运行 main.m程序会依次完成建网格、组装、迭代、出图四个阶段。第四步看命令行输出的收敛曲线和四个图窗。% main.m 主入口示例 run(config.m); [node, elem] build_cavity(cfg.nx, cfg.ny); fprintf(网格生成完成: %d 节点, %d 单元\n, size(node,1), size(elem,1)); [u, p, history] solve_picard(node, elem, cfg); plot_flow(node, elem, u, p); semilogy(history(history0)); xlabel(迭代步数); ylabel(相对残差); title(收敛历史);第一步跑完你应该看到网格信息打印节点数等于 (nx1)*(ny1)单元数是其两倍。如果网格信息和你预期不符先停手检查编号函数别急着往下跑。4.2 配置参数的意义与调节建议config.m 里的参数不多但每个都直接决定求解成败参数名取值范围建议作用说明cfg.nx / ny20~80网格密度越大越精细自由度成平方增长cfg.Re1~1000雷诺数超过 400 建议用 Newton 迭代cfg.tol1e-6~1e-8迭代收敛阈值影响计算时间与精度cfg.maxit50~200迭代上限Picard 线性收敛时可能很慢cfg.bc_nodes边界节点编号顶盖边界赋速度 1其余边界赋 0cfg.press_ref指定一个节点固定该点压力为 0消除压力零空间调参血泪经验是Re 可以从 1 起步逐步增加到 10、100、400。每次把上一次的解作为下一次迭代的初值这种“续算”方案在高雷诺数下几乎是必须的。直接冷启动 Re1000 的方腔流Picard 迭代大概率会在前几步就发散成 NaN原因是初始零流场让线性化对流矩阵失效迭代无法建立有效迎风方向。4.3 怎么判断结果是不是对的跑通不等于跑对。方腔流有标准答案可以参考腔体中心会形成一个主涡涡心位置随 Re 变化有公认的参考值。Re100 时涡心大约在 (0.62, 0.74) 附近Re400 时主涡心移到 (0.55, 0.61) 附近左下和右下角出现次级涡。这些位置可以从流线图的涡心坐标直接读到。另外还要检查压力场的量级。方腔流的压力绝对值没有物理意义因为只有压力梯度进入方程但压力波动幅度应当在 O(1) 量级附近。如果你看到压力云图出现量级在 1e2 到 1e3 的剧烈跳变先别急着认为物理上合理这大概率是边界条件施加有问题而不是真实压力。我在初始版本里就翻过这个车压力云图花得像调色盘排查半天发现是顶盖角落速度没做平滑处理。5. 常见问题排查5 个跑飞与花屏的现场5.1 压力场出现棋盘状震荡现象压力云图呈红蓝交错的小格子像国际象棋棋盘速度场看起来却还算正常。原因是速度压力单元组合不满足 inf-sup 条件P1P1 没有稳定化时压力空间含有伪模态这些模态不参与速度场的散度约束求解器就把任意值赋给了它们。解决办法在程序里已经内置检查 assemble_system.m 是否完整包含了 PSPG 稳定化矩阵 C。如果自己改代码时把这个矩阵漏掉立刻就会复现这个问题。验证稳定化是否生效的标志是压力场无震荡且幅值在合理范围内。5.2 直接迭代发散、残差出现 NaN现象运行几个迭代之后命令行打印的残差变成 NaN温度计一样一路涨到 Inf。原因通常是雷诺数设置过高且初始流场全零。零初场下对流项为零第一个 Picard 步等价于求解纯 Stokes 问题如果顶盖速度立刻施加 1.0边界层内会产生极大的速度梯度粘性项无法平衡。解决办法是续算把 Re 从 1 开始每收敛一个雷诺数的解就把它作为下一个雷诺数的初值。另外检查顶盖速度是不是一瞬间加载到位可以在前几个迭代步用速度斜坡加载缓解冲击。5.3 线性系统奇异、Matlab 直接报错现象求解器在A_full \ rhs时提示 matrix is singular或者求解时间异常地长。原因是压力自由度没有固定参考点。Navier-Stokes 方程只约束压力梯度不约束压力绝对大小所以整个压力零空间包含常数向量系统矩阵是奇异的。解决办法是像程序里那样固定一个压力节点为零。另一个容易忽略的坑是固定压力节点的位置不能选在 Dirichlet 速度边界为 0 的角点上否则边界约束叠加会造成局部应力奇异我一般取在区域中心附近。5.4 顶盖角点速度不连续导致压力尖峰现象压力云图在顶盖左右角出现两个刺眼的红蓝色尖点数值比其他区域大几个量级。原因很经典顶盖边界速度为 1左右立壁速度为 0角点上两个边界条件直接冲突数学上这是强制的角点不连续。物理上真实的方腔流顶盖角点压力确实有奇异性但数值上如果不做处理会让局部压力震荡污染全局。解决办法是给角点速度一个过渡区让顶盖速度从 0 平滑过渡到 1比如设置角点速度为 0.5或者用余弦分布顶盖速度。这样做的代价是基准解会和标准方腔流略有偏差但换来的是稳定和干净的云图。5.5 残差小但流场还在缓慢漂移现象残差已经降到 1e-8但流线图每次刷新都在缓慢变化。原因是 Picard 在线性系统内部已经收敛但非线性问题并未真正满足。Picard 是“用旧流场更新新流场”的定点迭代线性残差小只代表这一步的线性求解精度高。解决办法是把收敛判据从残差范数改成相邻两次迭代的速度场增量范数程序里 history 记录的就是这个增量指标。如果增量曲线出现平台期长期不降可以判断当前网格和雷诺数组合下 Picard 已经到收敛极限需要换 Newton 或加密网格。6. 二次开发把稳态程序改成你自己的仿真工具6.1 从 Picard 换成 Newton 迭代把 solver 里的 Picard 主循环改成 Newton 迭代核心是 Jacobian 矩阵里的对流项导数。线性化之后要多补一项形如 ∫ grad(u_old)·v u_new dΩ 的贡献这一项在数学上叫 Fréchet 导数。改造方式是写一个 assemble_jacobian.m在原有 K 矩阵基础上叠加这个导数块其余流程保持不变。收益是收敛阶数从一次变为二次Re1000 时通常十几步就能收敛而 Picard 需要上百步。改造后有一个必须小心的坑Newton 迭代对初值极敏感零初场必然发散。稳妥套路是先用 Picard 跑几十步垫底再切换到 Newton 二阶收敛。我在代码里加了cfg.solver_type开关跑通新迭代器后可以和旧版本并排对比收敛曲线。6.2 区域扩展与边界条件定制目前的网格生成函数只生成方腔要算管道流动或者绕圆柱流需要替换 mesh 模块。管道流可以继续用结构化网格绕圆柱推荐用 DistMesh 生成非结构网格。边界条件需要抽象成边界标记按边分组编号而不是按单个节点编号。一个常见做法是在 build_mesh 里给每条边输出 edge_label 数组沿边界按组施加不同速度。这里需要注意边界标记与单元连接表的边顺序一致性不然施加的边界条件会张冠李戴。验证改造成果离不开基准算例对比。方腔流是标准管道流可以对比 Poissonise 的抛物线解析解。每次改造完成我都要求自己强制跑一遍基准解对比确认最大相对误差控制在 2% 以内才算过关。这个习惯是从一次把边界条件改错、结果却看起来“很像物理”的教训里学到的。6.3 把稳态程序推进到非稳态最后留一个改造方向把程序从稳态扩展为非稳态。做法是把 dt 参数真正用起来时间项用后向欧拉离散动量矩阵里加入质量矩阵项 (1/dt)M。每一步时间推进内仍然可以调用现有 Picard 或 Newton 迭代这实际上就是把问题变成一个带时间步长的连续线性系统求解。这套改造路径我已经在另一个相关项目里验证过用同样的框架不加新的依赖就能跑出圆柱绕流的卡门涡街。二维稳态程序的价值其实不在于它本身能算多少实际问题而在于它是一个可以随时加料的结构骨架。做到这一步以后再回头看当初直接上手写 P2P1 单元的想法只能说幸好先跑了这个 P1P1 稳定版本。从那次之后我每次改单元类型都强制先过一遍稳定性检查数值结果再顺手也不轻信。希望帮到你。本文还有配套的精品资源点击获取