简介这份资源是面向流体力学、数值计算方向的学习者与工程师的二维稳态Navier-Stokes方程有限元求解程序基于Matlab实现适合已具备偏微分方程与有限元基础、希望深入理解FEM求解流程的中高级读者。程序围绕动量方程与连续性方程展开涵盖几何离散、函数空间定义、变分形式推导、矩阵组装、线性系统求解及速度压力后处理等完整环节并涉及Dirichlet、Robin、应力边界条件处理与Newton迭代初始化等细节。压缩包共45个文件全部为m脚本整体约19KB按功能可大致分为网格与有限元空间生成、局部与参考基函数、单元与边界积分、矩阵与向量组装、边界条件处理、解析解与误差计算等模块结构清晰便于按需查阅。目前已有312人学习下载。通过阅读与调试源码读者可掌握有限元法在流体问题中的落地方式并可将框架迁移至热传导、扩散等方程提升Matlab编程与数值算法实现能力。1. 二维稳态Navier-Stokes方程有限元求解从“算不动”到“算得准”的分水岭做 CFD 的人迟早会撞上同一个问题明明方程写对了网格也画了算出来的方腔驱动流却要么残差卡在 1e-3 不动要么速度场出现棋盘格振荡。二维稳态 Navier-Stokes 方程的有限元求解程序解决的正是这类“方程没错、结果不对”的落地困境。它面向的是需要在本地把不可压缩黏性流动算准的工程师——比如做微流道散热、翼型低速绕流、搅拌槽内流场分析的人。和瞬态求解不同稳态问题省掉了时间步进但压力-速度耦合带来的鞍点结构反而更棘手速度的试探空间和压力的试探空间必须满足 LBB 条件否则压力会像脱缰野马一样出现伪振荡。常见做法是采用 Taylor-Hood 单元速度二次、压力一次配合 Newton 迭代处理非线性对流项。这一章先把“为什么稳态 NS 有限元值得单独做一套程序”讲清楚后面几章再拆解弱形式、单元选型、Newton 收敛和避坑细节。2. 弱形式推导与 Taylor-Hood 单元为什么压力不能和速度同阶2.1 从强形式到弱形式分部积分到底改变了什么二维稳态不可压缩 Navier-Stokes 方程的强形式写成-ν ∇²u (u·∇)u ∇p f in Ω ∇·u 0 in Ω u g on Γ_D ν ∂u/∂n - p n h on Γ_N其中 u (u, v) 是速度p 是压力ν 是运动黏度。直接对二阶导数做离散有限元的 C¹ 连续性要求会让基函数构造变得非常麻烦。弱形式的核心动作是对黏性项做分部积分把 ∇²u 的导数负担转移一个到试探函数上∫Ω ν ∇u : ∇v dΩ ∫Ω (u·∇)u · v dΩ - ∫Ω p (∇·v) dΩ ∫Ω q (∇·u) dΩ ∫Ω f·v dΩ ∫Γ_N h·v dΓ这里 v 和 q 分别是速度和压力的试探函数。分部积分之后速度只需要 H¹ 连续性线性单元就能用。但代价是压力 p 不再有导数它变成了一个 Lagrange 乘子约束着 ∇·u 0。这就是鞍点问题的来源离散后的矩阵不是正定的而是对称不定的。注意弱形式里压力项 -∫ p (∇·v) dΩ 和连续性约束 ∫ q (∇·u) dΩ 必须同时存在少写任何一个都会导致压力解完全错误。2.2 LBB 条件与 Taylor-Hood 单元的选型理由离散之后速度空间 V_h 和压力空间 Q_h 不能随便选。它们必须满足 Ladyzhenskaya-Babuska-BrezziLBB条件也叫 inf-sup 条件inf_{q_h ∈ Q_h} sup_{v_h ∈ V_h} |∫ q_h (∇·v_h) dΩ| / (||v_h||_1 ||q_h||_0) ≥ β 0β 与网格尺寸无关。如果违反这个条件压力会出现棋盘格振荡而且加密网格也不会改善。常见的单元组合对比如下单元类型速度阶次压力阶次是否满足 LBB适用场景P1-P1线性线性否需额外稳定化不推荐直接用P2-P1 (Taylor-Hood)二次线性是稳态 NS 最常用稳健P2-P0二次常数否压力不连续需稳定化MINI 单元线性气泡线性是实现简单精度略低P3-P2三次二次是高精度计算量大我一般首选 Taylor-HoodP2-P1。原因很直接它满足 LBB 条件不需要额外的人工稳定项代码实现也不算复杂——速度用 6 节点三角形压力用 3 节点三角形。代价是每个单元的自由度从 3 个速度节点变成 6 个全局矩阵规模大约翻倍但换来的是压力场干净、Newton 迭代收敛稳定。2.3 单元刚度矩阵的组装从局部到全局Taylor-Hood 单元的速度形函数是二次的有 6 个节点3 个顶点 3 个边中点压力形函数是线性的有 3 个顶点。每个单元的局部自由度排列是 [u1, v1, u2, v2, ..., u6, v6, p1, p2, p3]共 15 个。组装的核心是计算以下几类积分import numpy as np from scipy.sparse import coo_matrix def assemble_stokes(nodes, elements, nu): 组装 Stokes 部分的全局矩阵忽略对流项 nodes: (N, 2) 节点坐标 elements: (M, 6) 三角形单元的速度节点索引二次 nu: 运动黏度 n_vel len(nodes) # 速度节点数 n_prs n_vel # 压力节点数线性节点是二次节点的子集 n_dof 2 * n_vel n_prs # 总自由度 rows, cols, vals [], [], [] for elem in elements: # 提取单元节点坐标 xy nodes[elem] # (6, 2) # 计算单元面积和形函数导数二次三角形 # 这里用 3 点高斯积分 area, dNdx, dNdy tri6_shape_deriv(xy) # 黏性项nu * ∫ ∇u : ∇v dΩ # 速度-速度耦合块 (12x12) K_vel np.zeros((12, 12)) for i in range(6): for j in range(6): K_vel[2*i, 2*j] nu * area * (dNdx[i]*dNdx[j] dNdy[i]*dNdy[j]) K_vel[2*i1, 2*j1] nu * area * (dNdx[i]*dNdx[j] dNdy[i]*dNdy[j]) # 压力-速度耦合块-∫ p (∇·v) dΩ # 速度节点 i (二次) 与压力节点 j (线性) 的耦合 B np.zeros((12, 3)) for i in range(6): for j in range(3): # 线性压力形函数在二次节点上的值 Np_j linear_shape_at_tri6(xy, j, i) B[2*i, j] -area * Np_j * dNdx[i] B[2*i1, j] -area * Np_j * dNdy[i] # 组装到全局坐标此处省略索引映射细节 # ... return coo_matrix((vals, (rows, cols)), shape(n_dof, n_dof)).tocsr()这段代码的关键点有三个。第一黏性项的积分用 3 点高斯积分对二次三角形已经足够精确因为 ∇u : ∇v 是二次多项式。第二压力-速度耦合块 B 的维度是 12×3因为速度有 6 个节点各 2 个分量压力有 3 个节点。第三B 的转置 B^T 对应连续性约束 ∫ q (∇·u) dΩ组装时直接复用 B 的转置即可不要重复计算。参数说明nu 是运动黏度对于空气约 1.5e-5 m²/s水约 1e-6 m²/s。area 是单元面积由节点坐标叉积得到。dNdx 和 dNdy 是二次形函数对 x 和 y 的偏导在等参变换下通过 Jacobian 矩阵的逆计算。3. Newton 迭代处理对流项从 Stokes 到完整 NS 的最后一公里3.1 对流项的线性化为什么 Picard 不够快Stokes 问题是线性的组装完矩阵直接求解就行。但完整的 NS 方程多了对流项 (u·∇)u这是非线性的。处理非线性有两种主流做法Picard 迭代也叫定点迭代和 Newton 迭代。Picard 迭代把对流项写成 (u_old·∇)u_new每次迭代解一个线性问题。它的收敛速度是线性的对于雷诺数稍大的问题可能需要几十次甚至上百次迭代而且有时候根本不收敛。Newton 迭代则对残差做一阶 Taylor 展开收敛速度是二次的——残差从 1e-2 降到 1e-10 可能只需要 3 到 4 次迭代。代价是 Newton 需要计算 Jacobian 矩阵也就是对流项对速度的导数。对于二维问题每个速度节点的 Jacobian 贡献是一个 2×2 块∂[(u·∇)u]/∂u [u_x ∂u/∂x, ∂u/∂y] [∂v/∂x, u_x ∂v/∂y]其中 u_x 是当前迭代步的 x 方向速度。这个 Jacobian 的组装比 Picard 复杂但收敛速度的提升通常值得。3.2 Newton 迭代的完整实现与收敛判据def newton_solve(nodes, elements, nu, f_body, bc, max_iter20, tol1e-8): Newton 迭代求解稳态 NS 方程 bc: 边界条件字典包含 Dirichlet 节点和值 # 初始猜测用 Stokes 解作为起点 U solve_stokes(nodes, elements, nu, f_body, bc) for k in range(max_iter): # 组装残差 R(U) 和 Jacobian J(U) R, J assemble_ns_residual_jacobian(nodes, elements, nu, U, f_body) # 施加 Dirichlet 边界条件置大数法或消行消列法 R, J apply_dirichlet(R, J, bc) # 求解 J * dU -R dU spsolve(J, -R) # 更新解 U U dU # 收敛判据残差范数和更新量范数同时检查 res_norm np.linalg.norm(R) du_norm np.linalg.norm(dU) print(fIter {k}: |R| {res_norm:.3e}, |dU| {du_norm:.3e}) if res_norm tol and du_norm tol: print(Newton 收敛) return U print(警告Newton 未在最大迭代次数内收敛) return U逻辑说明初始猜测用 Stokes 解因为 Stokes 问题是线性的求解代价低而且它的解已经满足了连续性方程和黏性项离 NS 的解不会太远。每次 Newton 迭代组装残差 R 和 Jacobian J然后解线性系统 J dU -R。收敛判据同时看残差范数和更新量范数只检查其中一个可能会误判——残差小但更新量大说明还在震荡更新量小但残差大说明可能卡在某个非物理解附近。参数说明max_iter 一般设 20 足够如果 20 次还不收敛要么是雷诺数太高超出了稳态解的存在范围要么是网格太粗导致离散误差主导。tol 设 1e-8 是残差的绝对范数对于无量纲化后速度量级为 1 的问题这个阈值对应大约 8 位有效数字。如果问题本身量级很大比如速度 100 m/stol 要相应放大到 1e-6 左右。3.3 雷诺数升高时 Newton 为什么不收敛这是稳态 NS 有限元最常翻车的地方。雷诺数超过某个临界值后稳态解可能根本不存在流动本质是瞬态的或者存在但 Newton 的收敛域很窄。我踩过的坑是Re100 的方腔驱动流用 Stokes 解做初始猜测Newton 直接发散换成 Picard 先迭代 10 次再切 Newton就能稳定收敛。常见做法是混合策略前 5 到 10 次用 Picard 迭代把解拉到 Newton 的收敛域内再切换到 Newton 加速收敛。代码上只需要在循环里加一个判断if k 10: # Picard 步忽略对流项对速度的导数 J assemble_picard_jacobian(nodes, elements, nu, U) else: # Newton 步完整 Jacobian J assemble_full_jacobian(nodes, elements, nu, U)另一个技巧是雷诺数延拓先算 Re10 的解把它作为 Re50 的初始猜测再作为 Re100 的初始猜测。每次增加的雷诺数不要超过 2 倍否则 Newton 照样发散。4. 避坑与排查稳态 NS 有限元最常见的 5 个翻车现场4.1 压力场出现棋盘格振荡现象速度场看起来合理但压力场在相邻节点之间正负交替像国际象棋棋盘。原因速度空间和压力空间不满足 LBB 条件。最常见的是用了 P1-P1 单元速度和压力都是线性或者虽然用了 Taylor-Hood 但压力节点的自由度编号和速度节点没有正确对应。解决确认速度用二次单元、压力用线性单元。检查组装时压力形函数在二次节点上的取值是否正确——线性压力形函数在边中点上的值是 0.5在顶点上是 1 或 0。如果压力振荡仍然存在检查边界条件是否给了压力一个参考值纯 Dirichlet 速度边界下压力只能确定到相差一个常数需要固定一个节点的压力。4.2 Newton 残差卡在 1e-3 不再下降现象Newton 迭代前几次残差下降很快到 1e-3 左右就几乎不动了继续迭代更新量也很小。原因通常是网格分辨率不够离散误差主导了残差。或者收敛判据的 tol 设得太小低于离散误差的底线。解决先检查网格是否足够细。对于方腔驱动流 Re100至少需要 64×64 的网格才能让残差降到 1e-8。如果网格已经够细把 tol 放宽到 1e-6 试试。另外检查残差的范数类型——用 L² 范数和用 L∞ 范数得到的数值可能差一个量级统一用 L² 范数。4.3 出口边界条件导致回流发散现象入口给速度、出口给压力或自由流出计算在出口附近出现回流Newton 发散。原因出口边界上如果只给“自然边界条件”即 ν ∂u/∂n - p n 0当流动出现回流时这个条件在数学上是不稳定的因为回流会把出口处的扰动带入计算域。解决把出口边界延长让回流区远离出口。或者在出口施加一个弱的 Dirichlet 条件比如只固定压力的参考值速度用“do-nothing”条件。如果回流严重考虑改用瞬态求解——稳态解可能根本不存在。4.4 黏性项积分精度不足导致解不收敛现象加密网格后残差反而上升或者 Newton 迭代出现震荡。原因黏性项的积分用了太低阶的高斯积分。对于二次三角形单元∇u : ∇v 是二次多项式至少需要 3 点高斯积分才能精确积分。如果用了 1 点积分形心积分误差会随网格加密而累积。解决确认高斯积分点数。二次三角形用 3 点或 7 点积分线性三角形用 1 点或 3 点。检查代码中积分点的权重和坐标是否与单元类型匹配。4.5 全局矩阵奇异或接近奇异现象spsolve 报奇异矩阵警告或者解出现 NaN。原因压力场的零空间没有被消除。纯速度边界条件下压力可以加上任意常数而不改变方程导致全局矩阵有一个零特征值。解决固定一个压力节点的值比如令 p[0] 0或者用 Lagrange 乘子法施加压力均值零约束。最简单的方法是在组装完成后把第一个压力自由度对应的行和列消去右端项设为 0。5. 验证与进阶用方腔驱动流标定你的求解器5.1 方腔驱动流的基准数据与验证方法方腔驱动流是稳态 NS 有限元最经典的验证案例单位正方形区域顶边速度 u1其他三边无滑移Re100 到 1000。Ghia 等人在 1982 年给出了高精度基准解常用的是沿垂直中线和水平中线的速度剖面。验证步骤在 64×64 的 Taylor-Hood 网格上算 Re100提取 x0.5 处的 u 速度沿 y 方向的分布和 Ghia 数据对比。如果最大偏差在 2% 以内求解器基本可信。Re1000 时需要更细的网格128×128 以上而且 Newton 可能需要延拓策略。# 提取 x0.5 处的 u 速度剖面 def extract_centerline_u(U, nodes, n_vel): u U[:n_vel] x_target 0.5 tol 1e-6 idx np.where(np.abs(nodes[:, 0] - x_target) tol)[0] y_vals nodes[idx, 1] u_vals u[idx] sort_idx np.argsort(y_vals) return y_vals[sort_idx], u_vals[sort_idx]这个函数返回 x0.5 处的 y 坐标和对应的 u 速度。和 Ghia 数据对比时注意 Ghia 的数据是有限差分算的网格分辨率和你的有限元网格不同插值误差是主要误差来源。5.2 从稳态到瞬态什么时候该放弃稳态求解稳态 NS 方程的解在雷诺数超过临界值后可能不存在。方腔驱动流的临界雷诺数大约在 8000 左右超过这个值流动本质是瞬态的。但实际计算中Re 超过 2000 后 Newton 就很难收敛了即使稳态解理论上存在。我的习惯是Re 1000 用稳态求解Re 1000 直接上瞬态。瞬态求解虽然计算量大但不需要处理 Newton 的收敛域问题时间步进本身就是一个天然的延拓。如果一定要算高雷诺数的稳态解用伪时间步进法在稳态方程里加一个虚拟时间导数项步进到残差足够小再关掉时间项用 Newton 收尾。5.3 一个容易被忽略的技巧用正规方程检查 Jacobian 的正确性Newton 迭代不收敛时第一件事是检查 Jacobian 组装是否正确。有限差分验证是最可靠的方法对每个自由度施加一个小扰动 ε计算残差的变化和 Jacobian 的对应列比较。def check_jacobian(nodes, elements, nu, U, f_body, eps1e-7): R0, J assemble_ns_residual_jacobian(nodes, elements, nu, U, f_body) n_dof len(U) max_err 0.0 for i in range(min(n_dof, 50)): # 只检查前 50 个自由度 U_pert U.copy() U_pert[i] eps R_pert, _ assemble_ns_residual_jacobian(nodes, elements, nu, U_pert, f_body) J_fd (R_pert - R0) / eps err np.linalg.norm(J_fd - J[:, i]) / (np.linalg.norm(J_fd) 1e-14) max_err max(max_err, err) print(fJacobian 最大相对误差: {max_err:.3e}) return max_err如果最大相对误差在 1e-5 量级Jacobian 基本正确。如果误差在 1e-2 以上说明对流项的导数推导或组装有误。这个检查花不了几分钟但能省掉几个小时的盲目调试。我做了这么多年有限元最大的教训是不要相信“看起来对”的代码。压力场没有振荡、速度场光滑不代表 Jacobian 是对的。每次改了单元类型或边界条件跑一遍 Jacobian 检查比事后 debug 划算得多。希望帮到你。本文还有配套的精品资源点击获取