简介：一套基于Matlab的二维Stokes问题有限元求解程序，面向学习偏微分方程数值解及流体力学仿真的初学者和科研人员，也适合数值计算课程备课参考。压缩包共135个文件，以133个m脚本和2个txt说明文件为主，整体仅93KB，代码精炼且模块划分清晰，完整覆盖网格生成、变分离散、矩阵组装、线性系统求解与流场可视化等核心环节，并包含瞬态求解及多种边界条件（如Robin、应力边界）的处理实现。相比纯理论推导，这份资料通过可直接运行的代码展示了不可压缩无粘流体稳态与瞬态问题的求解路径，有助于从代码层面理解有限元方法的实际落地。已有1557人学习下载，程序注释清楚、结构便于二次开发，可快速修改算例或拓展到类似偏微分方程问题，是入门到进阶阶段的高性价比参考资料。
1. 为什么是 Stokes 而不是 Navier-Stokes：二维有限元的第一道门槛
做低速流动仿真的人迟早会遇到 Stokes 方程。无论是蠕变流动、多孔介质渗流，还是 MEMS 里的微尺度流道，雷诺数低到惯性项可以忽略时，Navier-Stokes 直接退化成一组线性定常方程。这种场景下压力与速度强烈耦合，二维有限元程序的处理难度集中在混合元的稳定性上，而不是时间推进或湍流模型上。这份 stokes（matlab）二维有限元程序，把网格生成、P2P1 混合元装配、边界置入、求解与后处理串成了一条完整链路，适合正在学有限元方法、或需要快速拿到一个可改动的 Stokes 求解器做对比验证的从业者。我拆过不少这类代码，先提醒一句：别指望它能直接算高雷诺数流动，它的价值在于把混合元的骨架讲清楚，并且每一条矩阵组装逻辑都能在 MATLAB 里逐行调试。
2. 混合元与稳定条件：先搞懂 inf-sup 再写代码
2.1 速度-压力同阶插值为什么翻车
很多初学者拿到 Stokes 方程，第一反应是把速度和压力都用一阶线性三角形元离散，觉得这样自由度少、实现快。真跑起来就发现，压力场出现棋盘状分布，相邻单元压力一高一低交替排列，速度场看着平滑但压力完全不能用。这不是边界条件写错了，而是离散空间没满足 Ladyzhenskaya-Babuska-Brezzi（LBB）条件，也就是常说的 inf-sup 稳定性条件。
从代数角度看，Stokes 离散后得到一个鞍点系统：
% 鞍点系统结构：A 是速度刚度矩阵，B 是散度矩阵 % [A B'] [u] [f] % [B 0 ] [p] = [g]速度空间的维数必须明显大于压力空间的维数，且压力空间不能产生虚假的零能量模式。用同阶插值时，压力自由度相对速度超量，零空间里混入了非零的刚性模态，表现为棋盘压力。判断一个单元组合能不能用，最简单的方式是直接看速度节点数与压力节点数的比例，P1P1 在三角形网格上几乎必然翻车，P2P0 虽然压力分片常数，但总节点数上压力自由度太少，也容易锁死。工程上的常规做法是让压力插值比速度低一阶，同时保证同一条边上的速度节点能约束住压力变化。
2.2 Taylor-Hood 元（P2P1）是默认选择
这份资源里用的正是 P2P1 Taylor-Hood 元：速度单元采用二次三角形元（每个单元 6 个节点：3 个顶点加 3 个边中点），压力采用线性三角形元（每个单元 3 个顶点节点）。这个组合满足 inf-sup 条件，是二维 Stokes 有限元最常见的稳定组合，MATLAB 实现也有成熟的节点编号套路可以抄。
为什么不用更高阶的 P3P2？精度确实更好，但实现复杂度按单元形状函数个数翻倍。做教学验证和中型网格规模计算时，P2P1 的收敛阶已经足够：速度在 H1 范数下二阶收敛，压力在 L2 范数下二阶收敛，实际算顶盖驱动流时，64×64 网格下的涡心位置误差已经能压到百分位以后。对一份要兼顾可读性和可跑性的程序，P2P1 是性价比最高的选择。
2.3 这份资源里的网格与自由度组织
网格生成部分用的是结构化四边形剖分后按对角线拆成三角形的思路，好处是能复用现成的二维网格生成逻辑，坏处是边界附近没有加密能力。自由度编号顺序很关键：程序把速度自由度分成 x 分量和 y 分量两组分别编号，压力自由度单独一组。这样组装出来的矩阵不是按物理位置交织的，而是分块结构，A 矩阵内部还能拆出两个对称的标量刚度块。
% 自由度编号示意（伪代码逻辑） % u_dof = 1:2*Nvel; % 先排速度 x、y 分量 % p_dof = 2*Nvel+1 : 2*Nvel+Npres; % 压力排最后 % 这样 K = [A_xx 0 Bx' % 0 A_yy By' % Bx By 0 ]这种编号方式让矩阵装配循环非常清晰：遍历三角形单元，查单元顶点和中点的全局编号，逐个填充 13×13 的单元矩阵（速度 12 个自由度加压力 3 个自由度）。要注意的是，压力节点和速度顶点节点共用同一套几何顶点，但自由度编号不同。我刚开始拆这类程序时在这里吃过亏，以为编号对上了，结果装配时把压力自由度写进了速度刚度块里，矩阵不对称，求解器直接报错。正确做法是建一个全局节点到速度自由度和压力自由度的映射表，两个表分开建，装配时分别查表。
3. 装配与边界处理：把 Stokes 方程组离散成可解的代数系统
3.1 速度刚度矩阵 K 与散度矩阵 B 的装配
单元装配是整个程序的心脏。P2 单元有三个顶点和三个边中点，形状函数由六个二次基函数构成，压力 P1 单元则用三个线性基函数。对每个三角形单元，需要计算三个基本矩阵：速度刚度贡献、散度耦合项、以及压力-压力块（鞍点系统里是零块）。
% 单个三角形单元的 Stokes 装配（核心循环片段） % nodes_tri: 3x2 顶点坐标; mid_nodes: 3x2 边中点坐标 % 六节点 P2 速度元: 顶点(1,2,3) + 边中点(4,5,6) % 先在参考单元上建立 6 个 P2 基函数和 3 个 P1 基函数 % 再用仿射变换映射回物理单元：x = x1 + J * xi % 实际循环中： for i = 1:6 for j = 1:6 % K_vv(i,j) = ∫ ν * (∇φ_i · ∇φ_j) dΩ K_vv(i,j) = nu * (Dphi_x(:,i)' * wJ * Dphi_x(:,j) + ... Dphi_y(:,i)' * wJ * Dphi_y(:,j)); end for k = 1:3 % B(i,k) = ∫ ψ_k * div(φ_i) dΩ % 注意这里不是对称的，B 由速度基函数的散度投影到压力基函数 B(i,k) = Dphi_x(:,i)' * wJ * psi_k(:,k) + ... Dphi_y(:,i)' * wJ * psi_k(:,k); end end这段代码的逻辑分三层：第一层用参考单元上的高斯积分点计算 P2 基函数对 x、y 的偏导数，乘上雅可比行列式得到物理单元上的积分权重；第二层把基函数梯度两两相乘累加得到速度刚度块，粘性系数 nu 作为常数提到外面；第三层把速度基函数的散度和压力基函数相乘，得到散度耦合矩阵。参数上最需要注意的是 nu 的取值，程序里默认 nu=1，对应雷诺数为零的纯 Stokes 流动。如果你的问题带有特征速度 U 和特征长度 L，nu 应该取为 1/Re 而不是 1，否则解得的速度量级会对不上物理模型。另一个容易错的地方是雅可比行列式在每个高斯点上是不变的常量（仿射变换），但代码里 wJ 如果写成逐点变化的向量，积分权重计算就会重复累加，导致刚度矩阵整体偏大，压力场出现奇怪的放大。
3.2 压力自由度与零均值约束
Stokes 方程组的压力只出现在梯度项里，这意味着压力加减一个常数仍然是解。纯 Neumann 边界条件或者全 Dirichlet 速度边界条件下，压力解不唯一，代数系统里对应一个零特征值。不处理这个奇异模态，直接求解会得到 NaN 或者残差卡住不降。标准做法是给压力加零均值约束，程序里用的是拉格朗日乘子法，也就是在系统矩阵的右下角加一行一列。
% 零均值压力约束：∫ p dΩ = 0，用拉格朗日乘子 c 实现 % 扩展一行: [K B' 0; B 0 e; 0 e' 0] * [u; p; c] = [f; 0; 0] % 其中 e 是全 1 向量，作用在压力自由度上 % MATLAB 实现时： e = ones(Npres,1); K_total = [K_spd, B_spd'; B_spd, sparse(Npres,Npres)]; K_total = [K_total, sparse(size(K_total,1),1); sparse(1,size(K_total,2)), 0]; K_total(end-Npres:end-1, end) = e; K_total(end, end-Npres:end-1) = e'; f_ext = [f_vel; zeros(Npres,1); 0];这种处理方式比罚函数法更干净。罚函数法是把压力零均值约束用一个很大的罚数强行压进对角元，好处是矩阵规模不变，坏处是要调罚数大小，罚得不够约束不生效，罚得太大条件数恶化，直接法求解时舍入误差反而把压力精度毁掉。拉格朗日乘子法多一个自由度，但约束是精确满足的，我一般更推荐这个方案。注意 e 向量只作用在压力自由度块上，不能碰到速度块，否则会把零均值约束错加到速度上，导致速度场整体漂移，看起来像没加边界条件。
3.3 Dirichlet 边界条件的置入方式
速度边界条件有两类典型处理：完全置零法和置大数法。程序里默认用前者，对所有边界节点，先找到对应的速度自由度编号，再把系统矩阵相应的行和列做如下操作：对角线置 1，非对角线清零，右端项改成边界值，同时把该自由度的影响从其他方程的右端项中减去。
% 设第 i 个自由度的边界值为 bc_value % 做法：先取出原列对右端项的贡献再清零 f_rhs = f_rhs - K_total(:, i) * bc_value; K_total(i, :) = 0; K_total(:, i) = 0; K_total(i, i) = 1; f_rhs(i) = bc_value;这段代码的顺序有一点讲究：必须先做 f_rhs 的修正，再清空矩阵行列。因为 f_rhs 的修正利用的是未修改前的 K_total 第 i 列原始数据，如果把矩阵先清空了，这一列的贡献就丢了，边界节点上的值会传导不到相邻节点，结果是速度场在边界附近出现跳变。置大数法虽然实现只有两行（把对角线乘一个大数、右端项乘同样的大数），但在大规模稀疏矩阵里容易让预处理器的对角缩放失效，我见过 ILU 分解在这种矩阵上直接不收敛的情况。程序采用完全置零法，配合 MATLAB 的稀疏矩阵索引能保持矩阵对称性，这是我可以直接照搬的习惯。压力边界不需要额外处理，零均值约束已经解决了压力的非唯一性。
4. 求解器选择与参数调试：别一上来就偷懒用反斜杠
4.1 直接法 vs 迭代法：算到多大规模要换
MATLAB 里求解鞍点系统最省事的方式是直接用反斜杠运算符，程序默认在小网格下也这么干。但你要清楚反斜杠对稠密矩阵和稀疏矩阵的内部处理路径不一样：当系统矩阵是稀疏的且规模小于某个阈值时，它会走稀疏直接分解；一旦网格加密到 128×128 单元以上，自由度总数轻松超过三万，直接分解的速度会明显变慢，内存占用按 O(n^1.5) 增长。
% 直接法（小网格，自由度数 < 2万） u_sol = K_total \ f_rhs; % 迭代法（大网格，自由度数 > 3万） % 注意：必须用对称正定预处理器，否则鞍点系统收敛极慢 [u_sol, flag, relres, iter] = minres(K_total, f_rhs, tol, maxit, L, L');判断要不要换迭代法的临界点通常不在自由度数值本身，而在矩阵的非零元数量。P2P1 元每个节点关联的单元数多，非零元密度比线性元高很多，同样是两万自由度，P2P1 的矩阵分解时间可能是 P1P1 的三到五倍。我在实际中习惯的做法是：先用一种网格规模跑一次分解计时，如果分解时间超过求解时间的几十倍，或者内存接近物理内存一半，就立刻切迭代法直接法切换。细腻的做法是把直接法当作基准校验工具、迭代法当作生产工具，同一套网格两边都算一次，确认迭代解与直接解的最大相对误差在 1e-8 量级内。
4.2 常用迭代器与预处理器的参数设置
鞍点系统矩阵不定，不能用共轭梯度法直接处理，程序里实际可选的迭代器有两个方向：对非对称化后的系统用 GMRES，对原对称鞍点系统用 MINRES。两者在 MATLAB 里都有现成函数，但参数设置上有几个要点经常被忽略。
% MINRES 配合块对角预处理器：分别对 A 块和 Schur 补做 ILU % 这里 L 是预处理器因子，shape 需与 K_total 一致 tol = 1e-8; % 相对残差容差 maxit = 500; % 最大迭代步 % 实际项目中我常用的配置： % Schur 补 S = B * A^{-1} * B'，尽管 A^{-1} 没有显式求 % 常见做法是用压住对角的质量矩阵代替 S 的逆 % 对 P2P1，压力质量矩阵 Mp 的逆预处理效果可接受 % 若直接调用内置预处理器（MATLAB 部分版本支持） options = struct('type','ilutp','droptol',1e-3,'udiag',1); [L, U] = ilu(K_total(1:end-1,1:end-1), options);这里的逻辑分两步：第一步是选预处理子结构，对鞍点系统最有效的其实是块预处理器，它把 A 块和 Schur 补分开处理，但实现复杂。绝大多数现成代码会用单一 ILU 预处理整个扩展矩阵，参数上 droptol 是关键。我实测经验是 droptol 取 1e-3 到 1e-4 之间迭代步数差异不大，但取 0 保全部填充会导致预处理器构建时间和内存暴涨。如果矩阵对角元尺度差异过大（单元形状畸变容易导致），打开 udiag=1 作对角补偿，能避免 ILU 分解遇到零对角元直接报错。第二步是选迭代器，MINRES 比 GMRES 的内存占用稳定，因为不需要存储整个 Krylov 子空间，但 MINRES 要求矩阵对称，加入拉格朗日乘子扩展后矩阵仍然对称，可以放心用。
4.3 收敛判据与残差监控
程序里如果不主动监控残差，迭代法可能在一个看似合理的容差处停住，实际上压力场还没收敛。要注意 MATLAB 内置迭代函数的 tol 是相对残差还是绝对残差，不同版本默认行为有差异。我的通用做法是在迭代结束后手动计算一次真实残差，不信任求解器返回的 flag 为零。
% 求解后手动验证残差（强制步骤） res_u = norm(f_rhs - K_total * u_sol) / (norm(f_rhs) + 1e-12); fprintf('真实相对残差: %.2e\n', res_u); if res_u > 1e-6 warning('残差未达标，检查边界置入或预处理参数'); end这里有几个参数值得说明：分母里加的 1e-12 是防止右端项恰好为零时除零报错；阈值 1e-6 是工程可接受线，做基准验证时要收紧到 1e-8。另一个监控点是压力均值，迭代解由于舍入误差累积，压力均值可能偏离零到 1e-10 量级，这个量级无害，但如果偏离到 1e-3 以上，说明零均值约束的拉格朗日乘子那行在预处理过程中被压缩掉了，需要检查是否对扩展矩阵做了行缩放。程序里残差监控代码放置的位置也讲究：要在后处理脚本之前执行，这样一旦残差异常可以直接断点调试，而不是拿着一个错误的速度场去做流线图，白费几十分钟。
5. 避坑：二维 Stokes 有限元最常见的五个翻车现场
5.1 顶盖驱动流顶盖角点压力发散
现象：顶盖驱动流算出来，右上角和左上角附近压力出现大的正负震荡，速度场看起来正常，但压力等值线图上一片斑驳。
原因：顶盖速度边界条件在角点上同时有 x 方向和 y 方向的约束，而实际物理中角点处速度不连续（顶盖动、侧壁静），直接置入两个方向的 Dirichlet 条件会让角点速度产生奇异应力，压力解在角点处本来就弱奇异。网格不加密时这种奇异被线性插值放大。
解决：最简做法是角点速度不设顶盖值，让顶盖速度从角点向内缩一个节点，或者用渐变过渡：角点赋零速度，紧邻节点的速度按线性插值。程序里如果直接在全局边界数组里把角点写成 (1,0)，就会踩中这个坑。我在复现经典结果时习惯把顶盖速度写成一个函数，角点处返回 0，非角点处返回 1，这样既能保持边界处理逻辑统一，又不会破坏角点约束的相容性。
5.2 压力场整体偏移导致后处理积分异常
现象：算一个管道流动，速度剖面看起来对，但用力积分算压差时结果与理论值差一个常数倍，压力云图色标范围明显不对称。
原因：压力零均值约束只锁定了压力平均值，但某些后处理脚本直接取压力数组前几个值当作参考压力，忽略了压力常数自由度的存在。当边界条件全是 Neumann 或周期型时，压力基准只能靠零均值约束确定，任何单点取值都没有物理意义。
解决：后处理里的压力差计算必须用两个空间点的压力值相减，或者对压力做空间积分后再求差。我一般的做法是抛弃程序里给出的绝对压力云图，强制显示 p - mean(p)，这样色标中心被拉到零，能直接看到压力梯度的分布。另一个相关细节是零均值约束实现时如果 e 向量长度与压力自由度数量不一致，约束会静默失效，跑完后手动检查 mean(p_sol) 是否在 1e-12 量级，这个检查代码值得长期保留在脚本末尾。
5.3 单元编号顺序错乱导致散度矩阵转置错误
现象：速度场与压力场耦合不起来，压力几乎全为零，但速度场又满足边界条件，残差很小甚至为零。这个现象比明显报错更危险，容易让人误判程序正确。
原因：装配散度矩阵 B 时，循环里对单元局部自由度顺序的约定不一致。P2 单元如果先排三个顶点再排三个边中点，那么局部自由度 1~3 对应顶点，4~6 对应边中点；但如果在查全局编号时把边中点的编号映射到了顶点编号上，B 矩阵的列向量就会错位，散度方程约束不到正确的速度模式上，压力失去来源项，自然解成零。
解决：写一个单元连通性自检函数，对每个单元打印局部节点坐标与全局节点坐标，逐一核对边中点是否落在两个顶点坐标的平均位置。这种自检代码虽然简单，但我见过太多人在此翻车，因为网格生成脚本返回的节点数组顺序不是按先顶点后边中的约定来的，而是先全部顶点再全部边中点，装配时必须重新索引。程序里如果发现压力恒为零的异常，先查映射表，不要急着调求解器。
5.4 粘性系数 nu 取值导致速度量级失真
现象：速度场形状正确，但幅值比物理预期小几个数量级，或者无量纲化后对不上基准数据。
原因：Stokes 方程中 nu 不是随便填的。如果特征尺度是米、特征速度是米每秒，nu 取水的运动粘性 1e-6，那么速度解出来会接近零，因为力的量级和粘性量级没匹配。很多程序示例为了教学把 nu 默认设为 1，换到实际工况时忘记改。
解决：先对问题做无量纲化，确定特征长度 L 和特征速度 U，算有效雷诺数 Re = U*L/nu。Stokes 极限下 Re 趋于零，此时 nu 在实际代码中应该被替换为 1（如果力和速度都已无量纲化）或者是具体无量纲参数。我通常在程序配置文件里写清楚 nu 的定义域，并在求解后输出速度最大值做一次量级校验，速度最大值偏离边界驱动值不能超过一个量级，否则回去查 nu。
5.5 稀疏矩阵索引分配失败导致 OOM
现象：网格规模从 64×64 加密到 128×128 时，MATLAB 内存直接耗尽或提示 sparse 索引分配过多非零元。
原因：装配循环中如果使用动态增长的稀疏矩阵赋值方式，比如 K_total(i,j) = K_total(i,j) + val 逐点填充，每写一个元素 MATLAB 都要重新分配和排序索引，非零元数量上万后效率崩塌。更隐蔽的是，同一个位置被多个单元重复累加时，稀疏矩阵的累积计数没有预分配，最终索引数组长度失控。
解决：用稀疏矩阵的三元组构建方式，先开一个三元组数组记录行号、列号、值，装配全部完成后一次性用 sparse(I,J,V,m,n) 构造矩阵。同时统计每个单元的局部自由度数，预先估计非零元总数，把三元组数组预分配为估计值的两倍容量。程序里如果卡在装配阶段，把 sparse 调用改成预先累积 I、J、V 三个向量再做一次性构造，内存占用会下降一个量级以上。
6. 验证从 Lid-Driven Cavity 开始：一个能复现已发表基准的流程
6.1 基准设置
拿到程序后别急着算自己的工况，先在顶盖驱动流基准上跑通。计算域取单位正方形 [0,1]×[0,1]，顶盖速度设为水平方向 u=1、垂直方向 v=0，其余三壁速度全为零，压力零均值约束保留，nu 取 1 对应无穷大雷诺数下的 Stokes 极限。这个基准的好处是解析解不存在，但大量已发表的数值结果给了涡心位置、速度剖面等具体数字，可以直接对照误差。网格建议从 32×32 单元开始，逐步加密到 128×128，记录每个网格下的涡心坐标和中心线速度剖面。
6.2 看哪些量：涡心位置、速度剖面、压力水平
顶盖驱动流的流场呈一个大涡，涡心位置是检验程序可靠性的敏感指标。计算涡心时不能用单元中心的速度插值，那样误差太大，要基于节点速度用线性插值找出两个速度分量同时为零的点。常见做法是先在粗网格上定位大致区域，再在加密网格上精确化。我习惯输出三个验证量：主涡涡心坐标、通过几何中心竖直线上的 u 速度剖面、中心水平线上的 v 速度剖面，误差分别在 1e-3 量级才算通过。
% 验证脚本最后一段：检查涡心位置 % 节点速度场 u_v, v_v 已从解向量中恢复 % 用插值网格逼近 u=0、v=0 的点 [Uc, Vc] = meshgrid(linspace(0,1,200), linspace(0,1,200)); U_interp = griddata(node_x, node_y, u_v, Uc, Vc); V_interp = griddata(node_x, node_y, v_v, Uc, Vc); % 找最小速度模量位置 speed_mag = sqrt(U_interp.^2 + V_interp.^2); [~, idx] = min(speed_mag(:)); vortex_x = Uc(idx); vortex_y = Vc(idx); fprintf('涡心估计: (%.4f, %.4f)\n', vortex_x, vortex_y);这段脚本有一个隐含的逻辑陷阱：速度模量的最小点不一定就是两个分量同时为零的点，在涡心附近速度模量很小但两个分量可能都不到零。严格做法是把 U_interp 和 V_interp 各自穿零点的等值线求交，但工程上只要网格足够密，速度模量最小点与真实涡心的偏差已经小于网格尺寸，可以直接接受。我从某次项目经历里学到的教训是，验证不能只看一个网格，至少要三个逐渐加密的网格，观察涡心坐标是否随网格加密单调逼近稳定值。如果 64×64 和 128×128 的结果跳跃超过 0.01，说明装配有 bug 或者边界条件置入不一致。
最后说一个我自己的强迫症习惯：每次改完程序参数，都强制走一遍同一套基准脚本，并记录涡心坐标提交到备注里。从那以后，我只要看到涡心坐标偏离基准超过万分之五，就知道一定有人动了网格生成或边界置入代码，排查范围一下子缩小到两块。这套流程坚持下来，省掉了我大量把新代码从头翻到尾找 bug 的时间。希望帮到你。
本文还有配套的精品资源，点击获取
