简介面向流体力学数值计算方向的学习者与工程研究者这套 Matlab 程序以有限元方法求解二维稳态 Navier-Stokes 方程针对带复杂几何与非均匀边界条件的粘性流动给出速度场与压力场的数值结果适合具备偏微分方程与有限元初步基础的中高级读者研读。压缩包内共 45 个 .m 文件体积约 19KB全部为 Matlab 脚本与函数并按流程拆分为均匀网格生成、局部与参考基函数、高斯积分、矩阵及右端向量组装、Dirichlet 与应力、Robin 边界处理、Newton 迭代初值、压力单点固定、误差统计等模块。目前已有 315 人学习代码函数职责单一、命名对应清楚便于逐步跟踪离散、组装、求解与后处理全链路内置多组解析解与误差计算可用于检验收敛精度替换源项、边界条件或物性函数后还能迁移到热传导、扩散等同类椭圆型问题。1. 一个压力常数漂移引出的排查去年帮人看一份二维稳态 Navier-Stokes 方程有限元求解程序对方第一句话是「跑是能跑就是结果看着像随机生成的」。程序确实不长网格读取、稀疏组装、Matlab 的\求解都写全了速度场也画得出来问题出在压力场全场压力整体翘起来一个常数值出口附近的回流区被压成了正值看起来就像求解器在胡说。根因只有一行——压力泊松的零空间没有固定矩阵半正定直接解出来的压力只确定到一个常数。这类坑在这份程序里不止一处也正好说明它是个能读、能改、能拿来练手的完整实现而不是一段演示脚本。它做的事很明确把二维、稳态、不可压缩的 Navier-Stokes 方程用有限元离散解出速度场和压力场配套网格读取、边界条件施加、非线性迭代和后处理。适合三类人流体计算方向想从零搭一遍有限元的手里有 Matlab、想快速验证某个几何或 Re 数下流态的以及需要给后续的非稳态或湍流模型打底的。有限元、流体计算、Matlab 这三样凑在一起代码量可控改起来也不用重编译。2. 项目定位与技术选型它到底能解哪一类流动2.1 边界划清楚二维、稳态、不可压、层流为主程序求解的控制方程是定常不可压 Navier-Stokesρ (u·∇)u -∇p μ ∇²u f ∇·u 0无量纲化之后粘性项系数变成 1/Re对流项系数是 1整个问题的难度就压在 Re 上。稳态这一条限制最容易被忽略不是所有 Re 都有稳态解方腔这类几何在 Re 上千以后本身就进入非定常Newton 迭代再强也只会给你一个不存在的「定常解」或者干脆发散。经验上这套代码在 Re 一百到几百的量级最舒服跑方腔、后台阶、圆柱绕流、突扩管都合适Re 再高就得换伪瞬态延续或者直接上非稳态。不适用的场景也要说清楚免得改到一半才发现方向错了可压缩流密度不再常数连续性方程要换成质量守恒加状态方程、非牛顿流体粘度不再是标量常数、带自由面的流动网格要动、三维自由度直接翻几倍Matlab 的稀疏直接解扛不住。这些都不是加几行能解决的属于换一套求解框架的活。反过来它擅长的场景其实很具体教学复现、算法对比比如 Picard 和 Newton 的收敛半径对比、参数扫描批量跑不同 Re、不同几何尺寸、给代理模型造样本。这类任务要的是「跑得快、改得动、结果可信」而不是极致性能Matlab 加稀疏矩阵正好卡在这个位置上。2.2 有限元对上有限体积为什么这份程序用弱形式工程上算流体有限体积是主流局部守恒、通量好算、高 Re 下处理对流有成熟套路。这份程序偏偏用有限元理由不是「有限元更先进」而是几个很实际的点。第一是边界条件的处理。有限元的弱形式里Neumann 边界给定应力、给定压力是自然边界条件直接进右端项不用像有限体积那样单独构造边界通量。出口给「压力为零」这种条件在弱形式里就是一项边界积分删掉或者保留都行实现成本低。第二是精度。P2/P1 的 Taylor-Hood 单元速度是二次、压力是线性的同一个网格下速度的插值阶比一阶格式高一档。低 Re、小几何、要看清涡心位置和再附长度的时候这个差别挺明显。第三是压力-速度耦合的写法干净。有限元里压力天然是拉格朗日乘子混合单元速度 P2、压力 P1满足 LBB 条件不会出现棋盘格压力。等阶元P1/P1才需要额外的压力修正项。这份程序默认就走 Taylor-Hood省掉一堆调试。代价也有对流占优时 Galerkin 格式会震荡Re 高了必须加流线迎风项SUPG或者用更细的网格局部守恒性比有限体积弱算流量平衡时得用一致通量后处理再就是组装比有限体积繁琐单元矩阵要按高斯积分逐点累加。这些在代码里都能看到痕迹改的时候心里有数就行。2.3 依赖清单与运行环境技术栈非常轻基本只用 Matlab 自带的能力依赖用途是否必需Matlab R2018a 及以上主运行环境sparse、ilu、gmres都在基础包里必需sparse/accumarray全局矩阵组装与自由度映射必需ilugmres大规模时的迭代求解与预条件可选小规模用\就够Gmsh 或任意三角网格生成器生成.msh网格文件可选代码自带示例网格trisurf/patch/quiver/contour云图、流线、矢量场后处理必需画图不用 PDE Toolbox 也能跑网格解析是自己写的避免了对工具箱版本差异的依赖。这一点对二次开发挺关键你把它拷到另一台机器上只要 Matlab 版本不太老基本不会因为缺工具箱而报错。3. 目录结构与主流程从网格文件到速度压力场3.1 工程目录与文件职责代码组织按「网格 → 单元 → 组装 → 求解 → 后处理」拆一个文件干一件事不写巨型脚本路径职责main.m主流程读配置、读网格、迭代、出图cases/caseName/config.m算例配置Re、边界标签、迭代参数、输出路径src/mesh/read_msh.m解析 Gmsh 的.msh输出节点、单元、边界标签src/mesh/build_dof_map.m建立节点到全局自由度的映射src/element/shape_p2.m/shape_p1.m形函数及其对面积坐标的导数src/element/elem_ns.m单个单元的粘性、对流、压力耦合矩阵src/assembly/assemble_system.m全局残差与雅可比的稀疏组装src/bc/apply_bc.m强加 Dirichlet 条件、压力钉扎src/solver/solve_linear.m直接法或 GMRESILU 求解增量方程src/post/plot_fields.m速度幅值、压力、流线、涡量云图src/post/compute_forces.m壁面剪应力积分、压差、升阻力这种拆法有个好处想换求解器只动solve_linear.m想换单元只动src/element/主流程几乎不用改。第一次读的时候按main.m的顺序往下走就行。3.2 主脚本逐段拆解main.m是整份程序的骨架逻辑就是「初始化 → 迭代 → 收敛判断 → 后处理」% main.m —— 二维稳态 NS 有限元求解主流程 clear; clc; close all; addpath(genpath(src)); cfg setup_config(cavity); % 选算例cavity / step / cylinder mesh read_msh(cfg.meshFile); % 读网格得到节点坐标与单元连接 mesh tag_boundaries(mesh, cfg.bc); % 按几何位置给边界节点打标签 dof build_dof_map(mesh); % 建立 [u; v; p] 的全局编号 U init_field(mesh, dof, cfg); % 初值静止流体 零压力 U apply_initial_bc(U, mesh, dof, cfg); % 壁面速度先置成边界值 for k 1:cfg.maxIter [R, J] assemble_system(mesh, dof, U, cfg); % 残差 雅可比 res norm(R, inf); fprintf(iter %3d |R|_inf %.3e\n, k, res); if res cfg.tol, break; end % 收敛跳出 dU solve_linear(J, -R, cfg); % 解线性化方程 U U cfg.relax * dU; % 松弛前几步取 0.7 左右 if any(isnan(U)) % 发散的早期信号 error(迭代发散检查 Re、网格与边界条件); end end post_process(mesh, dof, U, cfg); % 云图、流线、壁面力几个参数值得单独说。cfg.relax是欠松弛因子Newton 在解附近收敛快但离解远的时候一步走满容易冲出去前 3 到 5 步取 0.5 到 0.8后面再放回 1.0是最省事的做法。cfg.tol一般取 1e-8 到 1e-10取太松会出现「残差降不下去但流场已经不怎么变」的假收敛取太紧则在机器精度附近空转。cfg.maxIter给 50 足够超过 50 步还没收敛基本是模型或网格有问题继续跑没意义。tag_boundaries这步容易被跳过但它决定了后面所有边界条件能不能正确落位。常见做法是按节点坐标判定x 1e-8归到左边界y 1-1e-8归到顶盖角点同时属于两条边时优先级要写死否则同一个自由度被赋两个不同的值矩阵会直接奇异。3.3 自由度编号与全局向量布局混合单元的节点数不一样速度和压力分开编号是最不容易出错的方式function dof build_dof_map(mesh) % 速度用 P26 节点三角形压力用 P13 节点三角形 nU mesh.nNodes; % 速度节点总数 nP mesh.nPressureNodes; % 压力节点总数P1 节点是 P2 的子集 dof.u (1:nU); % u 分量 dof.v nU (1:nU); % v 分量 dof.p 2*nU (1:nP); % 压力 dof.nTotal 2*nU nP; end全局向量就是U [u; v; p]一次线性求解把三个场同时更新。单元级的自由度收集靠一张索引表rows [ dof.u(mesh.elem(e).nodes); % 6 个速度节点 → u dof.v(mesh.elem(e).nodes); % 6 个速度节点 → v dof.p(mesh.elem(e).pnodes) ]; % 3 个压力节点这样拼出来的单元矩阵是 15×1512 个速度自由度加 3 个压力自由度全局矩阵是对称不定的鞍点结构写出来是分块形式J [ A B ] [ B 0 ]A是粘性加对流线性化后的速度块B是压力梯度B是散度约束。这个结构决定了后面线性求解器怎么选——正定对称的解法在这里都不适用。3.4 边界条件的三种施加方式Dirichlet 条件壁面速度、入口速度的施加方式直接影响矩阵性质代码里用的是消元法function [J, R] apply_dirichlet(J, R, dofs, values, U) % dofs : 被约束的自由度编号 % values : 对应的目标值 % U : 当前解向量用于计算增量 for k 1:numel(dofs) i dofs(k); delta values(k) - U(i); % 本步需要补上的增量 R R - J(:, i) * delta;% 已知增量挪到右端必须在清零之前算 R(i) delta; % 该行换成 u_i 目标值 J(i, :) 0; J(i, i) 1; J(:, i) 0; % 列也清零保持矩阵对称 end先算R - J(:,i)*delta再清零列顺序不能反反了就等于把已知位移漏掉了收敛会变成线性甚至不收敛。列清零是为了保持对称性规模大的时候对 ILU 预条件的质量影响很明显。另一种做法是罚函数法把对角元加大数实现简单但会污染条件数Re 高的时候不建议。压力场还差一步整个压力只有相对意义必须钉一个参考点否则矩阵奇异。% 压力钉扎把第 1 个压力自由度固定为 0 [J, R] apply_dirichlet(J, R, dof.p(1), 0, U);这一步就是第 1 章那个「压力整体漂移」的解法。钉哪个点都行出口中心、左下角都可以结果只差一个常数物理量压差、壁面剪应力不受影响。4. 核心模块拆解单元形函数、稀疏组装与非线性迭代4.1 Taylor-Hood 单元与形函数实现P2 三角形有 6 个节点3 个顶点加 3 个边中点。形函数用面积坐标写起来最简洁function [N, dNdL] shape_p2(L) % L [L1; L2; L3] 为面积坐标L1L2L3 1 % 节点顺序顶点 1,2,3边中点 12, 23, 31 N [ L(1)*(2*L(1)-1); L(2)*(2*L(2)-1); L(3)*(2*L(3)-1); 4*L(1)*L(2); 4*L(2)*L(3); 4*L(3)*L(1) ]; % 对面积坐标的偏导6×3 dNdL [ 4*L(1)-1, 0, 0; 0, 4*L(2)-1, 0; 0, 0, 4*L(3)-1; 4*L(2), 4*L(1), 0; 0, 4*L(3), 4*L(2); 4*L(3), 0, 4*L(1) ]; end压力用 P1形函数就是面积坐标本身N LdNdL eye(3)。两者配合满足 LBB 条件压力不会出现棋盘格。这条不是可选项如果图省事把压力和速度都用 P1等阶元散度约束太弱压力场会出现逐单元交替的正负跳变看着像网格噪声实质是单元不满足 inf-sup。真要用等阶元得额外加压力修正项代码量翻倍不如直接用 Taylor-Hood。从面积坐标到物理坐标的导数要走雅可比% Jm 为 (x,y) 对 (L1,L2) 的雅可比2×2 Jm [x2-x1, x3-x1; y2-y1, y3-y1]; detJ det(Jm); dNdx dNdL * inv(Jm); % 6×2第一列对 x第二列对 ydetJ的符号要检查一下。网格里如果有节点顺序反了的单元顺时针detJ为负积分出来的单元矩阵是负定的装到全局里就是一锅乱粥。常见做法是在读网格时统一翻转或者在elem_ns里加一句if detJ 0, error(...)兜底。4.2 单元矩阵到全局稀疏矩阵的组装单元级要算出三块速度块Ae12×12、压力耦合块Be12×3、右端项fe。核心是一个高斯积分循环function [Ae, Be, fe] elem_ns(coords, ue, ve, Re, ngauss) % coords: 6×2 速度节点坐标ue, ve: 当前单元上的速度自由度 Ae zeros(12,12); Be zeros(12,3); fe zeros(12,1); for g 1:ngauss [L, w] gauss_tri(g); % 面积坐标与权重 [N, dNdL] shape_p2(L); [Np, dNpL] shape_p1(L); [dNdx, dNpdx, detJ] map_grad(coords, dNdL, dNpL); if detJ 0, error(单元 %d 雅可比为负检查节点顺序, g); end ug N * ue; vg N * ve; % 高斯点速度 % 粘性项(1/Re) ∫ ∇N_i · ∇N_j Ae Ae (1/Re) * (dNdx * dNdx) * w * detJ; % 对流项∫ N_i (u·∇)N_j —— Picard 线性化u 取上一轮解 convj dNdx * [ug; vg]; % 6×1(u·∇)N_j Ae Ae (N * convj) * w * detJ; % 压力梯度耦合 B∫ N_i ∂Np_j/∂x Be(1:6, :) Be(1:6, :) - (N * dNpdx(:,1)) * w * detJ; Be(7:12, :) Be(7:12, :) - (N * dNpdx(:,2)) * w * detJ; endconvj dNdx * [ug; vg]这一步是对流项的关键把高斯点的速度和对流速度分开当前解算出来的ug, vg只用于线性化未知量还是N这一侧所以Ae对速度是线性的。这就是 Picard 迭代的数学本质也是它比 Newton 好实现的地方。组装到全局时有个性能坑。写法一最直观J(rows, rows) J(rows, rows) Ae; % 稀疏矩阵增量索引赋值Matlab 每次都要重新整理稀疏结构单元数上万以后会慢得离谱我的血泪经验是十几万自由度的模型能跑一晚上。正确做法是攒三元组最后一次性sparse% 预分配三元组单元数 × 15 × 15 的上界 nzmax mesh.nElem * 225; I zeros(nzmax,1); Jc zeros(nzmax,1); V zeros(nzmax,1); ptr 0; for e 1:mesh.nElem [Ae, Be, fe] elem_ns(...); Ke [Ae, Be; Be, zeros(3,3)]; % 15×15 单元矩阵 rows gather_dofs(mesh, dof, e); [ii, jj] ndgrid(rows, rows); idx ptr (1:numel(ii)); I(idx) ii(:); Jc(idx) jj(:); V(idx) Ke(:); ptr ptr numel(ii); end J sparse(I(1:ptr), Jc(1:ptr), V(1:ptr), dof.nTotal, dof.nTotal);同一位置被多个单元累加时sparse会自动求和不用自己写accumarray。这个改动能把组装时间压下去一个数量级值得一开始就写对。4.3 Picard 与 Newton迭代格式怎么选两种非线性迭代的差别在雅可比。Picard 用上一轮的速度线性化对流项雅可比里只有粘性加对流收敛是线性的每步残差大概降一个量级Newton 把对流项对速度的导数也算进去解附近二次收敛但收敛半径小初值差一点就发散。实际用的时候不建议二选一而是分段mode picard; for k 1:cfg.maxIter cfg.jacobian mode; % 组装时决定是否含对流导数项 [R, J] assemble_system(mesh, dof, U, cfg); res norm(R, inf); % 残差降到 1e-2 以下且已经走了足够步数切 Newton 加速 if strcmp(mode,picard) res cfg.switchTol k cfg.picardSteps mode newton; fprintf( - 切换到 Newton 迭代\n); end if res cfg.tol, break; end dU solve_linear(J, -R, cfg); U U cfg.relax * dU; endcfg.switchTol取 1e-2 到 1e-3cfg.picardSteps取 5 到 10。Re 三百以内这样跑通常 15 步内收敛Re 上千的时候 Picard 阶段就可能卡住这时候要上伪瞬态延续——把稳态问题当成dU/dt -R(U)的定常解加一个虚拟时间项改善条件数J J speye(dof.nTotal) / dt; % 简化处理严格写法应加质量矩阵 dU solve_linear(J, -R, cfg); dt min(dt * cfg.dtGrowth, cfg.dtMax); % 逐步放大时间步时间项一加矩阵的对角占优程度变好ILU 预条件也更给力。dt从小往大放比如 1e-3 起步每步乘 1.2上限 1e2是这类稳态求解器能不能跑到高 Re 的分水岭。4.4 线性求解直接法、GMRESILU 与内存权衡每步非线性迭代都要解一次J·dU -R这一步往往占总时间的七成以上。程序里留了两种模式function dU solve_linear(J, b, cfg) switch cfg.solver case direct dU J \ b; % 稀疏 LU小规模最省心 case gmres setup.type ilutp; % 带阈值丢弃的 ILU setup.droptol 1e-3; setup.udiag 1; [L, Umat] ilu(J, setup); [dU, flag, relres, iter] gmres(J, b, cfg.restart, ... cfg.linTol, cfg.linMaxIter, L, Umat); if flag ~ 0 warning(GMRES 未收敛relres %.2e, iter %d, relres, iter); end end end规模在十万自由度以内\基本够用而且不用调参数出问题也好排查。往上走稀疏 LU 的填充会吃掉大量内存二维问题还好三维直接爆。GMRES 的关键参数是restart重启步数20 到 50 之间和droptolILU 丢弃阈值1e-3 到 1e-4droptol放太松预条件不准迭代步数飙升放太紧 ILU 本身的内存又上去了得试。求解方式适用规模自由度内存调参成本稀疏直接解\十万以内中到高无需调参GMRES ILU(0)十万到百万低低收敛可能偏慢GMRES ILUT十万到百万中需要调droptol直接解 分块消元十万以内高需要理解鞍点结构鞍点结构还有个专门技巧对[A B; B 0]做块消元先解压力泊松BA⁻¹B p ...再回代速度。这样线性系统规模从2nUnP降到nP压力节点数比速度少一大截省下来的内存很可观。代价是要算A⁻¹B通常用A的 ILU 做近似实现复杂度上一个台阶属于二次开发阶段的优化项。5. 避坑与排查五类让结果「看着像随机数」的故障5.1 先跑最小算例三步验证法拿到代码别急着上自己的几何按这三步走一遍能提前暴露八成问题。第一步关掉对流项cfg.Re Inf或者直接把对流块置零退化成一维泊肃叶流。管道两端加压力差解析解是抛物线出口截面速度剖面可以直接对。P2 单元在L2范数下的收敛阶应该是三阶网格加密一次误差降八倍对不上就说明形函数或者雅可比有问题。第二步打开对流跑方腔Re 100。取竖直中心线上的u分量分布和公开的基准数据对比涡心位置误差在 1% 以内算通过。这一步能抓出对流项符号错、边界标签漏打、压力没钉扎这几类问题。第三步网格加密两次看再附长度或者涡心位置的相对变化。变化小于 2% 可以认为网格无关再加密就是浪费算力。5.2 五类高频故障的现象、原因与解决矩阵奇异\报 warning 或者解出 NaN。现象是第一次迭代就崩norm(R,inf)直接是Inf。原因通常是两件事之一压力参考点没固定鞍点矩阵的零空间没消掉或者某个边界节点被两个边界条件同时约束同一行被赋了两个不同的值。解决是把dof.p(1)钉成 0再检查tag_boundaries里角点的优先级是不是写死的。方腔的左下角同时属于左壁和下壁两个都是无滑移还好入口和壁面相交的角点就容易冲突。残差降到 1e-2 就再也不动迭代曲线拉平。现象是残差曲线像撞了墙。原因多半是线性求解精度不够——linTol设得比非线性tol还松外层根本收敛不了。内外层容差要匹配外层1e-8的话内层至少1e-10否则 Newton 的二次收敛性直接失效。另一个可能是 ILU 预条件太差droptol调小一档试试。速度场「漏」过壁面流线穿墙。现象很直观画流线的时候能看到线从壁面穿出去。原因是边界节点集不全——P2 单元的边中点节点是独立的自由度只按顶点坐标打标签会漏掉边中点。解决是在tag_boundaries里对边中点单独判定取两个端点的中点坐标再判或者干脆用几何距离判据所有到壁面距离小于1e-10的节点全打上标签。解对称、光滑但物理上明显不对。比如方腔的涡心跑到角落去了或者圆柱绕流没有脱涡稳态本来就没有但尾迹长度不对。这类最难查因为不报错。常见原因是雅可比和残差不一致——残差用了完整对流项雅可比用了个错的线性化Newton 会收敛到一个「雅可比为零」的点而不是残差为零的点。排查办法是数值验证雅可比随机取一个方向d比较(R(Uεd)-R(U-εd))/(2ε)和J·d相对误差应在1e-6量级。越算越慢内存一路涨。现象是前几步还行十步以后每次迭代要几分钟。原因在 4.2 提到的稀疏增量赋值每次都重建稀疏结构。改成三元组一次性sparse之后能快一个数量级。另一个隐性原因是每步都重新ilu如果矩阵变化不大可以隔几步重算一次预条件能省下不少时间。提示把这五类现象做成一张检查表贴在脚本旁边出问题按顺序过一遍比漫无目的地读代码快得多。5.3 收敛曲线与网格无关性检查跑通之后收敛历史要存下来别只 print 到屏幕。后面调参、对比不同 Re、判断是不是假收敛全靠这条曲线hist zeros(cfg.maxIter, 3); % [残差, 线性迭代数, 松弛因子] for k 1:cfg.maxIter [R, J] assemble_system(mesh, dof, U, cfg); res norm(R, inf); [dU, linIter] solve_linear(J, -R, cfg); U U cfg.relax * dU; hist(k,:) [res, linIter, cfg.relax]; if res cfg.tol, hist hist(1:k,:); break; end end semilogy(hist(:,1), o-); grid on; xlabel(迭代步); ylabel(|R|_{inf});健康的曲线是「前几步陡降Newton 阶段接近垂直下落」。如果残差每步只降 10% 左右说明一直在 Picard 阶段没切过去如果残差开始震荡多半是松弛因子太大或者网格太粗把relax降到 0.5 再试。网格无关性用一个简单循环就能做hList [0.04, 0.02, 0.01]; for i 1:numel(hList) mesh gen_mesh(hList(i), cfg); U solve_ns(mesh, cfg); metric(i) extract_recirculation_length(mesh, U); % 再附长度 end relChange abs(diff(metric)) ./ abs(metric(1:end-1));relChange连续两次小于 2%这个网格才敢拿来出结果。省掉这一步后面所有的物理解释都建立在流沙上。6. 二次开发换几何、换物性、换求解器拿到代码之后最常改的三件事按改动成本从低到高排一下。改造点需要动的地方风险换几何 / 换算例只改cases/name/config.m和网格文件低注意边界标签名要对上改 Re、改入口剖面config.m里的Re和inletProfile函数低Re 高了要配伪瞬态加 SUPG 迎风项elem_ns.m的对流块加流线方向权函数中tau参数要调换非牛顿粘度elem_ns.m里把1/Re换成mu(gamma_dot)/Re中要在 Picard 里更新粘度换线性求解器solve_linear.m中预条件参数要试改稳态为非稳态main.m加时间循环质量矩阵进组装高等于重写主流程换几何是最容易上手的用网格工具生成.msh在config.m里把边界标签映射写清楚主流程一行不用改。我一般会先在tag_boundaries里把每条边的节点数打印出来跟几何周长对一下数目明显不对就是标签判据写错了。改 Re 也简单但别一步从 100 跳到 2000中间没有延续步Newton 大概率发散。加 SUPG 是性价比最高的一档改造。对流占优时 Galerkin 格式的震荡不是数值噪声是格式本身的问题加一个流线方向的权函数就能压下去% SUPG 权函数W_i N_i tau * (u·∇N_i) h sqrt(2 * area_e); % 单元特征尺寸 tau h / (2 * norm([ug; vg])); % 简化取法够用 W N tau * convj; % convj (u·∇)N_j Ae Ae (W * convj) * w * detJ;tau的取法有很多版本这个h/(2|u|)是最粗的一档低 Re 下几乎不影响精度高 Re 下能把震荡压住。要注意它破坏了 Galerkin 的最优性网格足够细的时候反而是负担所以别默认开着按 Re 判断。验证方法我建议固定成三件套残差曲线、中心线剖面与基准数据的对比、网格加密后的相对变化。三样都过才敢说改对了。改非牛顿粘度的时候尤其要小心粘度依赖剪切率之后Picard 里必须同步更新mu忘了更新会收敛到一个「用旧粘度算出新速度」的不动点残差很小但物理是错的。从那以后我每次改完这类求解器都强制走一遍「泊肃叶 → 方腔 Re100 → 网格加密」这条线不管改动多小。省下的调试时间远比跑这三步多。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/2 5:31:59 · 阅读 1.2k Matlab二维稳态Navier-Stokes有限元求解器实现与避坑
简介面向流体力学数值计算方向的学习者与工程研究者这套 Matlab 程序以有限元方法求解二维稳态 Navier-Stokes 方程针对带复杂几何与非均匀边界条件的粘性流动给出速度场与压力场的数值结果适合具备偏微分方程与有限元初步基础的中高级读者研读。压缩包内共 45 个 .m 文件体积约 19KB全部为 Matlab 脚本与函数并按流程拆分为均匀网格生成、局部与参考基函数、高斯积分、矩阵及右端向量组装、Dirichlet 与应力、Robin 边界处理、Newton 迭代初值、压力单点固定、误差统计等模块。目前已有 315 人学习代码函数职责单一、命名对应清楚便于逐步跟踪离散、组装、求解与后处理全链路内置多组解析解与误差计算可用于检验收敛精度替换源项、边界条件或物性函数后还能迁移到热传导、扩散等同类椭圆型问题。1. 一个压力常数漂移引出的排查去年帮人看一份二维稳态 Navier-Stokes 方程有限元求解程序对方第一句话是「跑是能跑就是结果看着像随机生成的」。程序确实不长网格读取、稀疏组装、Matlab 的\求解都写全了速度场也画得出来问题出在压力场全场压力整体翘起来一个常数值出口附近的回流区被压成了正值看起来就像求解器在胡说。根因只有一行——压力泊松的零空间没有固定矩阵半正定直接解出来的压力只确定到一个常数。这类坑在这份程序里不止一处也正好说明它是个能读、能改、能拿来练手的完整实现而不是一段演示脚本。它做的事很明确把二维、稳态、不可压缩的 Navier-Stokes 方程用有限元离散解出速度场和压力场配套网格读取、边界条件施加、非线性迭代和后处理。适合三类人流体计算方向想从零搭一遍有限元的手里有 Matlab、想快速验证某个几何或 Re 数下流态的以及需要给后续的非稳态或湍流模型打底的。有限元、流体计算、Matlab 这三样凑在一起代码量可控改起来也不用重编译。2. 项目定位与技术选型它到底能解哪一类流动2.1 边界划清楚二维、稳态、不可压、层流为主程序求解的控制方程是定常不可压 Navier-Stokesρ (u·∇)u -∇p μ ∇²u f ∇·u 0无量纲化之后粘性项系数变成 1/Re对流项系数是 1整个问题的难度就压在 Re 上。稳态这一条限制最容易被忽略不是所有 Re 都有稳态解方腔这类几何在 Re 上千以后本身就进入非定常Newton 迭代再强也只会给你一个不存在的「定常解」或者干脆发散。经验上这套代码在 Re 一百到几百的量级最舒服跑方腔、后台阶、圆柱绕流、突扩管都合适Re 再高就得换伪瞬态延续或者直接上非稳态。不适用的场景也要说清楚免得改到一半才发现方向错了可压缩流密度不再常数连续性方程要换成质量守恒加状态方程、非牛顿流体粘度不再是标量常数、带自由面的流动网格要动、三维自由度直接翻几倍Matlab 的稀疏直接解扛不住。这些都不是加几行能解决的属于换一套求解框架的活。反过来它擅长的场景其实很具体教学复现、算法对比比如 Picard 和 Newton 的收敛半径对比、参数扫描批量跑不同 Re、不同几何尺寸、给代理模型造样本。这类任务要的是「跑得快、改得动、结果可信」而不是极致性能Matlab 加稀疏矩阵正好卡在这个位置上。2.2 有限元对上有限体积为什么这份程序用弱形式工程上算流体有限体积是主流局部守恒、通量好算、高 Re 下处理对流有成熟套路。这份程序偏偏用有限元理由不是「有限元更先进」而是几个很实际的点。第一是边界条件的处理。有限元的弱形式里Neumann 边界给定应力、给定压力是自然边界条件直接进右端项不用像有限体积那样单独构造边界通量。出口给「压力为零」这种条件在弱形式里就是一项边界积分删掉或者保留都行实现成本低。第二是精度。P2/P1 的 Taylor-Hood 单元速度是二次、压力是线性的同一个网格下速度的插值阶比一阶格式高一档。低 Re、小几何、要看清涡心位置和再附长度的时候这个差别挺明显。第三是压力-速度耦合的写法干净。有限元里压力天然是拉格朗日乘子混合单元速度 P2、压力 P1满足 LBB 条件不会出现棋盘格压力。等阶元P1/P1才需要额外的压力修正项。这份程序默认就走 Taylor-Hood省掉一堆调试。代价也有对流占优时 Galerkin 格式会震荡Re 高了必须加流线迎风项SUPG或者用更细的网格局部守恒性比有限体积弱算流量平衡时得用一致通量后处理再就是组装比有限体积繁琐单元矩阵要按高斯积分逐点累加。这些在代码里都能看到痕迹改的时候心里有数就行。2.3 依赖清单与运行环境技术栈非常轻基本只用 Matlab 自带的能力依赖用途是否必需Matlab R2018a 及以上主运行环境sparse、ilu、gmres都在基础包里必需sparse/accumarray全局矩阵组装与自由度映射必需ilugmres大规模时的迭代求解与预条件可选小规模用\就够Gmsh 或任意三角网格生成器生成.msh网格文件可选代码自带示例网格trisurf/patch/quiver/contour云图、流线、矢量场后处理必需画图不用 PDE Toolbox 也能跑网格解析是自己写的避免了对工具箱版本差异的依赖。这一点对二次开发挺关键你把它拷到另一台机器上只要 Matlab 版本不太老基本不会因为缺工具箱而报错。3. 目录结构与主流程从网格文件到速度压力场3.1 工程目录与文件职责代码组织按「网格 → 单元 → 组装 → 求解 → 后处理」拆一个文件干一件事不写巨型脚本路径职责main.m主流程读配置、读网格、迭代、出图cases/caseName/config.m算例配置Re、边界标签、迭代参数、输出路径src/mesh/read_msh.m解析 Gmsh 的.msh输出节点、单元、边界标签src/mesh/build_dof_map.m建立节点到全局自由度的映射src/element/shape_p2.m/shape_p1.m形函数及其对面积坐标的导数src/element/elem_ns.m单个单元的粘性、对流、压力耦合矩阵src/assembly/assemble_system.m全局残差与雅可比的稀疏组装src/bc/apply_bc.m强加 Dirichlet 条件、压力钉扎src/solver/solve_linear.m直接法或 GMRESILU 求解增量方程src/post/plot_fields.m速度幅值、压力、流线、涡量云图src/post/compute_forces.m壁面剪应力积分、压差、升阻力这种拆法有个好处想换求解器只动solve_linear.m想换单元只动src/element/主流程几乎不用改。第一次读的时候按main.m的顺序往下走就行。3.2 主脚本逐段拆解main.m是整份程序的骨架逻辑就是「初始化 → 迭代 → 收敛判断 → 后处理」% main.m —— 二维稳态 NS 有限元求解主流程 clear; clc; close all; addpath(genpath(src)); cfg setup_config(cavity); % 选算例cavity / step / cylinder mesh read_msh(cfg.meshFile); % 读网格得到节点坐标与单元连接 mesh tag_boundaries(mesh, cfg.bc); % 按几何位置给边界节点打标签 dof build_dof_map(mesh); % 建立 [u; v; p] 的全局编号 U init_field(mesh, dof, cfg); % 初值静止流体 零压力 U apply_initial_bc(U, mesh, dof, cfg); % 壁面速度先置成边界值 for k 1:cfg.maxIter [R, J] assemble_system(mesh, dof, U, cfg); % 残差 雅可比 res norm(R, inf); fprintf(iter %3d |R|_inf %.3e\n, k, res); if res cfg.tol, break; end % 收敛跳出 dU solve_linear(J, -R, cfg); % 解线性化方程 U U cfg.relax * dU; % 松弛前几步取 0.7 左右 if any(isnan(U)) % 发散的早期信号 error(迭代发散检查 Re、网格与边界条件); end end post_process(mesh, dof, U, cfg); % 云图、流线、壁面力几个参数值得单独说。cfg.relax是欠松弛因子Newton 在解附近收敛快但离解远的时候一步走满容易冲出去前 3 到 5 步取 0.5 到 0.8后面再放回 1.0是最省事的做法。cfg.tol一般取 1e-8 到 1e-10取太松会出现「残差降不下去但流场已经不怎么变」的假收敛取太紧则在机器精度附近空转。cfg.maxIter给 50 足够超过 50 步还没收敛基本是模型或网格有问题继续跑没意义。tag_boundaries这步容易被跳过但它决定了后面所有边界条件能不能正确落位。常见做法是按节点坐标判定x 1e-8归到左边界y 1-1e-8归到顶盖角点同时属于两条边时优先级要写死否则同一个自由度被赋两个不同的值矩阵会直接奇异。3.3 自由度编号与全局向量布局混合单元的节点数不一样速度和压力分开编号是最不容易出错的方式function dof build_dof_map(mesh) % 速度用 P26 节点三角形压力用 P13 节点三角形 nU mesh.nNodes; % 速度节点总数 nP mesh.nPressureNodes; % 压力节点总数P1 节点是 P2 的子集 dof.u (1:nU); % u 分量 dof.v nU (1:nU); % v 分量 dof.p 2*nU (1:nP); % 压力 dof.nTotal 2*nU nP; end全局向量就是U [u; v; p]一次线性求解把三个场同时更新。单元级的自由度收集靠一张索引表rows [ dof.u(mesh.elem(e).nodes); % 6 个速度节点 → u dof.v(mesh.elem(e).nodes); % 6 个速度节点 → v dof.p(mesh.elem(e).pnodes) ]; % 3 个压力节点这样拼出来的单元矩阵是 15×1512 个速度自由度加 3 个压力自由度全局矩阵是对称不定的鞍点结构写出来是分块形式J [ A B ] [ B 0 ]A是粘性加对流线性化后的速度块B是压力梯度B是散度约束。这个结构决定了后面线性求解器怎么选——正定对称的解法在这里都不适用。3.4 边界条件的三种施加方式Dirichlet 条件壁面速度、入口速度的施加方式直接影响矩阵性质代码里用的是消元法function [J, R] apply_dirichlet(J, R, dofs, values, U) % dofs : 被约束的自由度编号 % values : 对应的目标值 % U : 当前解向量用于计算增量 for k 1:numel(dofs) i dofs(k); delta values(k) - U(i); % 本步需要补上的增量 R R - J(:, i) * delta;% 已知增量挪到右端必须在清零之前算 R(i) delta; % 该行换成 u_i 目标值 J(i, :) 0; J(i, i) 1; J(:, i) 0; % 列也清零保持矩阵对称 end先算R - J(:,i)*delta再清零列顺序不能反反了就等于把已知位移漏掉了收敛会变成线性甚至不收敛。列清零是为了保持对称性规模大的时候对 ILU 预条件的质量影响很明显。另一种做法是罚函数法把对角元加大数实现简单但会污染条件数Re 高的时候不建议。压力场还差一步整个压力只有相对意义必须钉一个参考点否则矩阵奇异。% 压力钉扎把第 1 个压力自由度固定为 0 [J, R] apply_dirichlet(J, R, dof.p(1), 0, U);这一步就是第 1 章那个「压力整体漂移」的解法。钉哪个点都行出口中心、左下角都可以结果只差一个常数物理量压差、壁面剪应力不受影响。4. 核心模块拆解单元形函数、稀疏组装与非线性迭代4.1 Taylor-Hood 单元与形函数实现P2 三角形有 6 个节点3 个顶点加 3 个边中点。形函数用面积坐标写起来最简洁function [N, dNdL] shape_p2(L) % L [L1; L2; L3] 为面积坐标L1L2L3 1 % 节点顺序顶点 1,2,3边中点 12, 23, 31 N [ L(1)*(2*L(1)-1); L(2)*(2*L(2)-1); L(3)*(2*L(3)-1); 4*L(1)*L(2); 4*L(2)*L(3); 4*L(3)*L(1) ]; % 对面积坐标的偏导6×3 dNdL [ 4*L(1)-1, 0, 0; 0, 4*L(2)-1, 0; 0, 0, 4*L(3)-1; 4*L(2), 4*L(1), 0; 0, 4*L(3), 4*L(2); 4*L(3), 0, 4*L(1) ]; end压力用 P1形函数就是面积坐标本身N LdNdL eye(3)。两者配合满足 LBB 条件压力不会出现棋盘格。这条不是可选项如果图省事把压力和速度都用 P1等阶元散度约束太弱压力场会出现逐单元交替的正负跳变看着像网格噪声实质是单元不满足 inf-sup。真要用等阶元得额外加压力修正项代码量翻倍不如直接用 Taylor-Hood。从面积坐标到物理坐标的导数要走雅可比% Jm 为 (x,y) 对 (L1,L2) 的雅可比2×2 Jm [x2-x1, x3-x1; y2-y1, y3-y1]; detJ det(Jm); dNdx dNdL * inv(Jm); % 6×2第一列对 x第二列对 ydetJ的符号要检查一下。网格里如果有节点顺序反了的单元顺时针detJ为负积分出来的单元矩阵是负定的装到全局里就是一锅乱粥。常见做法是在读网格时统一翻转或者在elem_ns里加一句if detJ 0, error(...)兜底。4.2 单元矩阵到全局稀疏矩阵的组装单元级要算出三块速度块Ae12×12、压力耦合块Be12×3、右端项fe。核心是一个高斯积分循环function [Ae, Be, fe] elem_ns(coords, ue, ve, Re, ngauss) % coords: 6×2 速度节点坐标ue, ve: 当前单元上的速度自由度 Ae zeros(12,12); Be zeros(12,3); fe zeros(12,1); for g 1:ngauss [L, w] gauss_tri(g); % 面积坐标与权重 [N, dNdL] shape_p2(L); [Np, dNpL] shape_p1(L); [dNdx, dNpdx, detJ] map_grad(coords, dNdL, dNpL); if detJ 0, error(单元 %d 雅可比为负检查节点顺序, g); end ug N * ue; vg N * ve; % 高斯点速度 % 粘性项(1/Re) ∫ ∇N_i · ∇N_j Ae Ae (1/Re) * (dNdx * dNdx) * w * detJ; % 对流项∫ N_i (u·∇)N_j —— Picard 线性化u 取上一轮解 convj dNdx * [ug; vg]; % 6×1(u·∇)N_j Ae Ae (N * convj) * w * detJ; % 压力梯度耦合 B∫ N_i ∂Np_j/∂x Be(1:6, :) Be(1:6, :) - (N * dNpdx(:,1)) * w * detJ; Be(7:12, :) Be(7:12, :) - (N * dNpdx(:,2)) * w * detJ; endconvj dNdx * [ug; vg]这一步是对流项的关键把高斯点的速度和对流速度分开当前解算出来的ug, vg只用于线性化未知量还是N这一侧所以Ae对速度是线性的。这就是 Picard 迭代的数学本质也是它比 Newton 好实现的地方。组装到全局时有个性能坑。写法一最直观J(rows, rows) J(rows, rows) Ae; % 稀疏矩阵增量索引赋值Matlab 每次都要重新整理稀疏结构单元数上万以后会慢得离谱我的血泪经验是十几万自由度的模型能跑一晚上。正确做法是攒三元组最后一次性sparse% 预分配三元组单元数 × 15 × 15 的上界 nzmax mesh.nElem * 225; I zeros(nzmax,1); Jc zeros(nzmax,1); V zeros(nzmax,1); ptr 0; for e 1:mesh.nElem [Ae, Be, fe] elem_ns(...); Ke [Ae, Be; Be, zeros(3,3)]; % 15×15 单元矩阵 rows gather_dofs(mesh, dof, e); [ii, jj] ndgrid(rows, rows); idx ptr (1:numel(ii)); I(idx) ii(:); Jc(idx) jj(:); V(idx) Ke(:); ptr ptr numel(ii); end J sparse(I(1:ptr), Jc(1:ptr), V(1:ptr), dof.nTotal, dof.nTotal);同一位置被多个单元累加时sparse会自动求和不用自己写accumarray。这个改动能把组装时间压下去一个数量级值得一开始就写对。4.3 Picard 与 Newton迭代格式怎么选两种非线性迭代的差别在雅可比。Picard 用上一轮的速度线性化对流项雅可比里只有粘性加对流收敛是线性的每步残差大概降一个量级Newton 把对流项对速度的导数也算进去解附近二次收敛但收敛半径小初值差一点就发散。实际用的时候不建议二选一而是分段mode picard; for k 1:cfg.maxIter cfg.jacobian mode; % 组装时决定是否含对流导数项 [R, J] assemble_system(mesh, dof, U, cfg); res norm(R, inf); % 残差降到 1e-2 以下且已经走了足够步数切 Newton 加速 if strcmp(mode,picard) res cfg.switchTol k cfg.picardSteps mode newton; fprintf( - 切换到 Newton 迭代\n); end if res cfg.tol, break; end dU solve_linear(J, -R, cfg); U U cfg.relax * dU; endcfg.switchTol取 1e-2 到 1e-3cfg.picardSteps取 5 到 10。Re 三百以内这样跑通常 15 步内收敛Re 上千的时候 Picard 阶段就可能卡住这时候要上伪瞬态延续——把稳态问题当成dU/dt -R(U)的定常解加一个虚拟时间项改善条件数J J speye(dof.nTotal) / dt; % 简化处理严格写法应加质量矩阵 dU solve_linear(J, -R, cfg); dt min(dt * cfg.dtGrowth, cfg.dtMax); % 逐步放大时间步时间项一加矩阵的对角占优程度变好ILU 预条件也更给力。dt从小往大放比如 1e-3 起步每步乘 1.2上限 1e2是这类稳态求解器能不能跑到高 Re 的分水岭。4.4 线性求解直接法、GMRESILU 与内存权衡每步非线性迭代都要解一次J·dU -R这一步往往占总时间的七成以上。程序里留了两种模式function dU solve_linear(J, b, cfg) switch cfg.solver case direct dU J \ b; % 稀疏 LU小规模最省心 case gmres setup.type ilutp; % 带阈值丢弃的 ILU setup.droptol 1e-3; setup.udiag 1; [L, Umat] ilu(J, setup); [dU, flag, relres, iter] gmres(J, b, cfg.restart, ... cfg.linTol, cfg.linMaxIter, L, Umat); if flag ~ 0 warning(GMRES 未收敛relres %.2e, iter %d, relres, iter); end end end规模在十万自由度以内\基本够用而且不用调参数出问题也好排查。往上走稀疏 LU 的填充会吃掉大量内存二维问题还好三维直接爆。GMRES 的关键参数是restart重启步数20 到 50 之间和droptolILU 丢弃阈值1e-3 到 1e-4droptol放太松预条件不准迭代步数飙升放太紧 ILU 本身的内存又上去了得试。求解方式适用规模自由度内存调参成本稀疏直接解\十万以内中到高无需调参GMRES ILU(0)十万到百万低低收敛可能偏慢GMRES ILUT十万到百万中需要调droptol直接解 分块消元十万以内高需要理解鞍点结构鞍点结构还有个专门技巧对[A B; B 0]做块消元先解压力泊松BA⁻¹B p ...再回代速度。这样线性系统规模从2nUnP降到nP压力节点数比速度少一大截省下来的内存很可观。代价是要算A⁻¹B通常用A的 ILU 做近似实现复杂度上一个台阶属于二次开发阶段的优化项。5. 避坑与排查五类让结果「看着像随机数」的故障5.1 先跑最小算例三步验证法拿到代码别急着上自己的几何按这三步走一遍能提前暴露八成问题。第一步关掉对流项cfg.Re Inf或者直接把对流块置零退化成一维泊肃叶流。管道两端加压力差解析解是抛物线出口截面速度剖面可以直接对。P2 单元在L2范数下的收敛阶应该是三阶网格加密一次误差降八倍对不上就说明形函数或者雅可比有问题。第二步打开对流跑方腔Re 100。取竖直中心线上的u分量分布和公开的基准数据对比涡心位置误差在 1% 以内算通过。这一步能抓出对流项符号错、边界标签漏打、压力没钉扎这几类问题。第三步网格加密两次看再附长度或者涡心位置的相对变化。变化小于 2% 可以认为网格无关再加密就是浪费算力。5.2 五类高频故障的现象、原因与解决矩阵奇异\报 warning 或者解出 NaN。现象是第一次迭代就崩norm(R,inf)直接是Inf。原因通常是两件事之一压力参考点没固定鞍点矩阵的零空间没消掉或者某个边界节点被两个边界条件同时约束同一行被赋了两个不同的值。解决是把dof.p(1)钉成 0再检查tag_boundaries里角点的优先级是不是写死的。方腔的左下角同时属于左壁和下壁两个都是无滑移还好入口和壁面相交的角点就容易冲突。残差降到 1e-2 就再也不动迭代曲线拉平。现象是残差曲线像撞了墙。原因多半是线性求解精度不够——linTol设得比非线性tol还松外层根本收敛不了。内外层容差要匹配外层1e-8的话内层至少1e-10否则 Newton 的二次收敛性直接失效。另一个可能是 ILU 预条件太差droptol调小一档试试。速度场「漏」过壁面流线穿墙。现象很直观画流线的时候能看到线从壁面穿出去。原因是边界节点集不全——P2 单元的边中点节点是独立的自由度只按顶点坐标打标签会漏掉边中点。解决是在tag_boundaries里对边中点单独判定取两个端点的中点坐标再判或者干脆用几何距离判据所有到壁面距离小于1e-10的节点全打上标签。解对称、光滑但物理上明显不对。比如方腔的涡心跑到角落去了或者圆柱绕流没有脱涡稳态本来就没有但尾迹长度不对。这类最难查因为不报错。常见原因是雅可比和残差不一致——残差用了完整对流项雅可比用了个错的线性化Newton 会收敛到一个「雅可比为零」的点而不是残差为零的点。排查办法是数值验证雅可比随机取一个方向d比较(R(Uεd)-R(U-εd))/(2ε)和J·d相对误差应在1e-6量级。越算越慢内存一路涨。现象是前几步还行十步以后每次迭代要几分钟。原因在 4.2 提到的稀疏增量赋值每次都重建稀疏结构。改成三元组一次性sparse之后能快一个数量级。另一个隐性原因是每步都重新ilu如果矩阵变化不大可以隔几步重算一次预条件能省下不少时间。提示把这五类现象做成一张检查表贴在脚本旁边出问题按顺序过一遍比漫无目的地读代码快得多。5.3 收敛曲线与网格无关性检查跑通之后收敛历史要存下来别只 print 到屏幕。后面调参、对比不同 Re、判断是不是假收敛全靠这条曲线hist zeros(cfg.maxIter, 3); % [残差, 线性迭代数, 松弛因子] for k 1:cfg.maxIter [R, J] assemble_system(mesh, dof, U, cfg); res norm(R, inf); [dU, linIter] solve_linear(J, -R, cfg); U U cfg.relax * dU; hist(k,:) [res, linIter, cfg.relax]; if res cfg.tol, hist hist(1:k,:); break; end end semilogy(hist(:,1), o-); grid on; xlabel(迭代步); ylabel(|R|_{inf});健康的曲线是「前几步陡降Newton 阶段接近垂直下落」。如果残差每步只降 10% 左右说明一直在 Picard 阶段没切过去如果残差开始震荡多半是松弛因子太大或者网格太粗把relax降到 0.5 再试。网格无关性用一个简单循环就能做hList [0.04, 0.02, 0.01]; for i 1:numel(hList) mesh gen_mesh(hList(i), cfg); U solve_ns(mesh, cfg); metric(i) extract_recirculation_length(mesh, U); % 再附长度 end relChange abs(diff(metric)) ./ abs(metric(1:end-1));relChange连续两次小于 2%这个网格才敢拿来出结果。省掉这一步后面所有的物理解释都建立在流沙上。6. 二次开发换几何、换物性、换求解器拿到代码之后最常改的三件事按改动成本从低到高排一下。改造点需要动的地方风险换几何 / 换算例只改cases/name/config.m和网格文件低注意边界标签名要对上改 Re、改入口剖面config.m里的Re和inletProfile函数低Re 高了要配伪瞬态加 SUPG 迎风项elem_ns.m的对流块加流线方向权函数中tau参数要调换非牛顿粘度elem_ns.m里把1/Re换成mu(gamma_dot)/Re中要在 Picard 里更新粘度换线性求解器solve_linear.m中预条件参数要试改稳态为非稳态main.m加时间循环质量矩阵进组装高等于重写主流程换几何是最容易上手的用网格工具生成.msh在config.m里把边界标签映射写清楚主流程一行不用改。我一般会先在tag_boundaries里把每条边的节点数打印出来跟几何周长对一下数目明显不对就是标签判据写错了。改 Re 也简单但别一步从 100 跳到 2000中间没有延续步Newton 大概率发散。加 SUPG 是性价比最高的一档改造。对流占优时 Galerkin 格式的震荡不是数值噪声是格式本身的问题加一个流线方向的权函数就能压下去% SUPG 权函数W_i N_i tau * (u·∇N_i) h sqrt(2 * area_e); % 单元特征尺寸 tau h / (2 * norm([ug; vg])); % 简化取法够用 W N tau * convj; % convj (u·∇)N_j Ae Ae (W * convj) * w * detJ;tau的取法有很多版本这个h/(2|u|)是最粗的一档低 Re 下几乎不影响精度高 Re 下能把震荡压住。要注意它破坏了 Galerkin 的最优性网格足够细的时候反而是负担所以别默认开着按 Re 判断。验证方法我建议固定成三件套残差曲线、中心线剖面与基准数据的对比、网格加密后的相对变化。三样都过才敢说改对了。改非牛顿粘度的时候尤其要小心粘度依赖剪切率之后Picard 里必须同步更新mu忘了更新会收敛到一个「用旧粘度算出新速度」的不动点残差很小但物理是错的。从那以后我每次改完这类求解器都强制走一遍「泊肃叶 → 方腔 Re100 → 网格加密」这条线不管改动多小。省下的调试时间远比跑这三步多。希望帮到你。本文还有配套的精品资源点击获取