几何表示法研究（二十三）

第 11 章 有限元网格与等几何分析 IGA

本章摘要：本章围绕"为计算而生"的有限元网格展开。先厘清它与渲染用多边形网格的本质区别，并点出 CAD（NURBS）与仿真（线性网格）之间的“几何→网格”鸿沟；随后给出有限元网格与等几何分析（IGA）的定义，梳理常见单元类型；接着以泊松方程为例，从强形式推导弱形式（Galerkin 框架），引出刚度矩阵与 Strang 引理，并介绍等参映射、高斯求积及 IGA 用 NURBS 基同时充当几何与形函数的核心思想；再给出从 B-Rep 到 FEA 网格的生成流水线、IGA 实现要点与避坑清单；最后通过 FreeCAD 实例演示完整链路，并在小结中归纳要点、适用场景与常见坑。

1. 背景与动机

我们在 第 8 章 学过多边形网格——它是给 GPU 看的，目的是把形状画出来。但工程里还有一大类需求是算，而不是画：应力多大、会不会断、温度怎么分布、流体怎么绕流……这些都要求解偏微分方程（PDE）。

有限元网格（FEA 网格）就是“为计算而生”的几何表示：它和多边形网格长得像，但目的根本不同——后者为渲染，前者为求解。一个有限元网格把求解区域剖分成许多单元（element），在每个单元上用低次多项式近似未知场（位移、温度、电势……），再把所有单元“拼起来满足弱形式”逼近真实解。

这里还藏着一个行业痛点：CAD 用 NURBS（见第 3 章），仿真用线性网格，两者之间要断开重画一遍——这条“几何→网格”的鸿沟既费时又丢精度。2005 年 Hughes 等人提出的**等几何分析（IGA）**想一劳永逸地解决它。本章你就搞懂：有限元网格在算什么、怎么从 CAD 生成它、以及 IGA 为何是“几何即网格”。

2. 核心概念与原理

定义（有限元网格）：把求解域 $\Omega$ 剖分为有限个单元的集合；在每个单元上用分片多项式近似未知场 $u(\mathbf{x})$，再整体满足控制方程的弱形式。

常见单元（按维度）：

下面这张对比图直观展示了两种表示方式在几何保真度上的差异——左侧是经典 FEM 的线性三角形网格，用折边逼近光滑曲面；右侧是 IGA 的 NURBS 控制网格，直接以光滑曲面为几何：

几何保真度差异：经典 FEM 用分片线性三角形去逼近光滑曲面，单元越粗折边越明显，几何信息在"几何→网格"转换中有损；而 IGA 直接用 NURBS 控制点网格描述几何，曲面本身无损保留，控制点既是几何参数又是形函数系数，无需重新剖分。

关键判别：六面体（Hex）精度高、收敛快，但自动生成极难；四面体（Tet）易自动生成（Delaunay/前沿/八叉树），精度通常略逊，所以工程常用混合网格。

定义（等几何分析 IGA）：直接用 CAD 的 NURBS 基函数同时充当几何与有限元形函数，让"几何就是网格、网格就是几何"，消除几何→网格的转换误差。

3. 关键公式与推导

3.1 从强形式到弱形式（Galerkin 框架）

以泊松方程（椭圆边值问题）为例，强形式为：

$$-\nabla \cdot (k\nabla u)=f\text{in }\Omega ,u{|}_{\partial \Omega }=0$$

推导弱形式：两边乘以测试函数 $v$（满足齐次边界 $v{|}_{\partial \Omega }=0$），在 $\Omega$ 上积分：

$$-{\int }_{\Omega }v\nabla \cdot (k\nabla u)d\Omega ={\int }_{\Omega }fvd\Omega$$

对左端用分部积分（散度定理）：

$$-{\int }_{\Omega }v\nabla \cdot (k\nabla u)d\Omega =-{\int }_{\partial \Omega }v(k\nabla u)\cdot \mathbf{n}dS+{\int }_{\Omega }k\nabla u\cdot \nabla vd\Omega$$

由于 $v=0$ 在边界上，面积分项消失，得到弱形式：

$${\int }_{\Omega }k\nabla u\cdot \nabla vd\Omega ={\int }_{\Omega }fvd\Omega ,\forall v\in H_0^1(\Omega )$$

代入有限元近似：令 $u\approx {\sum }_i{u}_i{N}_i(\mathbf{x})$，并取测试函数 $v=N_j$，则

$$\sum _i{u}_i{\int }_{\Omega }k\nabla N_i\cdot \nabla N_{j}d\Omega ={\int }_{\Omega }f{N}_{j}d\Omega$$

写成线性系统 $\mathbf{Ku}=\mathbf{F}$，其中

$$K_{ij}={\int }_{\Omega }k\nabla N_i\cdot \nabla N_{j}d\Omega ,F_j={\int }_{\Omega }f{N}_{j}d\Omega$$

这就是刚度矩阵 $\mathbf{K}$ 的来源。收敛性由 Strang 引理 保证：总误差 = 逼近误差（单元次数 $p$） + 积分/几何误差（网格尺度 $h$），于是有了 $h$-加密、$p$-加密、$hp$-自适应的三大策略。

下面用一段 Python 代码，把一维杆单元（$u^{\prime \prime }=0$，两端点位移已知）的线性形函数与单元刚度矩阵完整推导并计算出来，验证它恰好等于弹簧刚度 $k/L$：

# 一维杆单元：u'' = 0，两端点位移 u0, u1 已知
# 单元长度 L，材料参数 k（这里 k = E*A，E 弹性模量，A 截面积）
import numpy as np

L = 2.0          # 单元长度
k = 100.0        # 轴向刚度 E*A

# 1) 线性形函数（自然坐标 xi ∈ [-1, 1]）
#    N0(xi) = (1 - xi)/2,  N1(xi) = (1 + xi)/2
#    在物理坐标 x ∈ [0, L] 下，xi = 2x/L - 1，代入得：
#    N0(x) = 1 - x/L,  N1(x) = x/L
def N0(x): return 1.0 - x / L
def N1(x): return x / L

# 2) 形函数对物理坐标的导数（B 矩阵）
#    dN0/dx = -1/L,  dN1/dx = +1/L
dN0 = -1.0 / L
dN1 = +1.0 / L
B = np.array([[dN0, dN1]])          # 形状 (1, 2)

# 3) 单元刚度矩阵：K^e = ∫_0^L B^T k B dx
#    被积函数为常数（B 与 x 无关），积分 = 长度 L 乘以被积函数
Ke = L * (B.T @ (k * B))            # 形状 (2, 2)
print("单元刚度矩阵 K^e =")
print(Ke)

# 4) 验证：结果恰好等于弹簧刚度 k/L
#    K^e = (k/L) * [[ 1, -1],
#                   [-1,  1]]
k_spring = k / L
expected = k_spring * np.array([[ 1.0, -1.0],
                                [-1.0,  1.0]])
print("\n弹簧刚度 k/L =", k_spring)
print("与弹簧刚度矩阵一致:", np.allclose(Ke, expected))

# 5) 组装并求解：给定两端位移，求节点力
u = np.array([0.0, 0.5])            # 左端固定，右端位移 0.5
F = Ke @ u
print("\n节点力 F = K^e @ u =", F)
# 物理意义：F1 = -k/L * u1（反力），F2 = +k/L * u1（外力），
# 与弹簧 F = k_spring * Δu 完全一致。

要点说明：一维杆单元的形函数 $N_0=1-x/L$、$N_1=x/L$ 是线性插值，其导数 $B=[-1/L,+1/L]$ 为常数，因此单元刚度矩阵 $K^e={\int }_0^L{B}^{T}kBdx=\frac{k}{L}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$，恰好就是弹簧刚度 $k/L$ 的矩阵形式——这正是“杆单元等价于弹簧”的数学来源，也呼应了本章练习第 1 题。

3.2 等参映射与刚度矩阵积分

物理单元与参考单元 $\hat{K}$ 之间用等参映射 $\mathbf{x}=F_K(\hat{\xi })={\sum }_i{x}_i^{node}{N}_i(\hat{\xi })$，其雅可比为

$$J_K=\frac{\partial \mathbf{x}}{\partial \hat{\xi }},d\mathbf{x}=|J_K|d\hat{\xi }$$

单元刚度矩阵用高斯求积近似（在参考单元取若干积分点 ${\xi }_q$、权重 $w_q$）：

$$K_{ij}^{(e)}=\sum _q{w}_q\nabla N_i({\xi }_q{)}^{T}k\nabla N_j({\xi }_q)|J_K({\xi }_q)|$$

注意：若雅可比行列式 $|J_K|\le 0$（退化/反转单元），积分会失效——这正是生成网格必须保证 $|J|>0$ 的原因（见第四节避坑）。

3.3 IGA：用 NURBS 基同时作几何与形函数

经典 FEM 的形函数 $N_i$ 是分片 Lagrange 多项式，而几何是 NURBS——表示不匹配。IGA 让形函数直接取自 CAD 的 B-spline/NURBS 基：

$$\text{几何：}x(\xi )=\sum _i{P}_i{B}_i^p(\xi ),\text{场：}u(\xi )=\sum _i{u}_i{B}_i^p(\xi )$$

同一组基！于是几何精度无损失、高阶连续性 $C^{p-1}$ 天然融入、细化不丢几何。IGA 的基函数就是第 3 章的 Cox–de Boor 递推：

$$B_i^p(u)=\frac{u-u_i}{u_{i+p}-u_i}{B}_i^{p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{B}_{i+1}^{p-1}(u)$$

则位移场 $u(\xi )={\sum }_i{u}_i{B}_i^p(\xi )$，空间梯度 ${\nabla }_{x}u=J^{-1}{\nabla }_{\xi }u$，$J={\sum }_i{x}_i(B_i^p{)}^{\prime }(\xi )$。

下表从几个关键维度对比经典 FEM 与 IGA，帮你快速抓住两者的本质差异：

点评：经典 FEM 胜在通用与成熟，但“几何→网格”的转换既费时又丢精度；IGA 用同一组 NURBS 基同时描述几何与场，换来几何无损与高阶连续性，代价是单单元计算量更大、生态尚不成熟。工程上两者并非替代关系——通用仿真仍以 FEM 为主，而薄壳、形状优化等对几何精度敏感的场景，IGA 的优势才真正凸显。

4. 具体实现步骤

4.1 网格生成流水线（从 B-Rep 到 FEA 网格）

1. 几何导入与清理：从 B-Rep（第 1 章）取面/边，修补缝隙、去除小特征（defeature）。
2. 表面网格：在参数域或直接在曲面上生成三角/四边形（前沿推进、Delaunay）。
3. 体网格生成：

• Delaunay / 约束 Delaunay（CDT）：最大化最小角，易实现、稳健，但边界质量差；
• 前沿推进（Advancing Front）：边界质量好但易自交；
• 八叉树（第 7 章）：稳健、并行友好（代表库：CGAL、TetGen、gmsh、MMG、TetWild）。

4. 质量优化：smoothing（Laplacian/优化）、边交换、重划分（remeshing）。
5. 各向异性自适应：基于误差估计调整度量场后重生成（边界层/应力集中处需要极扁的拉伸单元）。

4.2 IGA 实现要点

• 复用 CAD 内核（OCCT）的 NURBS 求值、求导、求积；
• 处理多片（multipatch）拼接：片间需 $C^0$（位移连续）或 $C^1$（法向连续），常用绑定（tying）/ Nitsche 法 / mortar 法；
• 修剪（trimmed）曲面仍是难题，需 THB-splines、LR B-splines 等局部细化参数化。

4.3 避坑清单

1. 雅可比为正：生成器必须保证每个单元 $|J|>0$，否则积分崩塌。
2. 边界层网格：CFD / 薄边界层需极扁单元，对求解器条件数极敏感。
3. 稀疏装配：刚度矩阵高度稀疏，用 CRS / CSR 存储 + 迭代求解器（AMG、CG、GMRES）。
4. 保形网格：多物理场耦合需一致网格，否则插值误差。
5. IGA 成本高：高阶基 + 密集高斯点，单单元计算量远大于线性 FEM，需权衡精度收益。

5. 实际应用示例

目标：在 FreeCAD 里走通“草图 → B-Rep 实体 → 网格（STL）→ FEA 网格”的完整链路，对应 第 14 章 的综合流程。

# FreeCAD Python 控制台（需已打开一个 PartDesign 实体）
import Mesh, Fem
# 1) B-Rep 实体已存在于 ActiveObject（如一个 Pad）
shape = FreeCAD.ActiveDocument.ActiveObject.Shape
# 2) 把 B-Rep 分面成网格（tessellation，弦差 0.1）
mesh = Mesh.Mesh()
mesh.addMesh(shape.tessellate(0.1))     # 返回 (nodes, facets)
Mesh.show(mesh)                          # 在视图显示三角网格
# 3) 交给 FEM 生成仿真网格（gmsh 后端）
from femmesh import gmshtools
gm = gmshtools.GmshTools()
gm.create_mesh(FreeCAD.ActiveDocument.ActiveObject,
               mesh_size=1.0, order=2)   # 二阶单元
print("FEA 单元数:", gm.mesh_object.FemMesh.VolumeCount)

这条链路的每一环都对应源码里的模块：草图在 Sketcher、实体在 PartDesign（特征树）→ Part::TopoShape（B-Rep，见 第 1 章）→ Mesh（分面）→ Fem/femmesh（调 gmsh/MMG 生成仿真网格，见上文 4.1）。对医学 CT 数据，则常见“体素 → 面网格 → Tet 网格”的路线（结合 第 6、7 章）。

上面的链路走的是“B-Rep → 三角网格 → FEA 网格”的传统路线。而 IGA 的思路是直接复用 CAD 的 NURBS 基，让几何本身就是网格。下面用 geomdl 库演示：构建一个双三次 B 样条曲面，把它的控制点网格直接当作 IGA 的初始网格，并说明如何把控制点坐标映射为形函数系数。

# pip install geomdl
from geomdl import BSpline, utilities

# 1) 构建一个 3 x 3 控制点的双三次 B 样条曲面（次数 p = q = 3）
surf = BSpline.Surface()
surf.degree_u = 3
surf.degree_v = 3
# 控制点网格：每个点 (x, y, z)，这里做一个轻微起伏的曲面
surf.ctrlpts2d = [
    [[0, 0, 0], [1, 0, 0.2], [2, 0, 0], [3, 0, 0.1]],
    [[0, 1, 0.3], [1, 1, 0.5], [2, 1, 0.4], [3, 1, 0.6]],
    [[0, 2, 0], [1, 2, 0.2], [2, 2, 0], [3, 2, 0.1]],
    [[0, 3, 0.1], [1, 3, 0.3], [2, 3, 0.2], [3, 3, 0.4]],
]
# 均匀节点向量（clamped，首尾重复 p+1 次）
surf.knotvector_u = utilities.generate_knot_vector(3, len(surf.ctrlpts2d))
surf.knotvector_v = utilities.generate_knot_vector(3, len(surf.ctrlpts2d[0]))

# 2) 提取控制点网格 —— 这就是 IGA 的“初始网格”
#    每个控制点 P_ij 对应一个基函数 B_i^p(u) * B_j^q(v)
ctrl_grid = surf.ctrlpts2d          # 形状 (nu+1) x (nv+1) x 3
nu, nv = len(ctrl_grid), len(ctrl_grid[0])
print(f"IGA 初始网格：{nu} x {nv} 个控制点（即 {nu*nv} 个形函数）")

# 3) 关键：控制点坐标 -> 形函数系数
#    在 IGA 中，几何与场共用同一组基：
#       几何:  x(u,v) = Σ_i Σ_j P_ij * B_i^p(u) * B_j^q(v)
#       场:    u(u,v) = Σ_i Σ_j u_ij * B_i^p(u) * B_j^q(v)
#    因此"控制点坐标 P_ij"就是几何的系数；而位移/温度等场量
#    的系数 u_ij 则存放在与 ctrlpts2d 同构的数组里，逐点对应。
#    下面把控制点坐标直接作为初始场系数（例如初始位移=几何坐标）：
field_coeffs = [[[p[0], p[1], p[2]] for p in row] for row in ctrl_grid]
print("形函数系数数组形状:", len(field_coeffs), "x", len(field_coeffs[0]))

# 4) 在参数域采样求值，验证曲面（可选）
#    surf.evaluate() 会按 knotvector 在 (u,v) 网格上求值
#    pts = surf.evaluate(lists=True)   # 形状 (nu_sample, nv_sample, 3)

要点说明：ctrlpts2d 的每个点 $P_{ij}$ 就是 IGA 里的一个"节点"（控制点），它同时承载几何坐标与场系数；generate_knot_vector 生成的 clamped 节点向量保证曲面精确插值边界。实际 IGA 求解时，只需把刚度矩阵 $K_{ij}$ 的积分换成对 B 样条基函数的高斯求积（见 3.2 节），自由度就是这些控制点系数，无需再生成三角网格。

六、小结与要点

• 记住什么：有限元网格是为求解 PDE 而生；弱形式由强形式乘以测试函数再分部积分得到；刚度矩阵来自 $\int \nabla N_i\cdot \nabla N_j$；IGA 用 NURBS 基同时当几何与形函数。
• 适用场景：结构力学（应力/模态/疲劳）、CFD、热电磁、生物力学；IGA 用于薄壳、形状优化、设计-分析一体化。
• 常见坑：① 经典 FEM 把光滑 CAD 压成折边，几何保真度低；② 六面体自动生成是开放难题；③ 上游 B-Rep 缝隙/小特征会毁掉网格；④ 高阶/密集单元计算量大；⑤ IGA 生态不成熟（多片拼接、修剪、局部细化难）；⑥ 它是“分析表示”不是“几何表示”，不能直接用于制造，必须与 B-Rep/网格并存。

练习与思考

1. 对一维杆单元（$u^{\prime \prime }=0$，两端点位移已知），手推其线性形函数 $N_0,N_1$ 与单元刚度矩阵，验证它恰好等于“弹簧刚度”。
2. 为什么 IGA 中片间连接通常用 Nitsche 法而不是简单“共享节点”？查阅其如何处理 $C^1$ 连续性。
3. 查阅 gmsh 的 Mesh.SizeFromCurvature 参数，思考“边界曲率大的地方自动加密”在工程上解决了什么问题。