简介二维Stokes问题有限元求解程序包基于Matlab实现面向计算流体力学与数值计算学习者尤其适合想通过具体代码掌握有限元方法从网格生成、形状函数、变分形式到矩阵组装与求解全流程的读者。压缩包共135个文件含133个m脚本与2个txt说明大小约93KB以求解器、边界处理及网格生成等模块化代码为主便于分步阅读与二次修改。已有1557人学习下载。程序覆盖二维定常/非定常Stokes方程的速度-压力求解附Robin与应力边界条件处理包含时间推进、索引组装等典型环节读者可对照源码理解拉普拉斯算子、连续性方程离散实现并借助quiver、contourf等可视化验证流场结果。适合作为课程设计或科研入门参考。1. Stokes 二维有限元程序低速流动仿真的地基先跑通再谈精度微流控芯片里的液滴、生物悬液中的小颗粒、地质蠕变甚至润滑薄膜流体运动的共同点是雷诺数低到惯性项可以忽略剩下粘性力和压力在较真。这里说的 Stokes 二维有限元程序就是用有限元求解这一极限下的不可压缩流动。它看起来比 Navier–Stokes 少了几项实际上反而更拧巴因为速度场和压力场被“不可压缩”这个约束绑在一起单元配对、压力唯一性、边界条件处处有坑。下面的路线用 MATLAB 搭一套可运行的二维 Stokes 有限元程序目标是让你用最小代价拿到可验证的结果。我默认读者已经会写最基本的泊松方程有限元比如知道单元装配和稀疏矩阵但不一定理解“混合元”。如果只想调用别人写好的工具箱帮助有限如果你想把 P2-P1 混合元从头写清楚或者正在做方腔流、低速绕流的基准测试跟着这套思路走至少能少翻几次压力震荡的车。后面四章会依次讲离散逻辑、装配代码、常见问题和收敛性验证。2. 二维 Stokes 方程的离散逻辑为什么速度压力必须混合且单元配对会翻车2.1 Stokes 是 Navier–Stokes 的极端简化但简化后的数学反而更“拧巴”低速流动从 Navier–Stokes 出发ρ(∂u/∂t u·∇u) -∇p μΔu f再加上 ∇·u 0。特征雷诺数很低时左侧的非线性项和时间项远小于粘性项可以剥掉稳态方程退化为-μ Δu ∇p f ∇·u 0。带上合适的边界条件就是二维 Stokes 问题。这里注意两点一是方程从二阶向量偏微分方程变成了“向量弹性方程”的极限二是压力没有出现在任何一个时间演化项里它本质上是一个拉格朗日乘数用来把速度场约束在无散度空间里。这个数学结构带来的直接后果是速度和压力必须作为两个未知量一起解。如果你想把压力消掉比如用流函数方法会得到一个四阶双调和方程虽然可行但边界条件表达很绕也不方便扩展到非牛顿模型。用混合元把压力保留下来才是二维 Stokes 程序最常见的稳定做法。从变分形式看问题变成找速度 u∈V 和压力 p∈L2_0使得a(u,v) b(v,p) F(v) 对所有 v∈V0 b(u,q) 0 对所有 q∈L2其中 a(u,v)μ∫Ω ∇u:∇v dxb(v,p)∫Ω p div v dx。离散后得到一个对称但不定号的线性系统[[K, B^T], [B, 0]] * [[U],[P]] [[F],[0]]这里的 K 是速度刚度块B 是散度约束块右下角是零。这个形式就是典型的鞍点系统不能用求解正定矩阵的直觉去套。很多初学者看到这里的右下角零块第一反应是“系统是不是奇异了”。说实话它确实奇异压力常数是 B 的零空间成员因为 ∫Ω div v dx 对闭域或纯 Dirichlet 边界在强意义下等于零。也就是说如果不加 ∫p0 或固定一个压力参考点线性求解器一定会报奇异。这一点不是 bug是 Stokes 的物理和数学常识必须先接受。2.2 不稳定单元配对为什么翻车P1-P1 棋盘压力与 P2-P1 的选择逻辑有限元理论里有一个不成文的“安全红线”如果速度空间和压力空间选择不匹配离散后的 inf-sup 条件不满足压力场就会高频振荡典型表现是棋盘格。最常见的错误配法就是速度 P1、压力 P1。两者都只有顶点自由度都线性看起来程序好写好调但算出来的压力场在一半单元上高一半低。在结构化网格上运气好可能看不出大问题一旦换三角形网格或者加密压力图立刻变成黑白格子。原因是压力空间“太强”在离散层面有大量的伪模态数值上等于你额外允许了一些不必要的压力变化。稳定配对的原则是速度比压力高或者等效地用稳定化项补足。三角形单元里最常用的稳定组合是 P2-P1速度用六节点二次三角形压力用三节点线性三角形。速度自由度分布在顶点和边中点压力只在顶点定义。这样压力空间被稀释inf-sup 就能在统一网格尺度上成立不会产生棋盘压力。它的实现成本在可接受范围内局部速度基函数数是6压力基函数数是3组装时多处理一组中点编号但边界条件、稀疏矩阵结构都还在可控范围。有人会问为什么不直接把压力取为分片常数 P0P2-P0 也是公认稳定组合之一但压力走得太平参数估计和结果显示都不友好压力在单元之间完全跳跃更容易出现“锯齿”后处理。P2-P1 更接近连续压力兼顾了稳定性和工程可读性所以下面装配阶段就按这个组合来做。在选择 P2-P1 时要意识到另一个隐藏点它只是“足够稳定”不是“无限稳定”。如果网格中单元的几何质量很差比如有接近退化的扁三角形inf-sup 常数会在数值上变小系统矩阵条件数变大压力可能依然出现局部振荡。真正动手前网格检查不能省至少把最小角或者最长边/最短边比打印出来发现畸形单元要么加密要么换网格。3. 在 MATLAB 里从零搭一套 P2-P1 二维 Stokes 程序最少代码跑通3.1 数据结构和自由度编号边中点怎么变成全局节点写有限元程序第一个习惯是“先把自由度编号定死再写装配”。不要一上来就写矩阵循环先想清楚数据结构。节点数组 nodes 按 Nn×2 存每行一个坐标 [x y]。单元数组 elems 按 Ne×3 存每行一个三角形的三个顶点编号。P2-P1 要求速度单元有六个节点其中三个顶点和三条边的中点。中点节点必须是全局唯一的否则同一位置算多个自由度矩阵就没法配对。常见做法是用一个边到新节点的映射表。function [elems6, nodes, edgeMap] expand_p2(nodes, elems) % elems6: Ne x 6, 前3列顶点, 后3列边中点 % nodes: 原节点基础上追加边中点 edgeMap containers.Map(KeyType,char,ValueType,double); elems6 zeros(size(elems,1), 6); nNode size(nodes,1); nextNode nNode 1; for e 1:size(elems,1) i1 elems(e,1); i2 elems(e,2); i3 elems(e,3); edges [i1 i2; i2 i3; i3 i1]; mid zeros(1,3); for k 1:3 key sprintf(%d_%d, min(edges(k,:)), max(edges(k,:))); if edgeMap.isKey(key) mid(k) edgeMap(key); else mid(k) nextNode; edgeMap(key) nextNode; nodes(nextNode,:) 0.5*(nodes(edges(k,1),:) nodes(edges(k,2),:)); nextNode nextNode 1; end end elems6(e,:) [i1 i2 i3 mid(1) mid(2) mid(3)]; end end这段代码先把原有顶点保留再对每条边查映射表如果这个边中点已经被别的单元建过就直接复用否则新建节点。注意节点的追加会让 nodes 数组不断变长示例里用动态追加是为了直观真实程序建议先统计一遍边数再预分配空位。关键点在于边的键用“小编号_大编号”来存这样两个相邻单元访问同一条边时能命中同一个键不会因为方向不同而重复建点。如果忘记统一方向矩阵会自动出现孤立压力自由度结果完全对不上。自由度编号我推荐按“先速度、后压力”排二维速度每个节点占两列自由度分别是 x 方向和 y 方向压力是三节点 P1 的三个自由度排在速度块最后。这样做的好处是施加 Dirichlet 条件时只操作速度自由度压力约束可以单独追加一行不用重构矩阵。3.2 参考单元基函数和装配循环现在写参考三角形上的 P2、P1 基函数。参考单元取 (0,0)、(1,0)、(0,1)三个归一化坐标 L11-xi-etaL2xiL3eta。二次节点按“顶点1-3、中点12/23/31”的顺序排列。function [Nv, dNv] shape_p2_ref(xi, eta) % 返回参考单元上6个P2基函数的函数值和梯度 % Nv: 6x1, dNv: 6x2, 每行是 [dN/dxi, dN/deta] L1 1 - xi - eta; L2 xi; L3 eta; Nv [L1*(2*L1-1); L2*(2*L2-1); L3*(2*L3-1); 4*L1*L2; 4*L2*L3; 4*L3*L1]; dNv zeros(6,2); dNv(1,:) -(4*L1-1)*[1 1]; dNv(2,:) [4*L2-1, 0]; dNv(3,:) [0, 4*L3-1]; dNv(4,:) [4*(L1-L2), -4*L2]; dNv(5,:) [4*L3, 4*L2]; dNv(6,:) [-4*L3, 4*(L1-L3)]; end function [Np, dNp] shape_p1_ref(xi, eta) L1 1 - xi - eta; L2 xi; L3 eta; Np [L1; L2; L3]; dNp [-1 -1; 1 0; 0 1]; endP1 的梯度在参考单元上是常数很好处理。P2 基函数在单元内部光滑变化所以每个高斯积分点都要重新计算一次函数值和梯度不能只在单元中心取一次。单元坐标映射是线性的所以雅可比矩阵 J 是常数可以提前算一次J [x(2)-x(1), y(2)-y(1); x(3)-x(1), y(3)-y(1)]; detJ abs(det(J));物理梯度用grad_phys dN_ref / J完成。注意这里不是普通的矩阵乘而是要理解为“参考梯度对雅可比求逆再变换”。MATLAB 里用矩阵右除直接写dN_ref / J等价于dN_ref * inv(J)对每行做梯度变换。装配时用三点高斯积分参考三角形上权重都是 1/6积分点是 (2/3,1/6,1/6)、(1/6,2/3,1/6)、(1/6,1/6,2/3)。这个套路的精度足够覆盖二次速度、线性压力乘积里出现的四阶多项式是单元刚度矩阵不会欠积分的标准做法。3.3 单元循环速度刚度块和压力散度块一起装接下来是正文中最核心的一段把每个 P2-P1 单元对总刚的贡献拆出来装。function [K, B, F] assemble_stokes_p2p1(nodes, elems6, mu, force) % nodes: 全部节点坐标, N x 2 % elems6: Ne x 6, 六个P2速度节点编号, 前3是顶点 % mu: 动力粘度, 可以是常数 % force: N x 2 的函数句柄, 在每个高斯点取体积力 nn size(nodes,1); K spalloc(2*nn, 2*nn, 12*nn); B spalloc(nn, 2*nn, 3*nn); F zeros(2*nn, 1); gpt [2/3 1/6; 1/6 2/3; 1/6 1/6]; gw [1/6 1/6 1/6]; for e 1:size(elems6,1) vloc elems6(e,:); x nodes(vloc,1); y nodes(vloc,2); J [x(2)-x(1), y(2)-y(1); x(3)-x(1), y(3)-y(1)]; detJ abs(det(J)); invJT inv(J); Kloc zeros(12,12); Bloc zeros(3,12); Floc zeros(12,1); for q 1:3 xi gpt(q,1); eta gpt(q,2); [Nv, dNv] shape_p2_ref(xi, eta); [Np, ~] shape_p1_ref(xi, eta); grad dNv * invJT; % 物理梯度, 6x2 wgh gw(q) * detJ; % 速度刚度块: K mu * grad*grad for i 1:6 for j 1:6 gij grad(i,1)*grad(j,1) grad(i,2)*grad(j,2); Kloc(2*i-1, 2*j-1) Kloc(2*i-1, 2*j-1) mu*wgh*gij; Kloc(2*i, 2*j) Kloc(2*i, 2*j) mu*wgh*gij; end end % 散度约束块: B(p, v) ∫ p div v for p 1:3 for v 1:6 Bloc(p, 2*v-1) Bloc(p, 2*v-1) wgh*Np(p)*grad(v,1); Bloc(p, 2*v) Bloc(p, 2*v) wgh*Np(p)*grad(v,2); end end % 体积力向量 fq force(nodes(vloc,:), xi, eta); % 按需要自行定义 for v 1:6 Floc(2*v-1) Floc(2*v-1) wgh*Nv(v)*fq(1); Floc(2*v) Floc(2*v) wgh*Nv(v)*fq(2); end end % 按自由度编号放入全局矩阵 dof zeros(1,12); for i 1:6 dof(2*i-1) 2*vloc(i)-1; dof(2*i) 2*vloc(i); end K(dof,dof) K(dof,dof) Kloc; pLoc vloc(1:3); B(pLoc, dof) B(pLoc, dof) Bloc; F(dof) F(dof) Floc; end K (KK)/2; end这里有两个必须说清楚的细节。一是 K 块组装时 x 和 y 两个速度分量各自都要乘梯度内积所以 Kloc 的对角小方阵里x 分量和 y 分量是同一套值但不要写错到交叉位置。二是 B 块里的 p 循环是压力节点编号v 循环是速度节点编号压力用的是 P1 基函数速度用的是 P2 基函数别用成长度统一的 6 节点基函数。体积力 force 那段在示例里是简化版本真实使用时我一般直接传(P) [fx(P(:,1),P(:,2)), fy(P(:,1),P(:,2))]的函数句柄再在每个高斯点取坐标求值。如果体积力为零就把 Floc 全部置零不做多余计算。全局矩阵推荐用spalloc先分配不然单元数上去以后zeros会把你内存瞬间打穿。上面的稀疏预分配是按每个速度节点最多和 6 个邻居耦合粗略估的够大多数二维网格用。3.4 边界条件与零均值压力最后一个容易被忽略的环节把 K 和 B 装好后得到鞍点系统nb size(nodes,1); % 这里指压力节点数取装配后的总节点数 S [K, B; B, zeros(nb,nb)]; rhs [F; zeros(nb,1)];但直接用这个矩阵解\会警告矩阵奇异。因为压力自由度里有一个常数自由度是自由漂浮的需要额外约束“平均压力为零”。常见做法是固定某一点的压力比如让 p(1)0这样会把边界附近压力定死对某些非光滑问题会放大该点附近的误差。我更推荐追加一整行拉格朗日乘数C zeros(1,nb); C(:) 1/nb; % 近似 ∫ p da0更精确可以按单元面积加权 S2 [K, B, C; B, zeros(nb), zeros(nb,1); C, zeros(1,nb), 0]; rhs2 [F; zeros(nb,1); 0]; sol S2 \ rhs2;这个 C 行实际是在强制所有压力节点求和平均为零。如果网格均匀、单元面积差异不大效果已经足够如果网格加密后单元面积差异很大建议把 C 改成对每个压力自由度按所属单元面积加权。加权方式是把压力自由度看成一个分片线性函数C 的每个分量就是 ∫Ω Np_i dx 的值这个值在每个单元上等于面积/3累加到三个顶点压力自由度上。Dirichlet 边界条件也要在求解前处理掉。我惯用的做法是先把边界速度自由度的号挑出来然后把对应方程行清零对角线置 1右端项填入已知速度值。注意如果同一行同时被 K 和 B 影响清零行不会破坏散度约束的“压力唯一性”因为 B 的零空间已经靠 C 行压住了。程序里可以这样写dirichletDofs find(bcFlag1d); for k 1:length(dirichletDofs) row dirichletDofs(k); val bcValue(k); S2(row,:) 0; S2(row,row) 1; rhs2(row) val; end这个方式最容易理解缺点是会保留大量已知行矩阵规模略变大。对千量级自由度以内的测试程序性能完全够。到了上万自由度以后更好的做法是把自由分量和约束分量分块只对一个缩小的 free 系统求解。对跑通目标来说不优化反而更容易排查错误。求解之后把 U 和 P 从 sol 里拆出来再按自由度编号散回节点坐标上就能画速度云图和压力云图。4. Stokes 有限元程序常见问题排查奇异矩阵、压力震荡、边界条件三座山4.1 系统矩阵报奇异压力核没有处理还是刚体模态没锁死现象用 MATLAB 的\求解时直接警告 “Matrix is singular to working precision”或者计算出的速度场明显漂移整个场向一个方向平移。原因Stokes 系统的零空间有两个来源。第一压力常数的零向量也就是所有压力自由度同时加一个常数方程仍然成立第二速度空间里如果 Dirichlet 边界加不够刚体平移和旋转会保留下来相当于结构问题里的悬臂梁没有加约束。很多人只注意到了压力常数忘了速度刚体模态。解决先做两步检查。第一步打印 Dirichlet 速度自由度数看它是否锁住每个连通边界上的独立刚体分量第二步给系统加一个零均值压力约束。注意顺序很重要如果边界条件是从外部读进来的最好先对自由边界做一次“连通检查”确认每个边界段都至少有一个速度法向或切向约束否则速度零空间依然残留。常见做法是固定速度场中某个内部节点的两个分量或者把边界闭合让整个区域四壁都是无滑移速度零空间自然消失。还有一个容易被忽略的元凶组装时 B 块散度约束可能少装了一列导致矩阵虽然不奇异但压力方程里有一部分速度自由度没被约束最后结果噪声极大。把 B 的非零模式单独画出来看看每一列是否都有“足够”的非零项这是一个很有效的诊断手段。4.2 压力场棋盘格单元配错还是局部自由度顺序没对上现象压力场不是光滑过渡而是像棋盘一样相邻单元一高一低交替加密网格也不消失反而更明显。原因最常见的是用了不稳定的速度、压力同阶单元比如 P1-P1其次就是明明打算用 P2-P1但装配时把压力自由度错误地接在了六个速度节点的某几个上例如让压力也在边中点定义。表面上看程序不报错实际上压力空间变大稳定性质已经悄悄退化成不稳定模式。解决确认单元配对的三个地方。第一压力局部数组取vloc(1:3)只取顶点不要按速度的六个节点长度去循环第二P1 基函数只有三个高斯点循环里没有Np(4)这种越界第三如果使用了约束或罚函数稳定化检查是否同时对速度剪应力做了调整。手头没有稳定化的 P1-P1 程序不要去调参数碰运气直接换 P2-P1 或 Q2-Q1 重装比在一个不稳定的骨架上打补丁省时得多。另外要注意非结构化网格下棋盘压力并不仅仅表现为“颜色交替”可能是一簇高频振荡模式。最直观的排查是把压力自由度按节点坐标画散点图看零均值附近的波动是否呈高频小尺度分布。如果是基本可以判断是 inf-sup 出了问题不是后处理污染。4.3 矩形角点压力尖峰网格越密压力越大不一定是程序 bug现象方腔驱动流里在靠近移动顶盖和固定侧壁的两个角点压力出现很高的尖峰细化网格后尖峰更高甚至局部等值线扭曲。原因这个角点上速度方向突变从无滑移侧壁变到切向运动的顶盖数学上存在角点奇异性。压力在角点附近的奇异强度取决于夹角和边界条件组合一般情况下是真奇异不是有限元误差。如果不用特殊奇异元网格越密能解析的高频成分越多压力峰值自然升高。解决工程判断上先别急着骂求解器。用一组加密网格观察角点附近速度场的收敛行为如果速度收敛正常压力只在角点一个很小的邻域异常就可以接受。需要更准确结果时可以采用角点奇异增强或局部加密但针对一个 Stokes 程序来说更务实的做法是在后处理里把角点附近压力云图的显示范围截断或者直接读远离角点的压力值。测试基准收敛阶时也要避开角点邻域取样否则误差范数会被局部奇异污染。4.4 迭代求解器残差不降或下降极慢鞍点矩阵的性质导致普通方法失效现象用小规模测试矩阵反斜杠没问题换成上万自由度后改用迭代法比如用普通 CG迭代几百步残差一直震荡甚至NaN。原因混合矩阵是鞍点系统既有正特征值也有负特征值普通要求正定的 CG 根本不该用。其次即使用了 GMRES没有预条件时条件数随网格细化变差速度块和 Schur 补的特征值分布可能非常散残差曲线长期不下降是常态。解决确定迭代方案时明确分层。规模不大直接用\MATLAB 对对称不定系统会走对称不定分解对调试期间够用。规模大了先做一个块对角预条件器主对角取速度刚度矩阵的近似逆压力部分取一个与压力质量矩阵成比例的量再套一个 GMRES 或 MINRES。如果接触到了工程级程序更彻底的常见做法是引入基于压力质量矩阵的 Schur 补近似但不要在第一版 Stokes 程序里做因为参数扰动会掩盖正确性判断。先把直接解的基准结果稳定住再谈迭代优化。5. 进阶验证技巧用网格收敛测试把“算对”变成能查的数字搞定一个能跑的程序只完成了三分之一。我给自己定的习惯是任何新的 Stokes 程序先用一个光滑解析解做收敛测试。制造解取流函数 ψx^2(1-x)^2 y^2(1-y)^2速度定义为 u(∂ψ/∂y, -∂ψ/∂x)发散为零压力随便取一个光滑函数比如 px^4-y^4代入 Stokes 方程反算出体积力 f再用有限元去解同一组 f。按不同网格尺寸跑统计误差。h [0.125, 0.0625, 0.03125]; errU zeros(size(h)); errP zeros(size(h)); for k 1:3 [nodes, elems, bc, f] make_mesh_and_source(h(k)); [u,p] solve_stokes(nodes, elems, bc, f); errU(k) compute_L2_error_u(u, u_exact, nodes, elems); errP(k) compute_L2_error_p(p, p_exact, nodes, elems); if k 1 rateU log(errU(k-1)/errU(k)) / log(2); rateP log(errP(k-1)/errP(k)) / log(2); fprintf(h%.5f rateU%.2f rateP%.2f\n, h(k), rateU, rateP); end end理论上P2-P1 元在网格均匀、解充分光滑时速度 L2 误差的收敛阶接近 3压力 L2 误差的收敛阶接近 2。实际网格里二阶误差已经足够说明程序没有低级错误。如果速度阶只有 1 左右优先检查雅可比变换和边界条件如果速度阶接近 3 但压力阶只有 1多半是压力零均值处理或者 B 块某列装配有错。把误差收敛阶做成表格贴在自己笔记里比单纯看一张云图可靠得多。云图好看不能排除整体偏差收敛阶却能暴露“看似合理但实际不对”的程序。这个习惯帮我挡下过至少三次组装错误每一次都是图上看不出问题、但收敛阶断崖式下跌才暴露的。如果时间允许再按相同思路跑一个方腔驱动流的基准题速度与经典数值结果对照。方腔流不需要解析解但需要对顶盖的切向速度赋值和角点奇异性做注释。有了这两道验证后续换非线性模型、换时间推进会踏实很多。希望帮到你。本文还有配套的精品资源点击获取