企业官网建设流程全解析

简介：这份资源是面向流体力学、数值计算方向的学习者与工程师的二维稳态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-Brezzi（LBB）条件，也叫 inf-sup 条件：

inf_{q_h ∈ Q_h} sup_{v_h ∈ V_h} |∫ q_h (∇·v_h) dΩ| / (||v_h||_1 ||q_h||_0) ≥ β > 0

β 与网格尺寸无关。如果违反这个条件，压力会出现棋盘格振荡，而且加密网格也不会改善。常见的单元组合对比如下：

我一般首选 Taylor-Hood（P2-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*i+1, 2*j+1] = 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*i+1, 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_iter=20, tol=1e-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(f"Iter {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/s），tol 要相应放大到 1e-6 左右。

3.3 雷诺数升高时 Newton 为什么不收敛

这是稳态 NS 有限元最常翻车的地方。雷诺数超过某个临界值后，稳态解可能根本不存在（流动本质是瞬态的），或者存在但 Newton 的收敛域很窄。我踩过的坑是：Re=100 的方腔驱动流，用 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)

另一个技巧是雷诺数延拓：先算 Re=10 的解，把它作为 Re=50 的初始猜测，再作为 Re=100 的初始猜测。每次增加的雷诺数不要超过 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 设得太小，低于离散误差的底线。

解决：先检查网格是否足够细。对于方腔驱动流 Re=100，至少需要 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 有限元最经典的验证案例：单位正方形区域，顶边速度 u=1，其他三边无滑移，Re=100 到 1000。Ghia 等人在 1982 年给出了高精度基准解，常用的是沿垂直中线和水平中线的速度剖面。

验证步骤：在 64×64 的 Taylor-Hood 网格上算 Re=100，提取 x=0.5 处的 u 速度沿 y 方向的分布，和 Ghia 数据对比。如果最大偏差在 2% 以内，求解器基本可信。Re=1000 时需要更细的网格（128×128 以上），而且 Newton 可能需要延拓策略。

# 提取 x=0.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]

这个函数返回 x=0.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, eps=1e-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(f"Jacobian 最大相对误差: {max_err:.3e}") return max_err

如果最大相对误差在 1e-5 量级，Jacobian 基本正确。如果误差在 1e-2 以上，说明对流项的导数推导或组装有误。这个检查花不了几分钟，但能省掉几个小时的盲目调试。

我做了这么多年有限元，最大的教训是：不要相信“看起来对”的代码。压力场没有振荡、速度场光滑，不代表 Jacobian 是对的。每次改了单元类型或边界条件，跑一遍 Jacobian 检查，比事后 debug 划算得多。希望帮到你。

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