二维稳态Navier-Stokes方程有限元程序：从零实现与避坑指南

简介：这份资源是面向流体力学、数值计算方向的学习者与工程师的二维稳态Navier-Stokes方程有限元求解程序，基于Matlab实现，适合已具备偏微分方程与有限元基础、希望深入理解FEM求解流程的中高级读者。程序围绕动量方程与连续性方程展开，涵盖几何离散、函数空间定义、变分形式推导、矩阵组装、线性系统求解及速度压力后处理等完整环节，并涉及Dirichlet、Robin、应力边界处理与Newton迭代初值生成等细节。压缩包共45个文件，全部为m脚本，约19KB，按功能可大致分为网格与有限元空间生成、局部与参考基函数、单元与边界积分、矩阵与向量组装、边界条件处理、解析解与误差计算等模块，结构紧凑便于逐层阅读。目前已有312人学习下载。通过研读源码，读者可掌握有限元法求解流体问题的组织方式，并较容易地将其迁移到热传导、扩散等同类偏微分方程问题中。

1. 从“算不对”到“算得稳”：二维稳态Navier-Stokes方程有限元程序到底能解决什么

做流体计算的人多半有过这种经历：网格画好了，边界条件也设了，求解器一跑，残差曲线要么平着不动，要么直接飙到发散。尤其是碰到二维稳态Navier-Stokes方程，非线性项和压力速度耦合搅在一起，新手很容易在“这程序到底能不能用”上反复翻车。这份二维稳态Navier-Stokes方程有限元求解程序，就是冲着这个场景来的——它用有限元方法把稳态不可压N-S方程离散成可迭代求解的代数系统，配套Matlab脚本，能直接跑出方腔驱动流、管道流这类经典算例的速度场和压力场。适合正在做流体仿真课程设计、想从零理解有限元组装过程、或者需要一个可改可调的轻量级求解器来验证自己网格和边界条件的人。它不是工业级CFD软件，但胜在代码透明、参数可查、每一步都能自己动手改。

2. 有限元离散的底层逻辑：为什么稳态N-S方程不能直接当线性问题解

2.1 从强形式到弱形式：加权残差到底在干什么

二维稳态不可压N-S方程写成向量形式是两项：动量方程和对流项、扩散项、压力梯度项耦合，加上连续性方程约束速度散度为零。强形式要求解在每个点上满足微分方程，但有限元走的是弱形式路线——把方程两边乘一个测试函数，在区域上积分，再用分部积分把二阶导数降阶。这一步的意义在于：它把对解的光滑性要求从二阶降到一阶，允许用分片多项式去逼近。

具体到二维情况，速度分量u、v和压力p是三个未知场。常见做法是用Taylor-Hood单元，速度用二次插值、压力用一次插值，这样满足LBB条件，避免压力出现棋盘格振荡。我一般会先确认单元类型再动手组装，因为一旦速度压力同阶，后面迭代再稳也白搭。弱形式里对流项是非线性的，因为速度既是被求解量又出现在对流系数里，这就是为什么整个系统不能一次性线性求解，必须走Newton迭代或者Picard迭代。

2.2 Newton迭代与Picard迭代的选型差异

Picard迭代把对流速度用上一步的解固定住，每次解一个线性化问题，收敛慢但稳健，适合初场离真解不远的情况。Newton迭代把非线性残差对未知量求导，得到Jacobian矩阵，二次收敛但初值敏感。这份程序里两种迭代都留了接口，默认用Picard起步，跑几十步后再切Newton加速。参数上，Picard的松弛因子一般取0.7到1.0，Newton的阻尼系数在残差下降不单调时降到0.5以下。

% 非线性迭代主循环（Picard + Newton混合策略） max_iter = 50; % 最大迭代步数 tol = 1e-8; % 残差收敛容差 relax = 0.8; % Picard松弛因子 for iter = 1:max_iter % 用上一步速度组装对流项系数 [K, F] = assemble_system(U_prev, V_prev, mesh, bc); % 求解线性化后的代数方程组 dU = K \ F; % 松弛更新，避免初值差时振荡 U_new = U_prev + relax * dU; % 计算残差范数 res = norm(F - K * U_new) / norm(F); if res < tol fprintf('收敛于第 %d 步，残差 %.2e\n', iter, res); break; end U_prev = U_new; end

这段代码里assemble_system负责把弱形式里的积分算成矩阵和右端项，K \ F是Matlab的直接求解。松弛因子relax是关键：取1.0就是标准Picard，取太小收敛慢，取太大容易震荡。残差用相对范数判断，比绝对范数更适应不同量级的问题。实际跑的时候，如果50步还没降到1e-8，先别急着加步数，检查网格质量和边界条件有没有矛盾。

2.3 压力速度耦合的处理方式

稳态N-S方程里压力没有自己的演化方程，只以梯度形式出现在动量方程中，这导致离散后的矩阵对角线上压力块是零。常见处理有两种：一是用分块求解策略，先解动量方程预测速度，再用压力泊松方程修正；二是直接用增广拉格朗日或罚函数法把压力消去。这份程序用的是分块迭代，好处是每一步物理意义明确，坏处是压力泊松方程的条件数通常比动量方程差，迭代次数多了容易积累误差。

提示：如果压力场出现高频振荡，优先检查单元类型是否满足LBB条件，而不是去调松弛因子。

3. 从零跑通一个方腔驱动流算例：网格、边界与求解参数怎么设

3.1 网格生成与单元编号

方腔驱动流是验证N-S求解器的经典算例：正方形区域，顶盖以单位速度拖动，其余三边无滑移。网格用结构化四边形单元最省事，Matlab里可以用meshgrid生成节点坐标，再按行优先编号。单元连接关系用四个节点索引表示，注意节点编号要逆时针排列，否则雅可比行列式为负，组装出来的刚度矩阵会出问题。

% 生成N×N结构化四边形网格 N = 32; % 每边单元数 x = linspace(0, 1, N+1); % x方向节点坐标 y = linspace(0, 1, N+1); % y方向节点坐标 [X, Y] = meshgrid(x, y); nodes = [X(:), Y(:)]; % 节点坐标矩阵，每行一个节点 % 单元连接：每个四边形四个节点逆时针排列 elems = zeros(N*N, 4); for i = 1:N for j = 1:N n1 = (i-1)*(N+1) + j; n2 = n1 + 1; n3 = n2 + (N+1); n4 = n1 + (N+1); elems((i-1)*N+j, :) = [n1, n2, n3, n4]; end end

N取32时总节点数1089，速度自由度2178，压力自由度1089，直接求解还能接受。如果N加到64，矩阵规模翻四倍，Matlab的\求解器开始吃内存，这时候要考虑迭代求解器或者稀疏矩阵存储。节点编号顺序影响带宽，行优先编号的带宽大约是N，列优先会更大，所以别随便换。

3.2 边界条件的施加方式

方腔驱动流的边界条件分两类：顶盖u=1、v=0，其余三边u=0、v=0。在有限元里施加Dirichlet边界条件有两种做法：一是把边界节点对应的自由度从方程组里消去，修改右端项；二是用罚函数法在大对角线上加大数。第一种精度高但改矩阵结构麻烦，第二种简单但罚数取不好会污染解。这份程序用的是消去法，先组装完整矩阵，再对边界自由度做行列置零和对角置一。

% 施加Dirichlet边界条件（消去法） fixed_dofs = []; % 存储被约束的自由度编号 for i = 1:size(nodes, 1) if abs(nodes(i,2) - 1) < 1e-12 % 顶盖节点 fixed_dofs = [fixed_dofs, 2*i-1]; % u自由度 elseif abs(nodes(i,1)) < 1e-12 || abs(nodes(i,1)-1) < 1e-12 ... || abs(nodes(i,2)) < 1e-12 % 其余三边 fixed_dofs = [fixed_dofs, 2*i-1, 2*i]; % u和v都约束 end end % 对K和F做边界处理 K(fixed_dofs, :) = 0; K(:, fixed_dofs) = 0; for d = fixed_dofs K(d, d) = 1.0; end F(fixed_dofs) = 0; % 顶盖u=1时改为对应值

这里fixed_dofs存的是全局自由度编号，u和v交替排列，所以节点i的u自由度是2i-1，v是2i。顶盖只约束u，其余边约束u和v。注意K(fixed_dofs,:)=0之后还要把列也置零，否则对称性被破坏，直接求解器可能报奇异。对角置1是保证矩阵非奇异，右端项对应位置设为目标值。

3.3 求解参数与收敛判断

跑方腔驱动流时，雷诺数通过顶盖速度和粘度控制。Re=100时用Picard迭代一般20步内收敛，Re=1000时可能要上百步，而且网格得加密到64×64以上才能捕捉角落涡。收敛判断除了看残差，还要看解的变化量：如果残差降下去了但速度场还在缓慢漂移，说明迭代矩阵的条件数太大，得换预处理或者加网格。

注意：Re=1000以上时，稳态解可能本身就不稳定，迭代不收敛未必是程序问题，可能是物理上进入了非稳态区间。

4. 避坑与排查：有限元求解N-S方程时最容易翻车的五个地方

4.1 压力场棋盘格振荡

现象：压力云图出现高频交替的棋盘格图案，速度场看起来正常但压力完全不可信。原因：速度压力同阶插值，不满足LBB条件，离散后的压力块秩亏。解决：换成Taylor-Hood单元，速度二次压力一次，或者加压力稳定项。这份程序默认用Taylor-Hood，如果自己改了单元类型要重新验证。

4.2 残差下降后突然反弹

现象：迭代前20步残差稳定下降，第21步突然跳高一个量级。原因：Picard迭代在雷诺数较高时对流项线性化误差积累，或者网格某处雅可比行列式接近零。解决：降低松弛因子到0.5，检查最小雅可比行列式是否小于1e-10，必要时对网格做光顺处理。

4.3 边界条件施加后矩阵奇异

现象：Matlab报“Matrix is singular to working precision”。原因：消去法施加边界时只置零了行没置零列，或者固定自由度列表里有重复编号。解决：确保行列同时置零，对角置非零值，用unique去重固定自由度列表。

4.4 收敛但解明显不对

现象：残差降到1e-10，但速度场对称性不对或者量级差几个数量级。原因：单元连接顺序错误导致雅可比为负，或者右端项在边界处理时被错误置零。解决：检查每个单元的节点排列是否逆时针，顶盖边界的目标值是否正确写入右端项。

4.5 网格加密后反而不收敛

现象：N从32加到64，同样的参数迭代发散了。原因：网格加密后条件数变差，直接求解器的舍入误差放大，或者时间步长（伪瞬态）没调整。解决：改用稀疏矩阵存储，加对角预处理，或者先用粗网格收敛解插值到细网格当初值。

5. 进阶用法：用伪瞬态延拓和网格自适应把Re推到更高

5.1 伪瞬态延拓：把稳态问题当时间推进来解

稳态N-S方程在Re较高时，Picard和Newton都容易发散。一个实用技巧是加伪时间项，把稳态问题改写成伪瞬态形式，每步推进一点，相当于用时间步的耗散来稳定非线性迭代。具体做法是在动量方程里加一项(U - U_prev)/dt，dt取0.1到1.0之间，随着残差下降逐步增大dt，最后dt趋于无穷就退化成稳态求解。

% 伪瞬态延拓主循环 dt = 0.1; % 初始伪时间步 dt_max = 1e6; % 最大伪时间步 for iter = 1:200 % 组装含伪时间项的矩阵 [K, F] = assemble_system_pseudo(U_prev, V_prev, mesh, bc, dt); dU = K \ F; U_new = U_prev + dU; res = norm(F - K * U_new) / norm(F); % 残差下降则增大dt，加速收敛 if res < res_prev dt = min(dt * 1.5, dt_max); else dt = max(dt * 0.5, 0.01); % 反弹则减小dt end res_prev = res; if res < 1e-8 break; end U_prev = U_new; end

这段代码的关键在dt的自适应调整：残差下降就放大步长加速，反弹就缩小步长增加耗散。assemble_system_pseudo比标准组装多了一个质量矩阵项除以dt。实际用的时候，Re=1000的方腔驱动流用这个方法通常100步内能收敛，比纯Picard稳得多。

5.2 网格自适应：在哪里加密比加多少更重要

均匀加密网格计算量增长太快，自适应加密只在速度梯度大的地方加节点。常见做法是用后验误差估计，比如基于速度梯度的跳量或者残差范数，标记误差大的单元进行二分加密。这份程序里留了误差估计的接口，可以按单元计算速度梯度模，超过阈值的单元一分为四。

5.3 验证方法：用解析解和网格收敛性检查

跑通算例不等于程序正确。我一般会做两步验证：第一步用有解析解的算例，比如Poiseuille流，检查数值解和解析解的L2误差是否随网格加密按理论阶收敛；第二步做网格收敛性分析，用三套网格算同一个问题，看解的变化是否趋于零。如果误差不降反升，多半是边界处理或者单元组装有bug。

从那以后我每次拿到一个新的有限元求解器，都强制走一遍解析解验证和网格收敛检查，不再只看残差曲线好不好看。希望帮到你。

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