简介二维稳态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把这个弱形式离散化之后你得到的是形如 [ AB Cᵀ; C 0 ] 的鞍点矩阵其中 A 来自扩散项B 来自对流项C 来自散度约束。关键点在于右下角那个 0 块——它意味着这个系统不是正定的不能用共轭梯度法直接套也意味着压力和速度的插值空间必须满足一个叫 inf-sup 条件的匹配关系。这个条件看起来像纯数学的产物实际上它直接决定你的压力场能不能解出来。2.2 速度压力空间的匹配为什么Q1-Q1必然翻车我见过不止一个新手用四个节点的双线性单元同时插值速度和压力也就是 Q1-Q1组装出来的矩阵能解速度场看着也像回事但把压力场画出来就露馅了压力呈红蓝棋盘状交替跳动相邻单元的p值锯齿形错开物理上完全扯淡。原因在于 Q1-Q1 不满足 inf-sup 条件也叫 LBB 条件离散连续方程的约束矩阵 C 存在大量虚假零空间除了物理上允许的“整体常数压力”之外还多出了高频振荡模态。解决办法有两类。第一类是换插值空间给速度用二阶单元、压力保持一阶典型就是本文要用的 Q2-Q1Taylor-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 后牛顿”的切换策略。低雷诺数比如 Re100直接用牛顿、初值给全零通常也扛得住。雷诺数超过一千之后我一般先跑十几步 Picard 把流场轮廓暖起来再切牛顿把残差压到机器精度附近。这个切换是这套程序里性价比最高的一个操作后面第3章的迭代循环会直接按这个结构写。3. 用MATLAB把有限元求解步骤落成最小骨架3.1 结构化网格生成与自由度编号先处理几何和编号。为了不引入网格生成器的黑匣子直接在 [0,1]×[0,1] 方腔上生成结构化四边形网格。Q2 单元每个方向上是二次映射所以每个大单元内部沿边长有2个节点间隔节点总数是 (2nx1)×(2ny1)。下面这段代码生成节点坐标并把每个单元的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) (j02-1)*npx i0 2; % 右上 e(4) (j02-1)*npx i0; % 左上 % 边中点和中心 e(5) (j0-1)*npx i0 1; % 下边中点 e(6) (j01-1)*npx i0 2; % 右边中点 e(7) (j02-1)*npx i0 1; % 上边中点 e(8) (j01-1)*npx i0; % 左边中点 e(9) (j01-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尺寸为 (184)x(184)顺序先速度后压力 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, 18b) Ke(2*a-1, 18b) C_ab*gw(i)*gw(j)*detJ; Ke(2*a, 18b) Ke(2*a, 18b) 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_v1: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顶盖y1水平速度 ux1垂直速度 uy0其余三条边速度全为零压力不设置边界条件由连续方程隐式确定但整个系统的压力只定义到相差一个常数。不要小看这个“相差常数”它是程序里最隐蔽的雷。下表是推荐的无量纲参数设置参数取值说明腔体尺寸1×1无量纲化后不再关注物理尺寸顶盖速度ux1, uy0方腔的唯一边界驱动雷诺数Re100 起步1000 进阶ReUL/νν1/Re网格16×16 至 64×64先16跑通再加密看变化压力参考点某个固定节点置0消除常数压力自由度收敛容差1e-8残差范数稳态问题值得压到机器精度附近无量纲化值得多说一句。如果你把程序写成带单位的参数后面想做无量纲基准解对比时会非常痛苦。我一般在一开始就把长度除以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(end1, :) 0; A_in(:, end1) 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 松弛因子、迭代初值与雷诺数续化越是新手越容易忽视迭代初值总觉得从全零开始迭代天经地义。对Re100的层流方腔全零初值确实没问题但Re超过1000后全零初值等于把牛顿迭代扔进一个完全陌生的区域前几步残差可能正常下降随后突然爆成Inf。这个现象的根源是对流项雅可比矩阵在远离解的地方不具备收敛性所以工程上常用“续化法”从小到大依次解Re100、400、1000、2000把上一个雷诺数的收敛解当作下一个雷诺数的迭代初值。每步之间Re的增幅一般不超过2倍加密网格时也要同步降低Re增幅。松弛因子和续化法配合起来才是一套完整的起手式。Step1用alpha0.5跑PicardStep2切牛顿时alpha回到1.0Step3在一切顺利的前提下把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*Re2 左右不满足就先加密别急着调求解器。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)再过细又变了。原因这种情况多半是边界条件被施加到了错误的自由度上。例如顶盖驱动流里顶盖的节点集合不止包含顶部一条线的角点还包括上边中点如果你只把角点设为ux1而漏了边中点Q2单元的二次边在顶部就会产生速度波动网格加密一换波动形态全变。另一个原因是高斯积分阶不足N-S矩阵里的对流项含速度梯度被积函数阶次比弹性问题高3×3高斯积分只是起步有些单元用3点积分在变形网格上会掉精度。解决检查边界自由度的覆盖范围把每个边界单元的全部相关节点列出来核对速度值随后做一次网格收敛性研究用16、32、64三套网格看某个特征量比如中心涡强度是否单调逼近一个值。5.5 超过一定网格规模后MATLAB直接法内存崩溃现象网格从32×32加密到64×64程序跑起来后内存占用飙升MATLAB直接报“Out of memory”或者把系统swap填满。原因N-S问题的自由度增长速度是二维的速度自由度为2×(2n1)²量级压力自由度为n²量级总自由度数约9n²对应的稀疏矩阵带宽也随n增大。直接法 A\b 在求解过程中会产生大量填充元矩阵非零元数量可能比原矩阵多几十倍这是直接法的固有属性不是代码写错了。解决换迭代求解器比如GMRES预处理加ILU分解但这类鞍点问题需要有效的预处理子简单ILU在Re高时收敛也慢。折中方案是先用一个中等网格把整个流程验证完毕再做最终大规模网格计算另一个实用的做法是把压力自由度重新排序到速度自由度之后并采用嵌套剖析排序能显著减小填充元数量。如果对迭代法不熟不要硬上用直接法跑一个40×40左右的网格完成论文级验证是完全足够的。6. 验证程序正确性用方腔基准解与网格加密判断可靠性6.1 对照Ghia基准解看中心涡位置和速度剖面程序写好之后第一件事不是跑新问题而是跑标准算例。方腔驱动流的对照基准是Ghia在1982年给出的Re100、400、1000的数值解这组结果被引用了四十多年是这类程序真正的试金石。对照方法很简单取经过腔体几何中心的两条垂线/水平线输出u和v速度分量与基准表中的数据点对比。Re100时中心涡的坐标大约在(0.6,0.73)Re1000时主涡中心移到约(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]; % 举例对应nx16,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升高而增大——如果Re1000的误差比Re100大一个数量级以上基本说明网格分辨率不够不要归咎于求解器。6.3 一个养成很久的习惯先按特征量验收再信局部细节我给这类程序做验收时有个固定动作先记录中心涡位置、次级涡强度和速度剖面峰值这三个全局特征量再去看流线图云图。全局特征量对边界条件、压力参考点的错误不敏感但正是这种不敏感导致它们不是万能。所以下一步我会挑一个局部细节——比如顶盖边界下方第一层单元的速度梯度——单独输出数值检查这个量的错误率远高于全局量适合做程序bug的照妖镜。这十几年写有限元程序的经验归纳成一个习惯就是永远假设自己的编号和边界条件有一次手滑所以任何程序的第一版验收都以利用网格加密的对比结果为准而不是以单套网格的漂亮云图为准。这样即便某一步疏忽了也能通过收敛趋势把它抓出来。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/30 10:17:15 · 阅读 1.2k 二维稳态Navier-Stokes方程有限元求解程序构建详解
简介二维稳态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把这个弱形式离散化之后你得到的是形如 [ AB Cᵀ; C 0 ] 的鞍点矩阵其中 A 来自扩散项B 来自对流项C 来自散度约束。关键点在于右下角那个 0 块——它意味着这个系统不是正定的不能用共轭梯度法直接套也意味着压力和速度的插值空间必须满足一个叫 inf-sup 条件的匹配关系。这个条件看起来像纯数学的产物实际上它直接决定你的压力场能不能解出来。2.2 速度压力空间的匹配为什么Q1-Q1必然翻车我见过不止一个新手用四个节点的双线性单元同时插值速度和压力也就是 Q1-Q1组装出来的矩阵能解速度场看着也像回事但把压力场画出来就露馅了压力呈红蓝棋盘状交替跳动相邻单元的p值锯齿形错开物理上完全扯淡。原因在于 Q1-Q1 不满足 inf-sup 条件也叫 LBB 条件离散连续方程的约束矩阵 C 存在大量虚假零空间除了物理上允许的“整体常数压力”之外还多出了高频振荡模态。解决办法有两类。第一类是换插值空间给速度用二阶单元、压力保持一阶典型就是本文要用的 Q2-Q1Taylor-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 后牛顿”的切换策略。低雷诺数比如 Re100直接用牛顿、初值给全零通常也扛得住。雷诺数超过一千之后我一般先跑十几步 Picard 把流场轮廓暖起来再切牛顿把残差压到机器精度附近。这个切换是这套程序里性价比最高的一个操作后面第3章的迭代循环会直接按这个结构写。3. 用MATLAB把有限元求解步骤落成最小骨架3.1 结构化网格生成与自由度编号先处理几何和编号。为了不引入网格生成器的黑匣子直接在 [0,1]×[0,1] 方腔上生成结构化四边形网格。Q2 单元每个方向上是二次映射所以每个大单元内部沿边长有2个节点间隔节点总数是 (2nx1)×(2ny1)。下面这段代码生成节点坐标并把每个单元的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) (j02-1)*npx i0 2; % 右上 e(4) (j02-1)*npx i0; % 左上 % 边中点和中心 e(5) (j0-1)*npx i0 1; % 下边中点 e(6) (j01-1)*npx i0 2; % 右边中点 e(7) (j02-1)*npx i0 1; % 上边中点 e(8) (j01-1)*npx i0; % 左边中点 e(9) (j01-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尺寸为 (184)x(184)顺序先速度后压力 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, 18b) Ke(2*a-1, 18b) C_ab*gw(i)*gw(j)*detJ; Ke(2*a, 18b) Ke(2*a, 18b) 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_v1: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顶盖y1水平速度 ux1垂直速度 uy0其余三条边速度全为零压力不设置边界条件由连续方程隐式确定但整个系统的压力只定义到相差一个常数。不要小看这个“相差常数”它是程序里最隐蔽的雷。下表是推荐的无量纲参数设置参数取值说明腔体尺寸1×1无量纲化后不再关注物理尺寸顶盖速度ux1, uy0方腔的唯一边界驱动雷诺数Re100 起步1000 进阶ReUL/νν1/Re网格16×16 至 64×64先16跑通再加密看变化压力参考点某个固定节点置0消除常数压力自由度收敛容差1e-8残差范数稳态问题值得压到机器精度附近无量纲化值得多说一句。如果你把程序写成带单位的参数后面想做无量纲基准解对比时会非常痛苦。我一般在一开始就把长度除以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(end1, :) 0; A_in(:, end1) 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 松弛因子、迭代初值与雷诺数续化越是新手越容易忽视迭代初值总觉得从全零开始迭代天经地义。对Re100的层流方腔全零初值确实没问题但Re超过1000后全零初值等于把牛顿迭代扔进一个完全陌生的区域前几步残差可能正常下降随后突然爆成Inf。这个现象的根源是对流项雅可比矩阵在远离解的地方不具备收敛性所以工程上常用“续化法”从小到大依次解Re100、400、1000、2000把上一个雷诺数的收敛解当作下一个雷诺数的迭代初值。每步之间Re的增幅一般不超过2倍加密网格时也要同步降低Re增幅。松弛因子和续化法配合起来才是一套完整的起手式。Step1用alpha0.5跑PicardStep2切牛顿时alpha回到1.0Step3在一切顺利的前提下把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*Re2 左右不满足就先加密别急着调求解器。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)再过细又变了。原因这种情况多半是边界条件被施加到了错误的自由度上。例如顶盖驱动流里顶盖的节点集合不止包含顶部一条线的角点还包括上边中点如果你只把角点设为ux1而漏了边中点Q2单元的二次边在顶部就会产生速度波动网格加密一换波动形态全变。另一个原因是高斯积分阶不足N-S矩阵里的对流项含速度梯度被积函数阶次比弹性问题高3×3高斯积分只是起步有些单元用3点积分在变形网格上会掉精度。解决检查边界自由度的覆盖范围把每个边界单元的全部相关节点列出来核对速度值随后做一次网格收敛性研究用16、32、64三套网格看某个特征量比如中心涡强度是否单调逼近一个值。5.5 超过一定网格规模后MATLAB直接法内存崩溃现象网格从32×32加密到64×64程序跑起来后内存占用飙升MATLAB直接报“Out of memory”或者把系统swap填满。原因N-S问题的自由度增长速度是二维的速度自由度为2×(2n1)²量级压力自由度为n²量级总自由度数约9n²对应的稀疏矩阵带宽也随n增大。直接法 A\b 在求解过程中会产生大量填充元矩阵非零元数量可能比原矩阵多几十倍这是直接法的固有属性不是代码写错了。解决换迭代求解器比如GMRES预处理加ILU分解但这类鞍点问题需要有效的预处理子简单ILU在Re高时收敛也慢。折中方案是先用一个中等网格把整个流程验证完毕再做最终大规模网格计算另一个实用的做法是把压力自由度重新排序到速度自由度之后并采用嵌套剖析排序能显著减小填充元数量。如果对迭代法不熟不要硬上用直接法跑一个40×40左右的网格完成论文级验证是完全足够的。6. 验证程序正确性用方腔基准解与网格加密判断可靠性6.1 对照Ghia基准解看中心涡位置和速度剖面程序写好之后第一件事不是跑新问题而是跑标准算例。方腔驱动流的对照基准是Ghia在1982年给出的Re100、400、1000的数值解这组结果被引用了四十多年是这类程序真正的试金石。对照方法很简单取经过腔体几何中心的两条垂线/水平线输出u和v速度分量与基准表中的数据点对比。Re100时中心涡的坐标大约在(0.6,0.73)Re1000时主涡中心移到约(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]; % 举例对应nx16,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升高而增大——如果Re1000的误差比Re100大一个数量级以上基本说明网格分辨率不够不要归咎于求解器。6.3 一个养成很久的习惯先按特征量验收再信局部细节我给这类程序做验收时有个固定动作先记录中心涡位置、次级涡强度和速度剖面峰值这三个全局特征量再去看流线图云图。全局特征量对边界条件、压力参考点的错误不敏感但正是这种不敏感导致它们不是万能。所以下一步我会挑一个局部细节——比如顶盖边界下方第一层单元的速度梯度——单独输出数值检查这个量的错误率远高于全局量适合做程序bug的照妖镜。这十几年写有限元程序的经验归纳成一个习惯就是永远假设自己的编号和边界条件有一次手滑所以任何程序的第一版验收都以利用网格加密的对比结果为准而不是以单套网格的漂亮云图为准。这样即便某一步疏忽了也能通过收敛趋势把它抓出来。希望帮到你。本文还有配套的精品资源点击获取