Stokes方程MATLAB二维有限元程序：从鞍点问题到稳定格式

简介：面向使用Matlab开展有限元数值计算的科研与工程人员，二维Stokes问题的完整求解程序以压缩包形式提供。Stokes方程描述忽略惯性后的不可压缩粘性流动，程序覆盖网格离散、形函数选取、变分形式转换、全局矩阵组装、线性系统求解及速度与压力后处理可视化等关键环节，适合初学者对照学习有限元实现流程。资源共135个文件，以133个m脚本为主，另有2个txt说明，压缩包整体仅93KB；其中边界处理函数单独封装，便于替换或扩展不同边界条件。已有1557人浏览学习，代码模块划分清晰，便于按步骤阅读与调试；多个核心步骤拆分为独立函数，便于二次开发与教学演示。通过运行和修改主程序，可掌握四边形或三角形网格下二维Stokes问题的建模技巧，并借助quiver、contourf等指令直接观察流场与压力分布结果。

1. stokes（matlab）二维有限元程序：先解决鞍点问题再考虑写代码

很多人拿到“stokes（matlab）二维有限元程序”这类需求时，第一反应是找现成代码改改标记，结果要么压力场走出来像棋盘，要么MATLAB直接提示矩阵奇异。Stokes方程是定常不可压N-S方程去掉对流项后留下的线性骨架，速度和压力在同一个线性系统里耦合，属于典型的鞍点问题；它能不能稳定跑通，取决于混合有限元里对LBB条件的理解，而不是盲目堆网格密度。本文按“选格式→建网格→组装→排错→验证”的顺序，把一类能落地的二维Stokes有限元程序讲透，适合有有限元基础、想把这类程序吃透的研究生或工程师。你如果只是想在MATLAB里快速出一张流线图，那工具箱可能更快；如果你需要控制每个自由度、改格式、上复杂边界条件，这篇路线能帮你少走几个月弯路。

2. 混合有限元格式怎么选：LBB 条件、Q2-P1 与 Q2-Q1 的取舍

2.1 Stokes 变分形式：速度、压力为什么拆不开

二维定常Stokes方程在不可压低速流动里常用这一组形式：

-μ Δu + ∇p = f
∇·u = 0

u是速度矢量，p是压力，μ是动力粘度，f是体积力。和热传导或弹性力学不同，这里有一个散度自由约束∇·u=0，它和动量方程耦合在同一个变分问题里。最自然的做法是把速度和压力同时当作未知数求解，这称为混合有限元。很多人习惯按弹性力学的思路，先解速度再单独求压力，放到Stokes问题里会得出完全错误的结果，因为速度空间任意“散度自由”的约束没有被强制，压力就会变成任意函数。

具体推导：对动量方程乘速度测试函数v，并在域上积分；对连续方程乘压力测试函数q并积分。于是得到弱形式：

a(u,v) + b(v,p) = (f,v)
b(u,q) = 0

这里a(u,v) = ∫ μ∇u:∇v dΩ，b(v,q) = ∫ q ∇·v dΩ，符号约定因教材而异。写成矩阵形式后，离散系统会呈现出鞍点结构：

[Kuu Kup] [U] [F]
[Kpu 0 ] [P] = [0]

其中Kuu是速度刚度矩阵，Kup和Kpu是速度-压力耦合块，右下是零块。零块的存在意味着整个系统不是对称正定，不能像弹性力学那样直接用共轭梯度法。如果程序跑出来的矩阵右下不是零块，那大概率你正在解的已经不是Stokes而是某个被改写的正则化问题。

这个变分形式还说明一个关键点：压力没有任何导数约束，只在弱形式里通过散度与速度耦合，所以压力空间不能选得比速度空间“稀疏”太远，否则系统会丧失稳定性。这正是下一小节要讲的LBB条件。

2.2 LBB 条件与形函数选型：为什么 Q2-P1 是默认牌型

LBB条件全称叫Ladyzhenskaya-Babuška-Brezzi条件，也有资料叫inf-sup条件。它的直观含义是：对任意可取的压力模式q，都存在一个速度模式v能“检验”出它，并且这个检验能力不依赖网格粗细。数学上写成：

inf_q sup_v [b(v,q) / (‖v‖_X ‖q‖_M)] ≥ β > 0

β是与网格尺寸无关的常数。如果这个条件不成立，压力离散空间里就会藏着一批零能模式，解出来的压力呈现出棋盘状或锯齿状振荡，而且网格越细振荡越明显。这不是求解器精度问题，而是有限元空间配对问题。

拿最常见的反例来说，Q1-Q1格式（速度双线性、压力双线性连续）不满足LBB，压力每隔一个单元正负交替，云图看起来像国际象棋棋盘。Q1-P0（速度双线性、压力分片常数）同样有此类spurious mode，只是振荡模式略有不同。这两个组合在快速原型里经常被新手选中，因为它们最容易实现，但结果几乎必然翻车。

真正稳定的常用组合有两套：

Q2-Q1是Taylor-Hood型格式的一种，速度节点取九节点四边形单元（或八节点serendipity单元），压力取四节点双线性连续节点。它实现简单，压力自由度和网格节点一一对应，云图可以直接画，出错时排查起来也容易。Q2-P1的压力在每个单元内部独立，跨单元允许跳跃，压力自由度需要单独管理，后处理时要先把不连续压力投影回连续节点再画图，否则看起来像马赛克。但Q2-P1在边界层和压力梯度较大的区域表现通常更锐利，很多教学型混合元程序库默认采用它。

在两个格式之间犹豫时，我一般会先选Q2-Q1：自由度少、后处理省心、满足LBB没有稳定性问题。只有在已经确定需要更高质量压力场且愿意处理投影细节时，才换成Q2-P1。

2.3 稳定化选项：什么时候需要给 Q1-Q1 加“后悔药”

如果你的网格结构已经固定，比如网格生成器只输出双线性四边形节点，不想为Q2-Q1重建一套中型节点，那么可以在Q1-Q1上做压力稳定化，让这个本不稳定的格式“后悔药”能救回来。常见的思路是PSPG/GLS稳定化：在压力方程里加入一个与压力梯度相关的耗散项，通常写成：

∑∫Ωe τ (∇p · ∇q) dΩ

这个稳定项让压力振荡被打压住，代价是原始Stokes方程被轻微修改，不再严格满足散度自由；但解会随网格细化一致收敛回原问题。稳定参数τ的选取是这里唯一的玄学：取小了没效果，取大了速度场会被过度扩散。常见做法是令τ = α h²/(4μ)，其中h是单元特征尺寸，α取0.1到1之间。我习惯先跑一次无稳定化的求解，观察压力棋盘振幅，再把α从1e-3开始逐步放大，每次只翻一倍，直到振荡消失。把这个“试参数”过程自动化成脚本，比在不动点迭代里手调快得多。

另一种稳定格式是mini-element：速度用线性三角形加中心泡函数，压力用连续线性三角形，满足LBB。但它只适合三角形网格，矩形结构网格上反而增加实现成本。决策优先级在我这里很明确：优先Q2-Q1，其次Q2-P1，只有网格结构无法重建时才考虑稳定化Q1-Q1。格式选对了，后面所有校验步骤才有意义。

3. 网格与形函数：Q2 单元的节点编号、雅可比和高斯积分

3.1 矩形网格与全局节点编号：让带宽尽量窄

为了把“物理正确”和“代码正确”分开验证，第一步要限定在矩形域[0, Lx]×[0, Ly]上，分别切nx段和ny段，一共nx*ny个单元。先不要碰四边形网格和曲边界，等矩形网上的解析解校验通过后，再改成等参元推广也不迟。Q2单元本身是任意四边形等参元实现的，但矩形网格让雅可比矩阵J变成常数对角阵，一旦程序出错，很快能判断是形函数问题还是网格几何问题。

矩形网格下，Q2速度节点总数等于细分后的节点格总数，按每边取二等分来算，就是(2nx+1)(2ny+1)，因为每条边被切成两段，中间多出一个边中点，单元内部还多一个中心节点。压力节点如果是Q1连续，总数是(nx+1)(ny+1)。如果压力量取Q2-P1不连续格式，压力自由度数是nxny3，每个单元3个三角形不连续线性自由度。两类格式的自由度管理方式差别很大，我推荐新手先把Q2-Q1跑通。

全局编号顺序会影响矩阵带宽。常见做法是按“先x方向后y方向”逐层扫描，速度节点编号结束后再编压力节点，总体自由度排序为先ux后uy再p。这种顺序下的分块结构非常清晰，spy(K)能直接看到左上角速度刚度、右上/左下耦合块、右下零块，很适合排查。如果你追求极致稀疏求解效率，后续可以用RCM排序重排自由度，但那会打散块结构，调试时不建议这么干。

一个实用的网格生成脚本骨架如下：

function [U_node, P_node, eU] = gen_mesh_rect(nx, ny, Lx, Ly) % 生成矩形域上 Q2-Q1 的网格数据 % U_node: 速度节点坐标，每个节点一行 [x, y] % P_node: 压力节点坐标（Q1连续） % eU: 单元速度节点连接表，每个单元9个节点 % % 速度节点按 (2nx+1) x (2ny+1) 网格排列 nxu = 2*nx; nyu = 2*ny; [xg, yg] = meshgrid(linspace(0, Lx, nxu+1), linspace(0, Ly, nyu+1)); U_node = [xg(:), yg(:)]; % 压力节点按 (nx+1) x (ny+1) 排列 [xp, yp] = meshgrid(linspace(0, Lx, nx+1), linspace(0, Ly, ny+1)); P_node = [xp(:), yp(:)]; % 单元连接表：每个单元9个速度节点 eU = zeros(nx*ny, 9); cnt = 0; for j = 1:ny for i = 1:nx cnt = cnt + 1; % 单元左下角速度节点编号 n = (j-1)*2*(nxu+1) + (i-1)*2 + 1; % 九节点顺序：4个角点 + 4个边中点 + 中心 eU(cnt,:) = [n, n+2, n+2*(nxu+1)+2, n+2*(nxu+1), ... n+1, n+2*(nxu+1)+3, n+2*(nxu+1)-1, n+2*(nxu+1)+1, ... n+(nxu+1)+1]; end end end

这段代码把九节点Q2单元明确落成坐标数组。角点和边中点的编号顺序必须和后面形函数求值的节点顺序一一对应，否则雅可比矩阵会算错。一个很常见的坑是：网格顶点和边中点编号顺序在图上看是“顺时钟”还是“逆时钟”不重要，只要形函数数组和连接表顺序一致就行；但单元节点顺序必须保证雅可比行列式detJ为正，否则积分权重会出现负面积，刚度矩阵符号全乱。顺带提一句，矩形网格下这个顺序天然满足detJ>0。

3.2 Q2 拉格朗日形函数与雅可比：在参考单元上算一切积分

九节点Q2单元参考坐标取[-1,1]×[-1,1]。它的形函数由一维二次拉格朗日插值拼出来。一维二次基函数在三个标准点{-1,0,1}上定义为：

L1(ξ) = 0.5 ξ(ξ-1)
L2(ξ) = 1 - ξ²
L3(ξ) = 0.5 ξ(ξ+1)

对二维Q2单元，每个形函数都是L_r(ξ) * L_s(η)的形式，其中r、s分别取1到3。节点在参考单元上对应的(ξ,η)位置决定你乘哪个一维基：角点对应两个端点基的组合，边中点对应“一个端点基+一个中点基”，中心节点则是L2(ξ)*L2(η)。这套写法比直接背公式更不容易记错，代码实现也简单。

物理坐标从参考坐标映射：

x = ∑ x_i N_i(ξ,η)
y = ∑ y_i N_i(ξ,η)

雅可比矩阵J是：

J = [dx/dξ dy/dξ
dx/dη dy/dη]

速度形函数关于物理坐标的梯度需要用J⁻¹变换：

dN/dx = J⁻¹₁₁ dN/dξ + J⁻¹₁₂ dN/dη
dN/dy = J⁻¹₂₁ dN/dξ + J⁻¹₂₂ dN/dη

在MATLAB里，我喜欢写一个函数返回所有9个形函数在单个积分点的值S和参考梯度dSdxi、dSdeta，物理梯度统一在外面做矩阵乘法。函数参数只要两个标量xi和eta，返回三个9×1数组。这样积分循环里调用一次就能取到全部节点信息，不会出现“形函数求值顺序和连接表顺序不一致”的问题。

网格是矩形时，J是一个对角常数矩阵，detJ等于单元实际面积除以4。这个值可以直接手算验证，比运行完整程序后再检查快得多。我通常在调试期加一句断言：如果任意单元detJ≤0，立即报错并打印单元编号。这能在早期拦截大部分由网格生成器造成的翻转单元问题。

3.3 高斯积分点选取：Q2 需要至少 3×3

二维Q2单元上速度刚度被积函数的最高多项式阶是4次，理论上每个方向用3个高斯点就足够精确积分。压力耦合项和体积力项的阶数可能更高，如果压力形函数是Q1，3×3对大多数情况已经可用；如果压力是P1不连续，某些实现会在单元内对压力做单独积分，这时也只多算几个交叉项，3×3仍然够初版使用。不要为了省计算量用2×2积分点处理Q2单元，那样会造成单元在某些变形模式下零能量，表现为稀疏矩阵的零特征值增多，速度场出现类似沙漏的伪模态。

3点高斯积分对应的坐标是±√(3/5)和0，权重分别是5/9、8/9、5/9。二维积分时每个方向遍历一遍，总权重是两方向权重相乘再乘以detJ。这套参数在大多数有限元教材里都有，但每年还是能看到有人把权重写成1，导致刚度矩阵整体偏小几个数量级，收敛阶却保持正确，这种问题最隐蔽。检查办法：用一个线性场或者二次场构造右端，直接对单单元做积分，对比解析体积。比对着代码发呆快得多。

4. 刚度矩阵组装与边界条件：MATLAB 里能直接跑的核心骨架

4.1 速度与压力的全局自由度对齐

全局自由度排序是组装时最容易出错的地方。以Q2-Q1为例，速度节点总数Nu=(2nx+1)(2ny+1)，每个速度节点有ux和uy两个分量，所以速度自由度是2Nu。压力节点总数Np=(nx+1)(ny+1)，每个压力节点一个标量自由度。我规定的排序是：

自由度1到Nu：所有速度节点的ux分量
自由度Nu+1到2Nu：所有速度节点的uy分量
自由度2Nu+1到2Nu+Np：所有压力节点的p值

这种排序下，全局矩阵自然分成左上速度刚度、右上与左下耦合、右下零压力块的鞍点结构。局部单元矩阵投射到全局前，先做局部自由度映射：

% 单元局部速度节点全局编号集合 vel_glob(1:9) vel_idof = [2*vel_glob-1, 2*vel_glob]; % 速度x和y自由度 pre_idof = 2*Nu + pre_glob; % 压力自由度 local_dofs = [vel_idof, pre_idof]; % 单元全部全局自由度

这个映射一旦写错，矩阵会乱成一团，而且不报错。建议在小网格nx=2、ny=2的规模下手动检查前几个单元的local_dofs，确保它和你想象中的节点编号顺序一致。

4.2 单元刚度组装：最小可运行骨架

下面给出一段可直接理解的结构化组装代码。它没有做向量化优化，但对单元数几千的算例仍然可运行，逻辑也最容易对照。

% 组装 Q2-Q1 Stokes 矩阵的核心循环 % 自由度排序：ux(1:Nu), uy(Nu+1:2Nu), p(2Nu+1:2Nu+Np) K = sparse(ndof, ndof); F = zeros(ndof, 1); % 高斯积分点与权重 gp = [-sqrt(3/5), 0, sqrt(3/5)]; gw = [5/9, 8/9, 5/9]; for e = 1:ne % 单元速度/压力节点全局编号 vel_glob = eU(e,:); % 1x9 pre_glob = eP(e,:); % 1x4 x_vel = U_node(vel_glob, 1); y_vel = U_node(vel_glob, 2); % 局部矩阵：速度9节点有18个速度自由度，压力4个 Ke = zeros(2*9 + 4, 2*9 + 4); for qi = 1:3 for qj = 1:3 xi = gp(qi); eta = gp(qj); w = gw(qi)*gw(qj); % 计算速度形函数和参考梯度 [S_u, dSdxi, dSdeta] = shape_q2(xi, eta); % 速度节点坐标插值得到雅可比 dXdxi = sum(x_vel .* dSdxi); dYdxi = sum(y_vel .* dSdxi); dXdeta = sum(x_vel .* dSdeta); dYdeta = sum(y_vel .* dSdeta); J = [dXdxi, dYdxi; dXdeta, dYdeta]; detJ = det(J); Jinv = inv(J); % 物理梯度 dSdx = Jinv(1,1)*dSdxi + Jinv(1,2)*dSdeta; dSdy = Jinv(2,1)*dSdxi + Jinv(2,2)*dSdeta; % 压力形函数（Q1双线性），用统一节点顺序求值 [S_p, dPdx, dPdy] = shape_q1(xi, eta); % 速度刚度块：mu*(dNi/dx*dNj/dx + dNi/dy*dNj/dy) for iu = 1:9 for ju = 1:9 a_ij = mu*(dSdx(iu)*dSdx(ju) + dSdy(iu)*dSdy(ju)); Ke(iu, ju) = Ke(iu, ju) + a_ij*detJ*w; Ke(iu+9, ju+9) = Ke(iu+9, ju+9) + a_ij*detJ*w; end end % 耦合块：Kup(q, i) = ∫ q * div(v_i) % 速度x分量对应第iu节点，散度为 dSdx(iu) % 速度y分量对应第iu+9节点，散度为 dSdy(iu) for ip = 1:4 for iu = 1:9 % Kup 列 = 速度自由度，行 = 压力自由度 Ke(2*9+ip, iu) = Ke(2*9+ip, iu) + S_p(ip)*dSdx(iu)*detJ*w; Ke(2*9+ip, iu+9) = Ke(2*9+ip, iu+9) + S_p(ip)*dSdy(iu)*detJ*w; % 对称块，注意这里的转置关系 Ke(iu, 2*9+ip) = Ke(iu, 2*9+ip) + S_p(ip)*dSdx(iu)*detJ*w; Ke(iu+9, 2*9+ip) = Ke(iu+9, 2*9+ip) + S_p(ip)*dSdy(iu)*detJ*w; end end % 体积力（此处按解析函数 f 求值，可留空） end end % 投射到全局 local_dofs = [2*vel_glob-1, 2*vel_glob, 2*Nu+pre_glob]; K(local_dofs, local_dofs) = K(local_dofs, local_dofs) + Ke; end

这段代码里，速度刚度块中ux自由度只和ux耦合、uy只和uy耦合，没有交叉项，因为拉普拉斯算子在笛卡尔坐标系下是分量独立的。耦合块的对称性来自b(v,q)的定义：Kup(q,i)=Kpu(i,q)，所以在两个下标的写入位置要保持一致，否则整个矩阵不对称，后面积分解出来结果偏差很难定位。volume force F向量我在这里留空了，实际计算时把f·v在积分点叠加到对应速度自由度上即可。

注意Kuu块中两个速度分量分开放置，每个9×9子块对应同一套速度节点顺序。这意味着如果你后来换成P1压力，只要把压力部分改成3个三角形内自由度并调整耦合循环就可以，速度部分完全不用动。把单元类型抽象成这种“局部自由度分块”结构，对后续扩展非常有用。

4.3 Dirichlet 边界与压力参考：顶盖驱动腔的最小实例

Stokes问题最适合入门验证的算例是顶盖驱动腔（lid-driven cavity）。矩形域全部边界速度给定：顶盖给水平速度1，其余三边速度为0，四条边上的所有节点都受速度本质条件约束。由于没有压力边界条件，纯Dirichlet下压力只能确定到差一个常数，所以必须施加一个压力参考约束，否则矩阵奇异。

Dirichlet边界处理的常见做法是置一法。先判断哪些全局自由度被约束，再对矩阵做行、列修正。一个稳妥的写法：

% u_bc: 被约束的全局自由度编号 % u_bc_val: 对应约束值 % b: 右端向量（组装完成后） for k = 1:length(u_bc) dof = u_bc(k); val = u_bc_val(k); % 先把该行清零并设置对角为1 K(dof, :) = 0; K(dof, dof) = 1; b(dof) = val; end % 修正其他方程：将边界自由度贡献从右端减去 for k = 1:length(u_bc) dof = u_bc(k); val = u_bc_val(k); col = K(:, dof); K(:, dof) = 0; b = b - col * val; b(dof) = val; end

这段代码必须先做列修正再清零那一列，否则col会在清零后变成0，修正无效。很多教程只写前半段行修正，导致远离边界节点上速度受边界影响缺失，表现为“边界速度对了但流场完全不对”。如果你在已经组装好的矩阵上反复加约束，建议每次运行时重新组装K，而不是在原矩阵上叠加修改，否则列修正会重复累计。

压力参考最简单的方案是选择一个压力节点，直接设置它的值为0。我通常取底部边中点：它不接近奇点，也不在角点上。如果压力用Q2-P1不连续自由度，参考约束要加在某个单元的压力分量上，否则整个系统仍然有一个对应“所有单元压力同时加常数”的零特征向量。

顶盖驱动腔的角点是另一个坑：左上角既是顶盖边界又是左壁边界，顶盖u=1和左壁u=0在这里冲突。物理上更合理的处理是所有角点速度设为0，顶盖驱动只加在除两个角点以外的顶盖节点上。这样避免了角点奇异性，压力振荡也小得多。如果硬把角点设为1，求解器不会报错，但角点附近速度梯度会异常大，后处理云图在角落出现明显的颜色断层。

4.4 求解与迭代法选择：先 A\b，再考虑 MINRES

单元规模不很大时，MATLAB反斜杠对稀疏矩阵的求解效率已经足够。以100×100网格的Q2-Q1为例，速度自由度约8万个，压力自由度约1万个，总自由度约9万，直接在稀疏矩阵上用A\b做一次分解在普通工作站上大约几十秒，内存占用也可控。反斜杠会自动判断矩阵结构，对鞍点矩阵选用合适分解器，不需要你手动指定。

网格继续加密时，直接法的时间和内存会快速上升。那时需要考虑迭代法。鞍点矩阵虽然对称但不正定，MINRES或GMRES是常见选择。预处理器建议单独构建速度块Kuu的不完全Cholesky因子，配合块对角结构做预处理。不要对整个K直接调用ichol，因为右下零块会让分解彻底失败。这一步属于性能优化，可以在程序验证通过后再加，不要一开始就让迭代法干扰正确性排查。

5. Stokes 程序常见翻车现场与排查：棋盘振荡、奇异矩阵与收敛阶

5.1 压力棋盘振荡：算得出速度，压力却像国际象棋

现象：速度云图看起来合理，流线也光滑，但压力云图呈现红蓝交替的棋盘状，且网格加密后振荡更明显，求解器却不报任何错误。这是Stokes程序里最经典的“不满足稳定性条件”的后果。原因几乎都是用了Q1-Q1或Q1-P0格式。这个组合的离散压力空间存在spurious mode，相邻压力正负交替，速度限制条件检测不到这个模式，于是压力场被零特征向量污染。解决的办法只有两条路：换Q2-Q1或Q2-P1格式；或者保留原格式并加入PSPG型稳定化项。不能靠增加网格精度救回来，因为这是一个纯代数性质的问题。

验证手段很直接：小网格上组装出全局压力-压力块，或者把整个K矩阵做奇异值分解，统计接近零的奇异值个数。Q2-Q1只应该有一个对应整体常数压力的零奇异值；Q1-Q1通常会多出若干个，而且特征向量呈棋盘分布。看到多余零奇异值就不用再查代码了，直接换格式。

5.2 MATLAB 提示矩阵奇异或警告矩阵接近奇异

现象：K\b执行时报错“Matrix is singular to working precision”，或者返回结果全为NaN，压力场尤其明显。原因有四类：纯Dirichlet边界下没有加压力参考，压力自由度整体漂移造成零空间；网格里存在“悬空节点”，即某个速度节点没有被任何单元的连接表引用，对应的整行为零；压力节点与速度节点耦合关系错误，压力自由度的某一行始终为零；Dirichlet边界处理时漏掉某条边的节点。解决步骤也按顺序来：先检查是否已钳制压力；再用spy(K)看零行零列，通常一眼就能发现悬空节点；最后逐条核对边界节点列表，看是否漏了边中点。我的习惯是在小网格nx=2、ny=2上先把自由度映射表、边界节点编号打印出来人工核对一遍，这一步看似笨拙，但能省掉后面全网格调参的几小时。

5.3 收敛阶不达标：速度二阶、压力二阶，但斜率怎么都不对

现象：解析解算例里，速度的L2误差随网格细化也在下降，但log-log斜率只有1.2左右，和理论上的2阶收敛差一截；压力误差更离谱，甚至不降。排查顺序有个固定套路：先检查高斯积分精度，Q2单元至少3×3，如果用了2×2就会出现这种问题；再检查雅可比物理梯度的变换公式，J⁻¹写反或转置错位是高频错误；再检查耦合块的符号和对称性，Kup和Kpu如果不对称，解的精度会在粗网格上凑合、细网格上原形毕露；最后检查Dirichlet列修正是否真的执行了。

有一个非常隐蔽的坑：压力钳制位置会引入压力场的整体偏移。如果误差计算直接比较p_h和p_exact，那么被钳制点处压力误差很小，远处则产生一个几乎常数偏差，导致压力误差不随网格细化下降。正确做法是先做压力均值修正，或者计算误差时把压力整体减去其均值再比较。很多收敛阶不达标的现象，最后都栽在这个细节上。

5.4 顶盖驱动算例流场不对称或角点速度异常

现象：顶盖驱动腔的流场中心没有落在几何中心，或左上角点附近速度显示明显异常，比如角点处速度矢量指向异常方向。原因大概率是角点边界条件冲突。顶盖速度为1，左壁和右壁速度为0，角点同时属于两条边。如果给角点赋值1，顶盖驱动在角点产生一个不物理的“钉扎效应”，角点速度梯度趋于无穷大；如果赋值0，顶盖的驱动速度在角点处强制降为0，这是物理上更合理的处理。解决：建立Dirichlet列表时，先用所有边界速度=0打底，再单独把顶盖上除去两个角点的节点设为1。这也是文献里通行的表述口径：cavity顶盖带角点无滑移处理。

5.5 网格变大后内存溢出或求解超时

现象：20×20网格一切正常，100×100直接Out of Memory，或A\b超时。原因通常是两个：一是用zeros(ndof,ndof)生成了全稠密矩阵再往里填，内存直接爆掉；二是边长数量增加后总自由度数增长二次方，直接分解计算量增长三次方，超时是自然结果。解决：组装全程用sparse预分配，绝不在中途转full；如果仍然慢，就把网格压回50×50验证正确性，然后上迭代法。迭代法阶段，先做Kuu的不完全Cholesky预处理，再跑MINRES。不要一开始就上并行或GPU，这类问题在几千单元规模下用迭代法配合预处理通常几十秒内就能收敛。

6. 用解析解验证程序：制造一个你提前知道答案的算例

6.1 构造散度自由的速度场与对应压力

算例验证是Stokes程序调试的临门一脚。最可靠的做法是设计一组速度场和压力场，代入Stokes方程反算体积力f，然后在有限元模型里把这个f作为右端载荷，让程序去还原原来的速度压力。一个常用的解析构造是选原始流函数：

ψ(x,y) = x²(1-x)² y²(1-y)²

速度取旋度：

u = ∂ψ/∂y
v = -∂ψ/∂x

这个速度场自动满足∇·u=0。把它代入矩形域[0,1]×[0,1]，四条边界上恰好全部速度为0，所以Dirichlet边界条件就是零向量，实现起来极其省事。压力直接取p=0或者p=x+y都可以，后者能让压力梯度进入右端体积力，检验耦合项是否正确。然后反算：

f = -μΔu + ∇p

把f按解析函数形式在积分点上求值，叠加到单元右端向量。这一步一旦程序跑通，误差对比立刻能暴露组装环节的所有问题。

6.2 误差度量与收敛阶：log-log曲线斜率才是验收标准

有限元解的误差度量要按域内积分来完成，不能只挑几个节点比较。对速度场，计算：

e_u = sqrt(∫ |u_h - u_exact|² dΩ)

压力场同理。在MATLAB里这个积分仍然用高斯点做，网格加密后分别记录误差。Q2-Q1的理论期望是速度L2误差收敛阶2，压力L2误差也至少收敛阶1到2。实测时在h从1/8细化到1/16、1/32的过程中，log-log斜率应当稳定接近2。如果斜率明显低于理论值，不要硬调参数，先回到第5章排查流程。

一个更细的验收标准是分块误差：单独看耦合项和速度刚度块的贡献可以通过“冻结压力”的方式解耦测试。比如只组装速度刚度矩阵，解一个已知解析解的Poisson问题；只组装压力耦合块，用已知散度场检验Kup的转置关系。每个块都独立通过后，再合在一起跑Stokes。这个分块测试习惯能帮你把排查范围缩小到某一个具体函数，而不是面对一整个矩阵怀疑人生。

6.3 程序结构整理与复用：从验证算例到工程骨架

我建议把程序拆成setup_mesh、assemble_stokes、apply_bc、solve_stokes四个独立函数，主脚本只负责串流程和写结果。这样换一个算例，只需要修改setup_mesh的边界条件和assemble_stokes里的右端载荷，组装核心代码完全不动。再进一步，把单元类型做成一个字符串参数，比如element='Q2Q1'或element='Q2P1'，同一个主程序可以切换格式，方便对比稳定化效果和收敛阶。

我最早写Stokes程序时也固执地选了Q1-Q1，压力棋盘振荡反复出现，我足足查了两天代码，一度怀疑是求解器的收敛容差问题。最后被逼着去读混合元稳定性的讨论，才意识到是格式选错，这个教训让我养成了习惯：任何新网格、新格式都先用解析算例做分块验证，再跑真实算例；在代码头部注释里固定写出弱形式的符号约定和自由度排序，防止隔几个月回来看不懂自己的矩阵。这套流程虽然看起来多花半天，但能省掉后续无数个“为什么压力又振荡”的深夜。希望帮到你。

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