简介：二维稳态Navier-Stokes方程有限元求解程序是一份基于Matlab的数值计算源码，面向流体力学、计算数学方向的学生与研究者，用于求解不可压缩黏性流动的速度场与压力场。资源共45个m文件，压缩包仅19KB，全部为Matlab源代码，模块覆盖均匀网格生成、有限元基函数与高斯积分、刚度矩阵组装、Dirichlet与Robin边界条件处理、牛顿线性化求解、压力点固定及Hs误差计算等，结构清晰，方便按需阅读与二次开发。作者Jack_Lue提供，已有312人学习/下载。程序实现了从几何离散、变分弱形式、矩阵组装到线性系统求解的完整流程，并附有多个解析解与数值误差算例；读者可通过阅读代码深入理解有限元求解流体力学的关键步骤，稍作修改即可扩展至热传导、扩散方程等类似偏微分方程，兼具教学演示与工程参考价值。

1. 二维稳态Navier-Stokes方程有限元求解程序到底在解什么

第一次把二维稳态Navier-Stokes方程有限元求解程序跑通时，你会清楚记得那个瞬间：屏幕上画出的方腔流场中心涡，终于和文献里的基准解重合了。这个程序听起来像个CFD课程的压轴作业，但它实际要跨三座坎——处理速度与压力的混合插值、对非线性对流项做迭代、维护鞍点矩阵的收敛可靠性，顺序和优先级错一个，结果都不会报错，只会安静地错得离谱。这篇文章写给需要亲手写一版、而且想把边界条件、压力参考点、迭代序列这些细节一次弄明白的工程师和研究生。你会看到为什么有限元在不可压缩流里不能照着弹性力学那套抄，也会拿到一个能跑通的最小骨架和一组照着填就能用的参数。

2. 从连续方程到离散系统：混合有限元为什么选Q2-Q1组合

2.1 二维稳态N-S方程的弱形式：三个算子同时要处理

二维不可压缩稳态流动的控制方程，写成强形式是这样的：连续方程 ∇·u = 0，动量方程 (u·∇)u − (1/Re)∇²u + ∇p = f。这套方程跟固体力学里的弹性问题有本质区别：弹性问题的刚度矩阵对称正定，而无粘对流项 (u·∇)u 带来了非对称性，压力与速度的耦合又让代数系统变成鞍点结构。换句话说，有限元在离散这套方程时，抽象矩阵里同时出现了三块性质完全不同的算子：速度梯度对速度梯度的扩散积分、速度对速度的对流输运、以及压力和速度之间的散度耦合块。如果还抱着弹性力学那套“组装刚度阵、求解、后处理”的惯性思维，第一轮组装出来的矩阵一定是奇异的，或者解出一片毫无意义的振荡。

把动量方程乘以速度测试函数 v，连续方程乘以压力测试函数 q，做分部积分之后，弱形式变成下面两个等式。速度测试函数空间是 H0^1，压力测试函数空间是 L0^2，注意压力只要求平方可积，不要求导数平方可积，这是混合元跟位移元最大的分界线：

∫(u·∇)u·v dx + (1/Re)∫∇u:∇v dx − ∫p(∇·v)dx = ∫f·v dx ∫q(∇·u)dx = 0

把这个弱形式离散化之后，你得到的是形如 [ A+B Cᵀ; C 0 ] 的鞍点矩阵，其中 A 来自扩散项，B 来自对流项，C 来自散度约束。关键点在于右下角那个 0 块——它意味着这个系统不是正定的，不能用共轭梯度法直接套，也意味着压力和速度的插值空间必须满足一个叫 inf-sup 条件的匹配关系。这个条件看起来像纯数学的产物，实际上它直接决定你的压力场能不能解出来。

2.2 速度压力空间的匹配：为什么Q1-Q1必然翻车

我见过不止一个新手用四个节点的双线性单元同时插值速度和压力，也就是 Q1-Q1，组装出来的矩阵能解，速度场看着也像回事，但把压力场画出来就露馅了：压力呈红蓝棋盘状交替跳动，相邻单元的p值锯齿形错开，物理上完全扯淡。原因在于 Q1-Q1 不满足 inf-sup 条件（也叫 LBB 条件），离散连续方程的约束矩阵 C 存在大量虚假零空间，除了物理上允许的“整体常数压力”之外，还多出了高频振荡模态。

解决办法有两类。第一类是换插值空间：给速度用二阶单元、压力保持一阶，典型就是本文要用的 Q2-Q1（Taylor-Hood 元），速度在每个大单元上有9个节点，压力在每个大单元上只取4个角点。这类配对满足 inf-sup 条件，不需要额外加任何修正项，压力场干净。第二类是用等阶单元再附加人工耗散项来消除振荡，工程软件里很常见，但调参数的玄学成分大，自己写代码时容易变成靠试。所以这种入门级稳态N-S程序，我一般直接锁定 Q2-Q1，少一个变量少一堆坑。

2.3 对流项线性化：Picard迭代与牛顿迭代的取舍

稳态N-S方程里的 (u·∇)u 是未知量自己乘自己的梯度，是个非线性项。有限元离散后得到 F(u,p) = 0 的方程组，要求解只能迭代。常见做法有两种：Picard 迭代是把上一轮的已知速度 uⁿ 当作对流速度，然后求解一个线性系统，配下方带 (uⁿ·∇)uⁿ⁺¹；牛顿迭代则是把整项做泰勒展开，保留未知量微增量的导数项 δu·∇uⁿ。Picard 的收敛相对更从容，初始解离真实解远的时候不容易发散，但收敛到解的速率是线性的；牛顿迭代是二阶收敛，一旦进入 converge radius 会非常快，但它娇气，初值不好会直接给你蹦出 NaN。

所以实践里更多用的是“先 Picard 后牛顿”的切换策略。低雷诺数比如 Re=100，直接用牛顿、初值给全零，通常也扛得住。雷诺数超过一千之后，我一般先跑十几步 Picard 把流场轮廓暖起来，再切牛顿把残差压到机器精度附近。这个切换是这套程序里性价比最高的一个操作，后面第3章的迭代循环会直接按这个结构写。

3. 用MATLAB把有限元求解步骤落成最小骨架

3.1 结构化网格生成与自由度编号

先处理几何和编号。为了不引入网格生成器的黑匣子，直接在 [0,1]×[0,1] 方腔上生成结构化四边形网格。Q2 单元每个方向上是二次映射，所以每个大单元内部沿边长有2个节点间隔，节点总数是 (2nx+1)×(2ny+1)。下面这段代码生成节点坐标，并把每个单元的9个节点号按固定顺序填入 elem 矩阵：

% 方腔网格生成：Q2 单元，每个大单元 9 个节点 nx = 16; ny = 16; % 期望的单元数，先取小值跑通，再加密 npx = 2*nx + 1; npy = 2*ny + 1; x = linspace(0, 1, npx); y = linspace(0, 1, npy); [X, Y] = meshgrid(x, y); node = [X(:), Y(:)]; % node(:,1)为x坐标，node(:,2)为y坐标 % 节点排布规则：meshgrid 展开后，沿x方向变化最快 % 每个大单元的9个节点顺序：左下、右下、右上、左上、下边中、右边中、上边中、左边中、中心 n_elem = nx * ny; elem = zeros(n_elem, 9); for ej = 1:ny for ei = 1:nx % 大单元角点在稀疏网格里的坐标索引（2倍关系） i0 = 2*(ei-1) + 1; % 左边界节点小坐标 j0 = 2*(ej-1) + 1; % 下边界节点小坐标 % 角点 e = zeros(1,9); e(1) = (j0-1)*npx + i0; % 左下 e(2) = (j0-1)*npx + i0 + 2; % 右下 e(3) = (j0+2-1)*npx + i0 + 2; % 右上 e(4) = (j0+2-1)*npx + i0; % 左上 % 边中点和中心 e(5) = (j0-1)*npx + i0 + 1; % 下边中点 e(6) = (j0+1-1)*npx + i0 + 2; % 右边中点 e(7) = (j0+2-1)*npx + i0 + 1; % 上边中点 e(8) = (j0+1-1)*npx + i0; % 左边中点 e(9) = (j0+1-1)*npx + i0 + 1; % 中心 elem((ej-1)*nx + ei, :) = e; end end

这里最关键的是编号约定。速度自由度直接挂在 node 上，每个节点两个分量，所以速度自由度数 ndof_v = 2*size(node,1)。压力采用 Q1，每个大单元只取四个角点，也就是一个单元只贡献一个内压力自由度，ndof_p = n_elem。如果你想把压力也按节点存，记得压力节点数少于速度节点数，组装时下标索引必须单独建一套，否则速度下标和压力下标混用一个编号表，矩阵装配必乱。这段生成代码跑完之后，建议先画一下 elem 里第1个单元的9个点坐标，确认顺序是逆时针，后面积分才不会出符号问题。

3.2 单元局部矩阵装配与高斯积分

局部矩阵的组装是所有有限元程序的核心体力活。Q2 速度单元在等参坐标系下形函数有9个，压力 Q1 形函数只有4个，它们的导数要通过雅可比矩阵从等参坐标变换到物理坐标。MATLAB 里最省事的写法是每个单元内循环 3×3 个高斯点，逐点累加局部矩阵。下面给的是单个单元局部矩阵的骨架函数，注意对流项里要传入当前迭代步的速度解：

function Ke = local_NS(elem, node, u_all, p_all, Re) % elem: 当前单元9个节点编号; u_all: 当前速度自由度值; p_all: 当前压力值 % 返回局部矩阵 Ke，尺寸为 (18+4)x(18+4)，顺序先速度后压力 ng = 3; [gx, gw] = gauss_rule(ng); % 3x3高斯点与权重 Ke = zeros(22, 22); % 当前单元速度在节点上的两个分量，用于组装对流项系数 ue = zeros(9,1); ve = zeros(9,1); for a = 1:9 ue(a) = u_all(2*(elem(a)-1)+1); ve(a) = u_all(2*(elem(a)-1)+2); end for i = 1:ng for j = 1:ng [N, dN_dxi] = shape_Q2(gx(i), gx(j)); % 速度形函数，1x9 [Np, dNp_dxi] = shape_Q1(gx(i), gx(j)); % 压力形函数，1x4 J = dN_dxi' * node(elem(1:9), :); % 2x2雅可比 dN_dx = dN_dxi / J; % 物理坐标下的速度形函数导数 dNp_dx = dNp_dxi / J; % 压力形函数导数 detJ = abs(det(J)); % 速度-速度块：粘性 + 对流 for a = 1:9 for b = 1:9 K_ab = (1/Re)*(dN_dx(a,:)*dN_dx(b,:)') ... + (ue'*N(:,a)+ve'*N(:,a)) * (N(b)*dN_dx(a,1) + N(b)*dN_dx(a,2)); % 这里涉及张量积，实际实现时按分量展开 % 速度自由度分量的索引为 a_u和 b_u Ke(2*a-1, 2*b-1) = Ke(2*a-1, 2*b-1) + K_ab*gw(i)*gw(j)*detJ; Ke(2*a, 2*b ) = Ke(2*a, 2*b ) + K_ab*gw(i)*gw(j)*detJ; end end % 压力-速度块 for a = 1:9 for b = 1:4 C_ab = -dNp_dx(b,1)*N(a) ; % -∫∂q/∂x * v_x D_ab = -dNp_dx(b,2)*N(a) ; Ke(2*a-1, 18+b) = Ke(2*a-1, 18+b) + C_ab*gw(i)*gw(j)*detJ; Ke(2*a, 18+b) = Ke(2*a, 18+b) + D_ab*gw(i)*gw(j)*detJ; end end end end % 对称部分由散度项转置补全，组装时 Ke(19:22,1:18) 直接取转置 end

这段代码为了可读性省略了张量积的完整展开，实际实现时要注意几个点。第一，高斯积分权重 gw 是一维权重，二维组合要用 gw(i)*gw(j)，如果漏了乘积，矩阵数值会差一个权重因子，收敛阶直接掉。第二，雅可比矩阵 J 在这里是一个对角矩阵（结构化矩形网格），但代码不要写成只对矩形有效的硬编码，等你换到非结构网格时雅可比会变成全矩阵，现在写通用形式剩得下以后的改造成本。第三，对流项的积分符号取决于你对双线性形式的定义，不同教科书约定不同，我这里的写法是 ∫(u·∇u)·v，对应到代码里需要把 u 和 v 的方向分量配对好，配对错一个正负号，Re大的时候矩阵就不占对角优势，迭代会莫名其妙发散。

3.3 Picard到牛顿切换的主迭代循环

外层迭代只需要做一件事：装配全局矩阵、解线性系统、更新解、判断残差。注意稳态N-S的残差不是常说的 F = A*u - b，而是 A(u)*u - b，因为A本身也依赖当前 u。这一点很多人会在残差判据上犯迷糊。下面是完整循环体的骨架：

% 主循环：先Picard后牛顿 u = zeros(ndof_v, 1); p = zeros(ndof_p, 1); alpha = 1.0; % 松弛因子，小于1时更稳但更慢 maxit = 100; tol = 1e-8; switch_step = 20; % 前20步用Picard，之后切牛顿 for iter = 1:maxit A = assemble_global(elem, node, u, p, Re); F = assemble_residual(elem, node, u, p, Re); % F = A(u)*u - b [A, F] = apply_bc(A, F); % 强加Dirichlet边界 du = A \ F; u = u + alpha * du(1:ndof_v); p = p + alpha * du(ndof_v+1:end); res = norm(F, 2); if res < tol break; end if iter == switch_step % 切到牛顿：此时雅可比矩阵里要额外加入对流项的导数项 newton_on = true; end end

这段代码里最容易被忽略的是 apply_bc 和 du 的维度。N-S问题里速度和压力都有约束：速度的Dirichlet条件用置大数法或划行划列法直接改矩阵和右端项，压力则要额外固定一个点的压力值，否则 A 在常数压力方向上是奇异的，MATLAB 解出来直接是无穷大。我的经验是每次迭代时把压力参考点那行的自由度和右侧同时清零，对角线置1，保证常数压力模态被钉住。第三个点是 alpha 取1在Re不太高的方腔上一般能走完，但Re到5000以上时，Picard阶段alpha取0.5左右更稳，切牛顿后再回到全步长。这个方法不是唯一标准答案，但它是这类问题里最不容易翻车的启动序列。

4. 边界条件与参数选择：方腔算例的雷诺数、松弛因子与收敛判据

4.1 方腔驱动流的标准边界设置与无量纲化

方腔驱动流是N-S有限元程序最合适的测试台，因为它的边界条件简单到只有一条带速度的边，却依然能激励出完整的涡结构。标准做法如下：腔体边长取1，顶盖（y=1）水平速度 ux=1，垂直速度 uy=0；其余三条边速度全为零；压力不设置边界条件，由连续方程隐式确定，但整个系统的压力只定义到相差一个常数。不要小看这个“相差常数”，它是程序里最隐蔽的雷。

下表是推荐的无量纲参数设置：

无量纲化值得多说一句。如果你把程序写成带单位的参数，后面想做无量纲基准解对比时会非常痛苦。我一般在一开始就把长度除以L、速度除以顶盖速度，这样Re成为方程里唯一的无量纲参数，边界条件也变成纯1和0。单位制清晰，后续判断网格分辨率是否足够也方便——高Re下边界层厚度量级约是1/√Re，网格最密处是否需要小于边界层厚度，一眼就能算出来。

4.2 施加边界条件与压力参考点的四种常见做法

强加Dirichlet边界条件的实现方式有直接消去法、置大数法、罚函数法、以及投影法。针对稳态N-S有限元，最稳的是“先组装后修改”的直接消去法思路：把速度边界节点对应的行列从全局矩阵中划出，再把对应的解分量直接赋值为已知边界值。实现起来麻烦些，但矩阵条件数不会改变。置大数法实现简单，可对内部节点和边界节点混在一起的鞍点系统来说，大数到底取多大本身是个玄学——取太小约束不紧，取太大把条件数搞坏，牛顿迭代尾期残差降不下去。

我在这套程序里更喜欢用划行划列，但对压力参考点用对角线置1的方式。原因很简单，压力参考点不是物理边界，它是消除零空间的手段，对角线置1既不改变矩阵稀疏结构，也不会引入额外误差。代码如下：

function [A, b] = apply_bc(A, b, dof_v, bc_nodes) % 速度边界：划行划列，直接替换 free_dofs = setdiff(1:size(A,1), bc_nodes); A_in = A(free_dofs, free_dofs); b_in = b(free_dofs); % 压力参考点：设压力自由度编号为 p_ref p_ref = ndof_v + 1; % 取第一个压力自由度作为参考 A_in(end+1, :) = 0; A_in(:, end+1) = 0; % 扩展矩阵示意 % 实际实现：把 p_ref 所在的行列置0，再令对角线=1、右端=0 A(p_ref, :) = 0; A(:, p_ref) = 0; A(p_ref, p_ref) = 1; b(p_ref) = 0; end

这段代码有一处需要特别留意：压力参考点不能选在边界速度约束与压力空间重叠的区域，否则会造成一个点是速度被钉死、压力也固定的双重约束。在Q2-Q1结构里，压力节点与速度角点位置重合，我通常选腔体左下角的速度角点对应的压力自由度作为参考点，这样它本身就是速度Dirichlet边界节点，双重约束反而不产生额外影响，因为在划行列时它已经被处理过。边界条件的施加顺序也很重要——必须先做速度约束，再做压力参考点约束，顺序反过来，速度边界条件会被压力参考点的置零操作污染。

4.3 松弛因子、迭代初值与雷诺数续化

越是新手越容易忽视迭代初值，总觉得从全零开始迭代天经地义。对Re=100的层流方腔，全零初值确实没问题，但Re超过1000后，全零初值等于把牛顿迭代扔进一个完全陌生的区域，前几步残差可能正常下降，随后突然爆成Inf。这个现象的根源是对流项雅可比矩阵在远离解的地方不具备收敛性，所以工程上常用“续化法”：从小到大依次解Re=100、400、1000、2000，把上一个雷诺数的收敛解当作下一个雷诺数的迭代初值。每步之间Re的增幅一般不超过2倍，加密网格时也要同步降低Re增幅。

松弛因子和续化法配合起来，才是一套完整的起手式。Step1用alpha=0.5跑Picard，Step2切牛顿时alpha回到1.0，Step3在一切顺利的前提下把Re翻倍。如果某一步发散，优先做两件事：降alpha到0.3，或者把Re增幅减半，而不是去调网格、调收敛容差。这些参数之间的耦合有个经验关系：Re越高，能承受的alpha越小，牛顿迭代的初值敏感度也越高。收敛判据这边，我习惯用“先看残差范数、再看增量范数”双条件，残差降到1e-6以下且增量比残差小一个量级，才认为收敛；只靠残差一个量容易过早停下，因为鞍点系统的残差在某些迭代阶段会有一条短暂的平台期。

5. 常见问题排查：压力振荡、对流占优与迭代发散

5.1 压力场出现棋盘振荡

现象：速度场看起来收敛了，顶盖驱动流的主涡也在，但是把压力场用云图画出来，会出现红蓝相间的棋盘状斑块，压力最大值和最小值交替出现，完全不符合物理。原因基本锁定为速度与压力插值不匹配，典型就是用了Q1-Q1等阶单元，离散连续方程不满足inf-sup条件；长度尺度上新单元里虚假压力模态与真实压力场的频率混在一起，肉眼看起来就是噪声。解决：要么改用Q2-Q1混合元，要么给等阶单元附加额外耗散项。自己写程序，我建议直接换Q2-Q1，因为那一步改动只在形函数和编号上，后面所有代码都能复用；加耗散项则牵涉稳定参数的选择，那个参数在不同网格和Re下都要调，教训是用等阶元调试压力场的夜晚永远比想象的多。

5.2 高雷诺数下牛顿迭代直接发散

现象：迭代前五步残差从1e-2降到1e-4，看起来一切正常，第六步残差突然跳到1e10，屏幕上出现NaN，或者解向量里出现一长串天文数字。原因：Re升高后对流项主导，雅可比矩阵非对称性增强，初值落在牛顿迭代收敛半径之外；更隐蔽的原因是对流项的双线性形式在离散后不满足离散对流占优条件，也就是网格Pe数过大。解决：第一层手段是前面讲的续化法，把Re一步步加上去；第二层手段是切到Picard迭代跑一段时间，Picard的收敛域比牛顿大得多，相当于先让速度场进入大致正确的轨道；第三层手段是加密网格，让Pe数降下来。一个实用的判断标准是看网格最细处的单元尺寸h是否满足 h*Re<2 左右，不满足就先加密，别急着调求解器。

5.3 线性系统奇异或条件数爆炸

现象：MATLAB 在解 A\b 时给出“Matrix is singular”的警告，或者解出来的压力场整体偏一个恒定常数，速度场却正确。原因：这是最典型的压力参考点缺失。速度边界约束已经把所有刚体模态消除了，但压力常数模态还没钉住，系统沿常数压力方向奇异。解决方法就是第4章讲过的压力参考点技术。另一个常见原因是速度边界节点重复约束，比如把顶盖边界的两个角点同时写进两条边的边界节点列表，导致方程组里有两个完全相同形式的方程，矩阵秩缺1。这个排查起来很费劲，但解决办法简单——在组装边界节点时先 unique 一下，再把节点坐标编号对应的边界速度值做一致性检查。

5.4 收敛结果随网格变化剧烈，加密后解却“漂移”

现象：用16×16网格算出中心涡的位置在(0.62,0.73)，加密到32×32后位置变成(0.55,0.68)，再过细又变了。原因：这种情况多半是边界条件被施加到了错误的自由度上。例如顶盖驱动流里，顶盖的节点集合不止包含顶部一条线的角点，还包括上边中点；如果你只把角点设为ux=1而漏了边中点，Q2单元的二次边在顶部就会产生速度波动，网格加密一换，波动形态全变。另一个原因是高斯积分阶不足，N-S矩阵里的对流项含速度梯度，被积函数阶次比弹性问题高，3×3高斯积分只是起步，有些单元用3点积分在变形网格上会掉精度。解决：检查边界自由度的覆盖范围，把每个边界单元的全部相关节点列出来核对速度值；随后做一次网格收敛性研究，用16、32、64三套网格看某个特征量（比如中心涡强度）是否单调逼近一个值。

5.5 超过一定网格规模后MATLAB直接法内存崩溃

现象：网格从32×32加密到64×64，程序跑起来后内存占用飙升，MATLAB直接报“Out of memory”或者把系统swap填满。原因：N-S问题的自由度增长速度是二维的，速度自由度为2×(2n+1)²量级，压力自由度为n²量级，总自由度数约9n²，对应的稀疏矩阵带宽也随n增大。直接法 A\b 在求解过程中会产生大量填充元，矩阵非零元数量可能比原矩阵多几十倍，这是直接法的固有属性，不是代码写错了。解决：换迭代求解器，比如GMRES预处理加ILU分解，但这类鞍点问题需要有效的预处理子，简单ILU在Re高时收敛也慢。折中方案是先用一个中等网格把整个流程验证完毕，再做最终大规模网格计算；另一个实用的做法是把压力自由度重新排序到速度自由度之后，并采用嵌套剖析排序，能显著减小填充元数量。如果对迭代法不熟，不要硬上，用直接法跑一个40×40左右的网格完成论文级验证是完全足够的。

6. 验证程序正确性：用方腔基准解与网格加密判断可靠性

6.1 对照Ghia基准解：看中心涡位置和速度剖面

程序写好之后，第一件事不是跑新问题，而是跑标准算例。方腔驱动流的对照基准是Ghia在1982年给出的Re=100、400、1000的数值解，这组结果被引用了四十多年，是这类程序真正的试金石。对照方法很简单：取经过腔体几何中心的两条垂线/水平线，输出u和v速度分量，与基准表中的数据点对比。Re=100时中心涡的坐标大约在(0.6,0.73)，Re=1000时主涡中心移到约(0.53,0.56)，同时左下角和右下角会出现两个次级涡。如果你的程序在这些定性特征上差了半点，先别急着调参，去检查边界条件列表和压力参考点——这两处的bug造成的涡位置偏差最大。

6.2 网格加密测收敛阶：误差指数的判断方法

验证程序的量化手段是看收敛阶。给定一个光滑解或者基准解 reference，三套网格的误差按 h 的幂次下降。Q2速度元理论收敛阶是3阶（L2模），Q1压力元是2阶，但由于鞍点系统耦合，实测通常混在2到3之间。计算方式如下：

% 三套加密网格的L2误差分析 err = [3.2e-3, 8.0e-4, 2.0e-4]; % 举例：对应nx=16,32,64 h = [1/16, 1/32, 1/64]; p_est = log(err(1:end-1)./err(2:end)) ./ log(h(1:end-1)./h(2:end)); disp(p_est); % 输出约为 [2.0, 2.0] 左右，接近理论阶次

这里有几个实践细节。误差计算时注意统一对比的变量——速度u的L2模、压力p的L2模要分开算，压力误差的收敛阶通常会比速度低半阶，这是混合元的正常现象；如果压力误差比速度误差还高，反而说明测试解选择不对或者边界处理有问题。网格加密时要保持Re不变，只变h，否则对流项与扩散项的相对占比被改变，收敛阶自然跑偏。我还习惯在同一套网格下跑两个Re值，比如100和1000，观察误差是否随Re升高而增大——如果Re=1000的误差比Re=100大一个数量级以上，基本说明网格分辨率不够，不要归咎于求解器。

6.3 一个养成很久的习惯：先按特征量验收，再信局部细节

我给这类程序做验收时，有个固定动作：先记录中心涡位置、次级涡强度和速度剖面峰值这三个全局特征量，再去看流线图云图。全局特征量对边界条件、压力参考点的错误不敏感，但正是这种不敏感导致它们不是万能。所以下一步我会挑一个局部细节——比如顶盖边界下方第一层单元的速度梯度——单独输出数值检查，这个量的错误率远高于全局量，适合做程序bug的照妖镜。这十几年写有限元程序的经验归纳成一个习惯就是：永远假设自己的编号和边界条件有一次手滑，所以任何程序的第一版验收都以利用网格加密的对比结果为准，而不是以单套网格的漂亮云图为准。这样即便某一步疏忽了，也能通过收敛趋势把它抓出来。希望帮到你。

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