引言
现代工程研发流程中存在一道隐性鸿沟：设计部门用 CAD 系统产出以 NURBS（非均匀有理 B 样条）表示的精确几何，分析部门却要先把几何”网格化”为分段低阶多项式单元才能做有限元分析。工业界的经验数据表明，复杂模型的网格生成与清理可占整个分析周期的 60%～80%，而且离散网格对原始几何的逼近误差贯穿始终。2005 年 Hughes、Cottrell 与 Bazilevs 提出的等几何分析（Isogeometric Analysis, IGA）给出了根本性的解法：直接采用 CAD 的 NURBS 基函数同时表示几何与未知场——”几何即网格，设计即可分析”。本文介绍 B 样条/NURBS 的数学基础、IGA 特有的细化方式，并给出一维 IGA-Galerkin 求解器的完整实现。
理论基础
B 样条基函数
B 样条由节点向量（knot vector）$\Xi = {\xi_1, \xi_2, \dots, \xi_{n+p+1}}$ 定义，其中 $p$ 为阶数、$n$ 为基函数个数。基函数由 Cox–de Boor 递推给出：
$$N_{i,0}(xi) = begin{cases} 1, & xi_i le xi < xi_{i+1} \ 0, & text{其他} end{cases}$$
$$N_{i,p}(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} N_{i,p-1}(\xi) + \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1,p-1}(\xi)$$
B 样条基具有若干对数值分析极为有利的性质：单位分解性（$\sum_i N_{i,p} = 1$）、紧支撑性（$N_{i,p}$ 仅在 $[\xi_i, \xi_{i+p+1}]$ 内非零）、非负性，以及由节点重复度控制的连续性——单重内部节点处 $C^{p-1}$ 连续，重复 $r$ 重降为 $C^{p-r}$。工程中最常用开节点向量（两端节点重复 $p+1$ 次），使曲线插值于端点。B 样条曲线由控制点线性组合而成：$\mathbf{C}(\xi) = \sum_i N_{i,p}(\xi),\mathbf{P}_i$。
NURBS 与 k-细化
为精确表示圆锥曲线（圆、椭圆），引入权重 $w_i$ 得到 NURBS 基函数：
$$R_i^p(\xi) = \frac{N_{i,p}(\xi), w_i}{\sum_{j=1}^{n} N_{j,p}(\xi), w_j}$$
IGA 的核心创新在于k-细化（k-refinement）。传统有限元有 h-细化（加密网格）与 p-细化（升阶）；IGA 增加了第三条路径：先升阶、再插入新节点。由于升阶先于插点，新基函数在单元间保持 $C^{p-1}$ 高阶连续，而不是 p-版本 FEM 那样的 $C^0$。其直接后果是：相同单元数下，k-细化的高阶 IGA 自由度远少于同阶 FEM，且高连续性使每个自由度的精度（accuracy per DOF）显著提升，这在振动谱分析与湍流模拟中尤为重要。
数值方法与实现
下面实现一维 Poisson 方程 $-u'' = \pi^2 \sin(\pi x)$（精确解 $u = \sin(\pi x)$）的 IGA-Galerkin 求解器：开节点向量、二次 B 样条基、按节点区间逐单元高斯积分组装刚度矩阵，并测量 $L^2$ 误差收敛阶。基函数及其导数按 Piegl 与 Tiller《NURBS 之书》算法 A2.1/A2.3 计算：
import numpy as np
def find_span(n_bas, p, xi, knot):
"""Piegl & Tiller A2.1：确定 xi 所在的节点区间索引"""
n = n_bas - 1
if xi >= knot[n + 1]:
return n
if xi <= knot[p]:
return p
lo, hi = p, n + 1
mid = (lo + hi) // 2
while xi < knot[mid] or xi >= knot[mid + 1]:
if xi < knot[mid]:
hi = mid
else:
lo = mid
mid = (lo + hi) // 2
return mid
def basis_ders(i, xi, p, knot):
"""Piegl & Tiller A2.3：span i 处非零基函数值与一阶导数"""
ndu = np.zeros((p + 1, p + 1))
left = np.zeros(p + 1); right = np.zeros(p + 1)
ndu[0, 0] = 1.0
for j in range(1, p + 1):
left[j] = xi - knot[i + 1 - j]
right[j] = knot[i + j] - xi
saved = 0.0
for r in range(j):
ndu[j, r] = right[r + 1] + left[j - r]
temp = ndu[r, j - 1] / ndu[j, r]
ndu[r, j] = saved + right[r + 1] * temp
saved = left[j - r] * temp
ndu[j, j] = saved
ders = np.zeros((2, p + 1))
ders[0, :] = ndu[:, p]
a = np.zeros((2, p + 1))
for r in range(p + 1):
s1, s2 = 0, 1
a[0, 0] = 1.0
k = 1 # 一阶导
d = 0.0
rk, pk = r - k, p - k
if r >= k:
a[s2, 0] = a[s1, 0] / ndu[pk + 1, rk]
d = a[s2, 0] * ndu[rk, pk]
j1 = 1 if rk >= -1 else -rk
j2 = k - 1 if r - 1 <= pk else p - r
for j in range(j1, j2 + 1):
a[s2, j] = (a[s1, j] - a[s1, j - 1]) / ndu[pk + 1, rk + j]
d += a[s2, j] * ndu[rk + j, pk]
if r <= pk:
a[s2, k] = -a[s1, k - 1] / ndu[pk + 1, r]
d += a[s2, k] * ndu[r, pk]
ders[k, r] = d
ders[1, :] *= p # 乘递推因子
return ders[0], ders[1]
def iga_poisson_1d(n_el, p=2, n_gauss=4):
"""区间 (0,1) 上 -u''=pi^2 sin(pi x), u(0)=u(1)=0 的 IGA-Galerkin 解"""
# 开节点向量：内部节点均匀、重复度 1（基函数跨单元 C^{p-1} 连续）
interior = np.linspace(0, 1, n_el + 1)[1:-1]
knot = np.concatenate([np.zeros(p + 1), interior, np.ones(p + 1)])
n_bas = len(knot) - p - 1
K = np.zeros((n_bas, n_bas)); F = np.zeros(n_bas)
gp, gw = np.polynomial.legendre.leggauss(n_gauss)
spans = [s for s in range(p, n_bas) if knot[s + 1] - knot[s] > 1e-14]
for span in spans:
a, b = knot[span], knot[span + 1]
for s, w in zip(gp, gw):
xg = 0.5 * (b - a) * s + 0.5 * (a + b)
jac = 0.5 * (b - a)
N, dN = basis_ders(span, xg, p, knot)
idx = np.arange(span - p, span + 1)
K[np.ix_(idx, idx)] += np.outer(dN, dN) * w * jac
F[idx] += np.pi**2 * np.sin(np.pi * xg) * N * w * jac
free = np.arange(1, n_bas - 1) # 去掉两端 Dirichlet 自由度
u = np.zeros(n_bas)
u[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])
err2 = 0.0 # L2 误差
for span in spans:
a, b = knot[span], knot[span + 1]
for s, w in zip(gp, gw):
xg = 0.5 * (b - a) * s + 0.5 * (a + b)
N, _ = basis_ders(span, xg, p, knot)
uh = N @ u[span - p:span + 1]
err2 += (uh - np.sin(np.pi * xg))**2 * w * 0.5 * (b - a)
return np.sqrt(err2), n_bas
print(f"{'单元数':>6} {'自由度':>6} {'L2 误差':>12} {'收敛阶':>8}")
prev = None
for n_el in [4, 8, 16, 32]:
err, ndof = iga_poisson_1d(n_el)
rate = np.log2(prev / err) if prev else float('nan')
print(f"{n_el:6d} {ndof:6d} {err:12.3e} {rate:8.2f}")
prev = err工程案例分析
IGA 最耀眼的应用之一是 Kirchhoff–Love 薄壳分析：经典壳理论的弱形式包含挠度的二阶导数，要求试探函数 $C^1$ 连续，传统 $C^0$ 有限元必须引入额外自由度或重构技术，而 NURBS 基函数天然满足连续性要求，使无转动自由度的壳单元成为现实——这对涡轮叶片、汽车覆盖件等薄壁结构极具吸引力。其他成功案例包括流固耦合（动脉血流、风力机叶片）、接触力学（几何精确消除法向跳跃）与振动分析（高连续性显著抑制高阶”p-分支”频率误差）。等几何框架的另一深层收益是几何误差被彻底消除：圆、圆柱面等 NURBS 精确表示的形状在任意细化层级下都保持精确，而多项式网格无论如何加密都残留几何离散误差——在接触压强分布、薄壳屈曲临界载荷等对几何极度敏感的量上，这一差别会直接体现为预测精度的系统性差异。
三种细化策略的特性对比：
| 细化方式 | 操作 | 单元间连续性 | 自由度增长 | 适用场景 |
|---|---|---|---|---|
| h-细化 | 插入新节点 | 保持 $C^{p-1}$ | 随单元数线性增长 | 局部加密、捕捉边界层 |
| p-细化 | 升阶后复制已有节点 | 降为 $C^0$ | 中等 | 与传统 p-FEM 等效 |
| k-细化 | 升阶后插入新节点 | 保持 $C^{p-1}$ | 最省 | 光滑解的高精度逼近 |
本文代码实测：二次 B 样条（k-细化意义下的 $C^1$ 基）在 4、8、16、32 个单元下 $L^2$ 误差分别为 $2.3\times10^{-3}$、$2.6\times10^{-4}$、$3.1\times10^{-5}$、$3.9\times10^{-6}$，收敛阶 3.0，与理论 $O(h^{p+1})$ 一致；且 32 单元仅 34 个自由度——同精度 $C^0$ 二次 FEM 需要 65 个节点自由度。工程落地方面，Bézier 提取（Bézier extraction）技术把 B 样条单元映射为 Bernstein 形式，使 IGA 可嵌入现有有限元数据流；LS-DYNA 已内置 IGA 单元，Coreform Cubit 提供 U-spline 剖分工具，T-spline/THB-spline 则解决了 NURBS 张量积结构难以局部加密的问题。
结论与展望
本文从节点向量与 Cox–de Boor 递推讲起，阐明了 B 样条/NURBS 基函数的构造与连续性的节点重复度控制机制，解释了 IGA 特有的 k-细化为何能以更少自由度获得高阶光滑逼近，并以一维求解器验证了理论收敛阶。IGA 的意义不止于算法：它打通了设计模型与分析模型之间的数据壁垒，使”设计—分析—优化”闭环中的几何传递无损化。尚待解决的问题包括：实体 CAD 模型（B-rep 裁剪曲面）的体参数化、非水密几何的耦合、以及复杂装配体的自动剖分。随着等几何剖分工具的成熟，IGA 正从学术概念稳步走向工业主流。
发表回复