★ FEATURED ARTICLE 简介面向学习有限元流体计算的高年级本科生、研究生与工程技术人员这套 Matlab 程序围绕二维稳态 Navier-Stokes 方程将动量方程与连续性方程离散为代数系统覆盖网格生成、基函数与高斯积分、弱形式组装、Dirichlet 与应力边界处理、Newton 迭代及压力定点修正、误差计算等环节并配有多个解析解函数用于验证收敛阶。全部 45 个文件均为 .m 脚本与函数压缩包仅 19KB按单元计算、矩阵装配、边界处理、求解与后处理分模块拆分命名清晰便于逐个阅读、调用与改写。目前已有 315 人学习下载适合作为课程作业、科研入门或自编有限元代码的参考底稿读者可从中掌握刚度矩阵与载荷向量的组装流程、边界条件的施加方式以及 Newton 求解的组织思路稍作修改即可迁移到热传导、扩散等同类偏微分方程问题。1. 二维稳态 Navier-Stokes 方程有限元求解一份能改的 Matlab 程序比一个收敛解更值钱第一次拿到二维稳态 Navier-Stokes 方程的有限元程序时我把它当黑匣子用跑方腔算例一分钟出彩色流线图改个入口条件却卡了三天。后来才想明白Matlab 做流体计算真正的价值不在算得快而在于整条装配链路——网格编号、形函数、单元矩阵、非线性迭代——全摊在编辑器里能断点、能单步、能把中间矩阵拉出来看。这套程序干的就是这件事用 Q2-Q1 混合单元离散速度与压力Picard 或 Newton 迭代推到稳态输出速度场、压力场和残差历史。它适合两类人。一类是刚学完有限元理论、想看看弱形式到底怎么变成稀疏矩阵的学生另一类是手里有算例、需要一个能随手改边界条件和雷诺数的原型工具的开发。商用软件当然更稳但你想在动量方程里加一项自定义源项、或者换一种稳定化格式做对比改别人的黑盒不如改自己的三百行代码。2. 程序骨架与技术选型六个文件撑起一条混合有限元流水线2.1 目录结构主脚本只负责串联算法都沉到函数里我见过太多把八百行塞进一个 .m 文件的版本改一个边界条件要翻半小时找变量。这套程序按“数据流”切文件网格归网格形函数归形函数装配归装配主脚本只做参数登记和结果展示。目录大致如下。ns2d/ ├── run_cavity.m % 主入口登记参数、调求解、出图 ├── buildMesh.m % 结构化四边形网格 Q2/Q1 两套节点映射 ├── shapeQ2Q1.m % 参考单元上的形函数及其导数 ├── gaussQuad.m % Gauss-Legendre 点与权 ├── assembleSystem.m % 单元循环三对角块装配成全局矩阵 ├── applyBC.m % Dirichlet 消元 压力零均值约束 ├── newtonSolve.m % Picard / Newton 非线性迭代主循环 └── postProcess.m % 涡量、流线、中心线剖面导出对应各文件职责如下表读代码之前先建立这张映射后面看任何一个函数都不会迷路。文件输入输出关键依赖buildMesh.mLx, Ly, nx, ny节点坐标、单元连接表、自由度计数无shapeQ2Q1.m参考坐标 (xi, eta)9 个速度形函数 4 个压力形函数及导数无assembleSystem.mmesh、当前解、nu全局稀疏矩阵 A 与右端项 bshapeQ2Q1、gaussQuadapplyBC.mA、b、边界自由度与值施加边界后的 A、b无newtonSolve.mmesh、nu、边界结构收敛解与残差历史上面全部postProcess.m速度、压力、网格涡量场、剖面数据、图窗无这种切法的好处是每一层都能单独验证。你怀疑形函数写错了写个xi0, eta0的调用打印一下 N 是否等于 1/4 那一组数就行不用跑整个求解器。我一般会先这么干一遍再往下走省很多事后排查。2.2 单元选型为什么速度用二次、压力用一次速度、压力同阶插值是最容易翻车的做法。不可压约束 ∇·u0 相当于给速度场加了一个逐点约束如果压力和速度用同一套节点离散后的鞍点系统不满足 inf-sup 条件压力会冒出周期为两倍网格步长的棋盘格振荡——图看起来挺花哨实际上完全是数值伪影。新手最常见的疑问就是“为什么我的压力场一格一格的”答案基本都在这里。Q2-Q1 配对二次速度、线性压力有些教材按提出者叫它 Taylor-Hood 元是工程上最稳妥的选择速度 9 节点、压力 4 节点压力自由度嵌在速度自由度里满足 LBB 稳定性条件不需要额外加稳定化项。代价是每个单元的自由度从 8 涨到 22矩阵规模大概翻三倍但对二维问题这个量级完全扛得住。另一个选型是离散方式。结构化四边形网格配等参变换雅可比在矩形单元上是常数省掉一半计算三角网格更贴复杂边界但 Q2-Q1 在三角元上的对应版本P2-P1形函数推导更长。这套程序选四边形是因为它面向的第一个算例就是规则区域先把链路跑通再谈复杂几何。提示如果你的算例有曲面边界圆柱、翼型四边形网格的贴体性会变差这时候要么换三角元要么做边界层加密不要指望粗网格硬扛。2.3 自由度排布速度与压力共享一个全局向量混合有限元最后要解的是一个鞍点系统速度未知量和压力未知量必须放进同一个列向量否则没法用一次A\b搞定。常见排布是“速度在前、压力在后”且每个速度节点占两个槽位x 分量、y 分量。% ndofU速度标量自由度总数 2 * 速度节点数 % ndof 全局自由度总数 ndofU 压力节点数 ndofU 2 * mesh.nu; ndof ndofU mesh.np; % 给定单元 e 的 9 个速度节点编号 eU9×1与 4 个压力节点编号 eP4×1 % 展开成 22 个全局自由度位置先 x y 交替的 18 个再 4 个压力 idxU reshape([2*mesh.elemU(e,:) - 1; 2*mesh.elemU(e,:)], [], 1); % 18×1 idxP ndofU mesh.elemP(e,:); % 4×1 loc [idxU; idxP]; % 22×12*node-1和2*node分别拿到该节点的 x、y 分量位置reshape之后顺序是[u1x; u1y; u2x; u2y; ...]正好和后面局部矩阵的行列顺序对上。压力从ndofU1开始编号保证整块矩阵的右下方是连续性方程那几行。这个排布不是唯一的但一旦定下来就必须全程一致局部矩阵、边界条件列表、初始解、后处理的取数下标都要按同一套来。我在两个项目之间搬代码时踩过这个坑一边按“所有 x 分量连续、所有 y 分量连续”排一边按“按节点交替”排结果速度场前一半是 u 后一半是 v残差收敛得漂亮图画出来像被打乱的拼图。3. 网格与形函数从节点坐标到单元矩阵的完整链路3.1 结构化四边形网格与两套节点编号网格生成看起来是体力活实际上决定了后面所有索引的正确性。Q2 速度节点要落在单元的角点、边中点和中心一共 9 个Q1 压力节点只落角点4 个。最省事的做法是先生成一张“细网格”2nx1 乘 2ny1速度节点全部落在这张细网格上压力节点取细网格中行列都为奇数的那些点天然就是 Q1 的角点。function mesh buildMesh(Lx, Ly, nx, ny) % 生成结构化四边形网格nx×ny 个单元 % 速度节点细网格 (2nx1)×(2ny1)压力节点(nx1)×(ny1) mesh.nx nx; mesh.ny ny; nrx 2*nx 1; nry 2*ny 1; % 细网格行列数 xv linspace(0, Lx, nrx); yv linspace(0, Ly, nry); [Xg, Yg] meshgrid(xv, yv); % 注意 meshgrid 的 Xg 是 nry×nrx mesh.X Xg; mesh.Y Yg; mesh.nu nrx * nry; % 速度标量自由度的节点数 mesh.np (nx 1) * (ny 1); % 压力节点数 mesh.elemU zeros(nx*ny, 9); % 每个单元 9 个速度节点 mesh.elemP zeros(nx*ny, 4); % 每个单元 4 个压力节点 c 0; for j 1:ny for i 1:nx c c 1; r0 2*j - 1; c0 2*i - 1; % 单元左下角在细网格上的位置 rU [r0 r0 r0 r01 r01 r01 r02 r02 r02]; cU [c0 c01 c02 c0 c01 c02 c0 c01 c02]; mesh.elemU(c,:) rU (cU - 1) * nry; % 线性编号 行 (列-1)*行数 % 压力角点左下、右下、右上、左上逆时针 mesh.elemP(c,:) [ j (i-1)*(ny1), ... j i *(ny1), ... j 1 i *(ny1), ... j 1 (i-1)*(ny1) ]; end end endelemU的排列顺序是“先沿 xi 方向、再沿 eta 方向”的 3×3 张量积扫描这一点必须和形函数文件里的循环顺序严格一致否则单元矩阵会被装到错误的节点对上。这是整套代码里最容易出错、又最难从结果上看出来的一处我在 3.2 会再强调一遍。elemP的四个角点按逆时针排列对应参考单元上 (-1,-1)、(1,-1)、(1,1)、(-1,1)。逆时针的好处是雅可比行列式恒为正不会出现单元翻转导致的符号错误。3.2 Q2-Q1 形函数一次写对后面全靠它参考单元取 [-1,1]²。一维二次 Lagrange 基函数在 -1、0、1 三个点上的表达式是现成的张量积之后就是 9 个 Q2 形函数Q1 就是四个双线性函数。function [N, dN, Nq] shapeQ2Q1(xi, eta) % 参考单元上的 Q29 个速度形函数与 Q14 个压力形函数 % 输出 % N : 9×1 速度形函数值 % dN : 9×2 速度形函数对 (xi, eta) 的偏导 % Nq : 4×1 压力形函数值 L [0.5*xi*(xi-1); 1-xi^2; 0.5*xi*(xi1)]; % 一维二次基节点 -1/0/1 dL [xi-0.5; -2*xi; xi0.5]; % 对应导数 M [0.5*eta*(eta-1); 1-eta^2; 0.5*eta*(eta1)]; dM [eta-0.5; -2*eta; eta0.5]; N zeros(9,1); dN zeros(9,2); k 0; for b 1:3 % eta 方向 for a 1:3 % xi 方向内层变化最快必须与 elemU 的列优先顺序一致 k k 1; N(k,1) L(a) * M(b); dN(k,1) dL(a) * M(b); dN(k,2) L(a) * dM(b); end end % Q1 压力形函数顺序对应 elemP 的逆时针四角 Nq 0.25 * [ (1-xi)*(1-eta); (1xi)*(1-eta); (1xi)*(1eta); (1-xi)*(1eta) ]; enda是内层循环所以 k 从 1 到 3 时 xi 从 -1 变到 1、eta 固定在 -1这和elemU里rU [r0 r0 r0 ...]、cU [c0 c01 c02 ...]的写法完全对得上。如果哪天你改成外层 a 内层 b所有单元矩阵都会错位但程序照样能跑出数只是解会变得莫名其妙——这类问题最难查因为残差曲线看起来很正常。形函数写完后建议做两个手工检查在(xi,eta)(0,0)处9 个 N 值应该只有中心那个等于 1其余全为 0四个角点处 Q1 的 Nq 应该只有一个等于 1。两个检查都过这套基函数就可以放心用了。3.3 等参变换与高斯积分把参考单元映射到物理单元参考单元上的积分要乘雅可比行列式搬到物理单元。四边形等参变换写出来就是J dN * [x_node, y_node]其中dN是 9×2 的参考坐标导数[x_node, y_node]是 9×2 的物理坐标。function [gp, gw] gaussQuad(n) % n 点 Gauss-Legendre 积分返回 n^2 个二维点与权张量积 [x, w] gaussLegendre1D(n); gp zeros(n*n, 2); gw zeros(n*n, 1); k 0; for j 1:n for i 1:n k k 1; gp(k,:) [x(i), x(j)]; gw(k) w(i) * w(j); end end end function [x, w] gaussLegendre1D(n) % 1~3 点 Gauss-Legendre够 Q2 单元用 switch n case 2 x [-1; 1] / sqrt(3); w [1; 1]; case 3 x [-sqrt(3/5); 0; sqrt(3/5)]; w [5/9; 8/9; 5/9]; otherwise error(只实现了 2 点和 3 点); end endQ2 形函数是二次的被积函数里出现形函数导数的乘积是二次、三个形函数相乘是六次严格说 3×3 高斯点不够精确。但实际算下来 3×3 对 Q2-Q1 的精度已经足够这是有限元教材里的经典结论而 4×4 会让装配时间翻倍性价比不高。我一般用 3×3只有在做收敛阶验证看误差随网格按 h³ 下降时才换成 4×4 排除积分误差的干扰。物理坐标下的形函数导数用dNdx dN * inv(J)得到注意维度dN是 9×2J是 2×2乘出来仍是 9×2第 1 列是 ∂N/∂x第 2 列是 ∂N/∂y。这一点和很多教材里的转置写法不同抄公式的时候一定要对着自己的维度检查一遍。注意矩形单元上inv(J)每次高斯点都一样可以提到循环外算一次省时间一旦换成非结构网格或者曲线边界雅可比随积分点变化就必须规规矩矩逐点求逆。4. 非线性方程组装配Picard 与 Newton 到底选哪条路4.1 弱形式逐项对应到局部矩阵稳态不可压 Navier-Stokes 写成弱形式找 (u, p) 使对任意测试函数 (v, q) 有粘性项∫ ν ∇u : ∇v dΩ对流项∫ (u·∇u)·v dΩ压力项-∫ p (∇·v) dΩ连续性∫ q (∇·u) dΩ 0这四项在代码里就是局部 22×22 矩阵的四个块粘性和对流进 (u,u) 块压力项进 (u,p) 块连续性进 (p,u) 块(p,p) 块是零——这正是鞍点系统的来源。局部自由度顺序是[u1x u1y ... u9x u9y p1 p2 p3 p4]前 18 个是速度后 4 个是压力。4.2 Picard 迭代稳但收敛慢的那条路Picard 的做法是把对流项里的对流速度“冻结”成上一步的解 u^k只对 u^{k1} 线性化。好处是线性系统始终对称性较好、收敛域宽坏处是收敛速度只有线性高雷诺数下要迭代几十次。function [A, b] assembleSystem(mesh, U, P, nu, f) % 组装 Picard 线性化后的全局矩阵A * [dU; dP] b % U: nu×2 当前速度P: np×1 当前压力f: 体力项可为空 ndofU 2*mesh.nu; ndof ndofU mesh.np; [gp, gw] gaussQuad(3); rows []; cols []; vals []; rhs zeros(ndof,1); % 三元组方式装配 for e 1:size(mesh.elemU,1) eU mesh.elemU(e,:); eP mesh.elemP(e,:); xy [mesh.X(eU), mesh.Y(eU)]; % 9×2 物理坐标 ue U(eU, :); % 9×2 单元速度 pe P(eP); % 4×1 单元压力 Ae zeros(22,22); for g 1:size(gp,1) [N, dN, Nq] shapeQ2Q1(gp(g,1), gp(g,2)); J dN * xy; % 2×2 雅可比 dJ det(J); if dJ 0, error(单元 %d 出现负雅可比检查节点顺序, e); end dNdx dN / J; % 9×2 物理导数 w gw(g) * dJ; ugp N * ue; % 高斯点处速度 1×2用于冻结对流 for a 1:9 for bIdx 1:9 % 粘性项 ν ∇N_a : ∇N_b kv nu * (dNdx(a,1)*dNdx(bIdx,1) dNdx(a,2)*dNdx(bIdx,2)); % 对流项 (u_gp · ∇N_b) N_au_gp 用上一步解冻结 cv (ugp(1)*dNdx(bIdx,1) ugp(2)*dNdx(bIdx,2)) * N(a); val (kv cv) * w; Ae(2*a-1, 2*bIdx-1) Ae(2*a-1, 2*bIdx-1) val; % u-u Ae(2*a, 2*bIdx ) Ae(2*a, 2*bIdx ) val; % v-v end for c 1:4 % 压力梯度项 -p (∇·v) Ae(2*a-1, 18c) Ae(2*a-1, 18c) - dNdx(a,1) * Nq(c) * w; Ae(2*a, 18c) Ae(2*a, 18c) - dNdx(a,2) * Nq(c) * w; % 连续性 ∫ q (∇·u) Ae(18c, 2*a-1) Ae(18c, 2*a-1) Nq(c) * dNdx(a,1) * w; Ae(18c, 2*a ) Ae(18c, 2*a ) Nq(c) * dNdx(a,2) * w; end end end % 局部自由度映射到全局 idxU reshape([2*eU-1; 2*eU], [], 1); loc [idxU; ndofU eP]; [ii, jj, vv] find(Ae); rows [rows; loc(ii)]; cols [cols; loc(jj)]; vals [vals; vv]; end A sparse(rows, cols, vals, ndof, ndof); % 重复三元组会被自动求和 b rhs; end几个参数值得单独说。nu是运动粘度和雷诺数的关系是nu U_ref * L_ref / Re方腔算例里通常取U_ref1、L_ref1所以nu 1/Re。ugp是从上一步解插值出来的高斯点速度它在这一轮里是常数不参与求导——这就是 Picard 和 Newton 的唯一区别。find(Ae)只取非零元局部 22×22 里绝大多数位置是零这样能省一半三元组。sparse会自动把重复的 (i,j) 累加所以单元循环里不需要手动做加法这个特性是装配代码能写得这么短的关键。4.3 Newton 修正二次收敛的代价是 JacobianNewton 是在残差 R(U,P)0 上做的需要在 Picard 那三项之外再补一项对流项的线性化多出来的∫ (δu·∇u^k)·v。也就是说 (u,u) 块里除了粘性和冻结对流还要加一块非对称贡献。% 在 4.2 的双重循环内与 cv 并列追加 % Newton 附加项(δu · ∇u_gp) · N_a其中 δu 是被求导的那一支 cn N(bIdx) * (dNdx(a,1)*ugp(1) dNdx(a,2)*ugp(2)); val (kv cv cn) * w;三个字符的差别换来的是从线性收敛到二次收敛。方腔算例 Re100 时Picard 大概 15 到 25 步收敛Newton 一般 5 到 8 步就压到 1e-10但 Newton 对初值敏感得多如果一开始速度全是零对流项 Jacobian 几乎为零第一步等于在解 Stokes问题不大如果初值是上一组错参数留下来的乱场Newton 很容易一步跳到负压力发散。我自己的习惯是“混合策略”前 5 步用 Picard 把解拉到收敛域附近后面切 Newton 收尾。实现上就是一个 if 判断if k 5, cn 0; end几行代码的事但能省掉大量调参时间。这个血泪经验是 Re1000 的算例教会我的纯 Newton 从零初值出发几乎必炸。4.4 边界条件、压力零均值与收敛判据方腔算例的边界很直观上边界 u(1,0)其余三边 u(0,0)。但压力没有边界条件——纯 Neumann 的压力场有一个常数零空间不解掉的话线性系统是奇异的Matlab 会给你一个警告然后返回一个压力整体漂了几千的解。消除零空间的常见做法有两种。一是锚定随便挑一个压力自由度把对应行列设为单位行右端设为 0等价于强制 p0。二是加约束行把某个方程换成∫p dΩ 0。前者简单但会让矩阵失去对称性后者更“正规”但要多写几十行。做原型我选第一种只要能收敛压力相对分布是对的。function [A, b] applyBC(A, b, dof, val, pAnchor) % dof: 边界自由度索引列向量val: 对应值 % pAnchor: 用来锚定压力的那个自由度编号 b b - A(:, dof) * val; % 先把已知量搬到右端 A(:, dof) 0; A(dof, :) 0; A(dof, dof) speye(numel(dof)); % 边界行只留对角 1 b(dof) val; % 压力锚定把 pAnchor 那一行一列清空强制 p(pAnchor) 0 A(pAnchor, :) 0; A(:, pAnchor) 0; A(pAnchor, pAnchor) 1; b(pAnchor) 0; endb b - A(:,dof)*val这一行的顺序不能反。先把边界自由度对应的列乘上已知值减到右端再去清零行列——反过来做右端项就丢掉了边界贡献解出来的速度在壁面附近会有一个明显的错位。这是新手最常犯的一个顺序错误而且程序不报错只是解不对。收敛判据用相对增量||Δx|| / ||x||比较稳阈值取 1e-8 到 1e-10。只看||Δx||的绝对值在高雷诺数下会过早停止解的量级本身很大时绝对增量看着小其实离收敛还远。5. 常见问题排查五个让残差停在 1e-2 的坑5.1 线性系统层面的三个坑现象一A\b报“矩阵接近奇异”或者压力解整体漂到 ±1e4 量级。原因几乎一定是压力零空间没消掉——要么忘了锚定要么applyBC里锚定的那一行又被后面的循环覆盖了。检查方法很简单算一下rank(full(A))和size(A,1)差 1 就是零空间还在。解决就是在applyBC末尾加锚定并且确认它被执行在最后一次修改 A 之后。现象二残差第一步就到 1e-13解全零。这种情况通常是边界条件把驱动力一起清掉了。比如顶盖速度边界设在全局自由度 1 到 2n 上结果和某个消元步骤撞车右端项被清零。排查时打印nnz(b)如果边界施加后右端项里非零元数量比预期少很多就是这里出了问题。另一个可能是雷诺数设成了Re nu而不是nu 1/Re粘度大到 100流动直接被粘死。现象三网格加密后解反而变差壁面附近速度剖面明显偏低。八成是dNdx dN / J里dN的维度搞反了或者J用了转置。做一次 patch test 就能定位给一个线性速度场 ux、v-y散度为零且满足动量方程看离散残差是不是机器精度量级。如果不是问题一定在形函数导数或雅可比。5.2 迭代与后处理层面的两个坑现象四Re100 能收敛调到 Re400 残差开始来回震荡甚至发散。Picard 的收敛域随雷诺数上升而收窄这是方法本身的限制不是代码 bug。解决办法是参数延拓先用 Re100 算到收敛把解存下来当 Re200 的初值逐步往上加每次加 100。这套程序里就是一个外层循环改十几行。加完之后 Re1000 也能稳稳推到稳态比死磕初值省事得多。现象五压力场出现周期为两倍网格步长的棋盘格。要么是单元配对不满足 LBB比如你手滑把shapeQ2Q1里的 Q1 换成了双线性速度要么是节点顺序错位。后者更隐蔽elemU按“列优先”排形函数按“行优先”排两者错一位压力就会出现规律的振荡而速度场看起来居然还行。排查方式是把elemU第一行和shapeQ2Q1在 9 个参考节点上取到的值打印出来对照一一对应才算过关。提示所有这五类问题里只有第一类会明确报错其余四类都会“正常跑完”。所以每次改完代码先跑一个已知答案的小算例比如 Re1 的方腔解接近 Stokes 解再去碰大算例。6. 二次开发换算例、换单元、换求解器的三条改造路线这套代码的可改性主要来自三层分离几何在buildMesh离散在shapeQ2Q1物理在assembleSystem。想改什么就只动对应的那一层。改造目标需要动的文件工作量注意事项换几何后台阶、管道扩张buildMesh 边界自由度列表小非矩形区域要重写单元连接逻辑雅可比不再是常数加体力项/源项assembleSystem的rhs累加很小记得乘N(a)*w做积分换单元Q1-Q1 加稳定化shapeQ2Q1 局部矩阵维度中必须补压力稳定化项否则棋盘格立刻出现换求解器GMRES ILUnewtonSolve里一行本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
