Navier-Stokes方程有限元Matlab程序：Q2-Q1混合单元与迭代求解

简介：面向学习有限元流体计算的高年级本科生、研究生与工程技术人员，这套 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 % 涡量、流线、中心线剖面导出

对应各文件职责如下表，读代码之前先建立这张映射，后面看任何一个函数都不会迷路。

这种切法的好处是每一层都能单独验证。你怀疑形函数写错了，写个xi=0, eta=0的调用打印一下 N 是否等于 1/4 那一组数就行，不用跑整个求解器。我一般会先这么干一遍再往下走，省很多事后排查。

2.2 单元选型：为什么速度用二次、压力用一次

速度、压力同阶插值是最容易翻车的做法。不可压约束 ∇·u=0 相当于给速度场加了一个逐点约束，如果压力和速度用同一套节点，离散后的鞍点系统不满足 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 个速度节点编号 eU（9×1）与 4 个压力节点编号 eP（4×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×1

2*node-1和2*node分别拿到该节点的 x、y 分量位置，reshape之后顺序是[u1x; u1y; u2x; u2y; ...]，正好和后面局部矩阵的行列顺序对上。压力从ndofU+1开始编号，保证整块矩阵的右下方是连续性方程那几行。

这个排布不是唯一的，但一旦定下来就必须全程一致：局部矩阵、边界条件列表、初始解、后处理的取数下标都要按同一套来。我在两个项目之间搬代码时踩过这个坑，一边按“所有 x 分量连续、所有 y 分量连续”排，一边按“按节点交替”排，结果速度场前一半是 u 后一半是 v，残差收敛得漂亮，图画出来像被打乱的拼图。

3. 网格与形函数：从节点坐标到单元矩阵的完整链路

3.1 结构化四边形网格与两套节点编号

网格生成看起来是体力活，实际上决定了后面所有索引的正确性。Q2 速度节点要落在单元的角点、边中点和中心，一共 9 个；Q1 压力节点只落角点，4 个。最省事的做法是先生成一张“细网格”（2nx+1 乘 2ny+1），速度节点全部落在这张细网格上；压力节点取细网格中行列都为奇数的那些点，天然就是 Q1 的角点。

function mesh = buildMesh(Lx, Ly, nx, ny) % 生成结构化四边形网格：nx×ny 个单元 % 速度节点：细网格 (2nx+1)×(2ny+1)；压力节点：(nx+1)×(ny+1) 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 r0+1 r0+1 r0+1 r0+2 r0+2 r0+2]; cU = [c0 c0+1 c0+2 c0 c0+1 c0+2 c0 c0+1 c0+2]; mesh.elemU(c,:) = rU + (cU - 1) * nry; % 线性编号 = 行 + (列-1)*行数 % 压力角点：左下、右下、右上、左上（逆时针） mesh.elemP(c,:) = [ j + (i-1)*(ny+1), ... j + i *(ny+1), ... j + 1 + i *(ny+1), ... j + 1 + (i-1)*(ny+1) ]; end end end

elemU的排列顺序是“先沿 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) % 参考单元上的 Q2（9 个速度形函数）与 Q1（4 个压力形函数） % 输出： % N : 9×1 速度形函数值 % dN : 9×2 速度形函数对 (xi, eta) 的偏导 % Nq : 4×1 压力形函数值 L = [0.5*xi*(xi-1); 1-xi^2; 0.5*xi*(xi+1)]; % 一维二次基，节点 -1/0/1 dL = [xi-0.5; -2*xi; xi+0.5]; % 对应导数 M = [0.5*eta*(eta-1); 1-eta^2; 0.5*eta*(eta+1)]; dM = [eta-0.5; -2*eta; eta+0.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); (1+xi)*(1-eta); (1+xi)*(1+eta); (1-xi)*(1+eta) ]; end

a是内层循环，所以 k 从 1 到 3 时 xi 从 -1 变到 1、eta 固定在 -1，这和elemU里rU = [r0 r0 r0 ...]、cU = [c0 c0+1 c0+2 ...]的写法完全对得上。如果哪天你改成外层 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 end

Q2 形函数是二次的，被积函数里出现形函数导数的乘积是二次、三个形函数相乘是六次，严格说 3×3 高斯点不够精确。但实际算下来 3×3 对 Q2-Q1 的精度已经足够（这是有限元教材里的经典结论），而 4×4 会让装配时间翻倍，性价比不高。我一般用 3×3，只有在做收敛阶验证（看误差随网格按 h³ 下降）时才换成 4×4 排除积分误差的干扰。

物理坐标下的形函数导数用dNdx = dN * inv(J)得到，注意维度：dN是 9×2，J是 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^{k+1} 线性化。好处是线性系统始终对称性较好、收敛域宽；坏处是收敛速度只有线性，高雷诺数下要迭代几十次。

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_a，u_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, 18+c) = Ae(2*a-1, 18+c) - dNdx(a,1) * Nq(c) * w; Ae(2*a, 18+c) = Ae(2*a, 18+c) - dNdx(a,2) * Nq(c) * w; % 连续性 ∫ q (∇·u) Ae(18+c, 2*a-1) = Ae(18+c, 2*a-1) + Nq(c) * dNdx(a,1) * w; Ae(18+c, 2*a ) = Ae(18+c, 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_ref=1、L_ref=1，所以nu = 1/Re。ugp是从上一步解插值出来的高斯点速度，它在这一轮里是常数，不参与求导——这就是 Picard 和 Newton 的唯一区别。find(Ae)只取非零元，局部 22×22 里绝大多数位置是零，这样能省一半三元组。sparse会自动把重复的 (i,j) 累加，所以单元循环里不需要手动做加法，这个特性是装配代码能写得这么短的关键。

4.3 Newton 修正：二次收敛的代价是 Jacobian

Newton 是在残差 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;

三个字符的差别，换来的是从线性收敛到二次收敛。方腔算例 Re=100 时，Picard 大概 15 到 25 步收敛，Newton 一般 5 到 8 步就压到 1e-10，但 Newton 对初值敏感得多：如果一开始速度全是零，对流项 Jacobian 几乎为零，第一步等于在解 Stokes，问题不大；如果初值是上一组错参数留下来的乱场，Newton 很容易一步跳到负压力发散。

我自己的习惯是“混合策略”：前 5 步用 Picard 把解拉到收敛域附近，后面切 Newton 收尾。实现上就是一个 if 判断，if k <= 5, cn = 0; end，几行代码的事，但能省掉大量调参时间。这个血泪经验是 Re=1000 的算例教会我的，纯 Newton 从零初值出发几乎必炸。

4.4 边界条件、压力零均值与收敛判据

方腔算例的边界很直观：上边界 u=(1,0)，其余三边 u=(0,0)。但压力没有边界条件——纯 Neumann 的压力场有一个常数零空间，不解掉的话线性系统是奇异的，Matlab 会给你一个警告然后返回一个压力整体漂了几千的解。

消除零空间的常见做法有两种。一是锚定：随便挑一个压力自由度，把对应行列设为单位行，右端设为 0，等价于强制 p=0。二是加约束行：把某个方程换成∫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; end

b = 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 就能定位：给一个线性速度场 u=x、v=-y（散度为零，且满足动量方程），看离散残差是不是机器精度量级。如果不是，问题一定在形函数导数或雅可比。

5.2 迭代与后处理层面的两个坑

现象四：Re=100 能收敛，调到 Re=400 残差开始来回震荡甚至发散。Picard 的收敛域随雷诺数上升而收窄，这是方法本身的限制，不是代码 bug。解决办法是参数延拓：先用 Re=100 算到收敛，把解存下来当 Re=200 的初值，逐步往上加，每次加 100。这套程序里就是一个外层循环，改十几行。加完之后 Re=1000 也能稳稳推到稳态，比死磕初值省事得多。

现象五：压力场出现周期为两倍网格步长的棋盘格。要么是单元配对不满足 LBB（比如你手滑把shapeQ2Q1里的 Q1 换成了双线性速度），要么是节点顺序错位。后者更隐蔽：elemU按“列优先”排，形函数按“行优先”排，两者错一位，压力就会出现规律的振荡，而速度场看起来居然还行。排查方式是把elemU第一行和shapeQ2Q1在 9 个参考节点上取到的值打印出来对照，一一对应才算过关。

提示：所有这五类问题里，只有第一类会明确报错，其余四类都会“正常跑完”。所以每次改完代码，先跑一个已知答案的小算例（比如 Re=1 的方腔，解接近 Stokes 解），再去碰大算例。

6. 二次开发：换算例、换单元、换求解器的三条改造路线

这套代码的可改性主要来自三层分离：几何在buildMesh，离散在shapeQ2Q1，物理在assembleSystem。想改什么，就只动对应的那一层。

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