基于等几何分析的泊松问题求解

2026/05/20 14:00:00

课程等几何分析·16 min read

等几何分析泊松方程 Galerkin NURBS 课程笔记

GAMES302 · 高级图形学与增强现实

授课教师：徐岗 · 杭州电子科技大学

引言：为什么用 IGA 求解偏微分方程

传统有限元方法（FEM）求解偏微分方程（PDE）时，几何表示（网格）与物理场近似（基函数）通常来自两套相互独立的模型——这导致两个根本性问题：

几何精度受限于网格分辨率：加密网格可以提升精度，但网格生成本身成本极高。

几何-物理不一致误差：CAD 几何与离散网格之间的差异会在边界条件处理上引入额外误差。

等几何分析（Isogeometric Analysis, IGA） 提出了一个统一方案：用 NURBS 基函数同时表示几何体和物理场——几何建模用 NURBS，物理场近似也用相同的 NURBS 基函数展开。几何与物理共享同一种数学表示，彻底消除了两套模型之间的不一致性。

本讲以泊松方程（Poisson Equation）——热传导方程的稳态形式——为模型问题，完整推导 IGA 求解流程，包括强形式 → 弱形式 → Galerkin 离散 → 矩阵装配 → 代码实现，最后介绍 $$h$$ 、 $$p$$ 、 $$r$$ 细化策略。

泊松方程的强形式

设区域 $\Omega \subset \mathbb{R}^2$ 由 NURBS 曲面定义，边界 $\Gamma = \partial\Omega$ 分三类：

$\Gamma_D$ （迪利克雷边界）：温度场 $$u$$ 固定， $u|_{\Gamma_D} = g$

$\Gamma_N$ （纽曼边界）：热通量固定， $K\nabla u \cdot n|_{\Gamma_N} = h$

$\Gamma_R$ （罗宾边界）：对流换热条件， $K\nabla u \cdot n + \beta(u - U_\infty) = 0$

强形式为：找到 $u : \Omega \to \mathbb{R}$ ，使得

\begin{cases} -\nabla \cdot (K \nabla u) = f, & \text{in } \Omega \\ u = g, & \text{on } \Gamma_D \\ K\nabla u \cdot n = h, & \text{on } \Gamma_N \\ K\nabla u \cdot n + \beta(u - U_\infty) = 0, & \text{on } \Gamma_R \end{cases}

其中：

$$K > 0$$ ：热传导率（标量或张量）

$$f$$ ：热源密度函数

$g, h, \beta, U_\infty$ ：边界数据，均为已知

从强形式到弱形式

1 试探函数空间与测试函数空间

定义试探空间（trial space）：

\Phi = \{ v \in H^1(\Omega) \mid v = g \text{ on } \Gamma_D \}

其中 $H^1(\Omega)$ 是索伯列夫空间，即一阶偏导数平方可积的函数空间。

定义测试空间（test space）：

\Psi = \{ w \in H^1(\Omega) \mid w = 0 \text{ on } \Gamma_D \}

注意测试函数在迪利克雷边界上必须强制为零，这是弱形式推导的要点。

2 乘以测试函数并分部积分

将 PDE 两边同乘 $w \in \Psi$ ，在 $\Omega$ 上积分：

\int_\Omega w(-\nabla \cdot (K\nabla u)) d\Omega = \int_\Omega wf d\Omega

对左端做分部积分（Green 第一公式）：

\int_\Omega K\nabla w \cdot \nabla u \, d\Omega - \int_{\partial\Omega} w K\nabla u \cdot n \, d\Gamma = \int_\Omega wf d\Omega

3 代入三类边界条件

对三类边界分别处理：

罗宾边界（ $\Gamma_R$ ）：代入 $K\nabla u \cdot n = -\beta(u - U_\infty)$ ，得到

\int_{\Gamma_R} w \cdot [-\beta(u - U_\infty)] d\Gamma = -\beta \int_{\Gamma_R} w u \, d\Gamma + \beta U_\infty \int_{\Gamma_R} w \, d\Gamma

纽曼边界（ $\Gamma_N$ ）：已知 $K\nabla u \cdot n = h$ ，保留

\int_{\Gamma_N} w h \, d\Gamma

迪利克雷边界（ $\Gamma_D$ ）： $$w = 0$$ ，自动消除

4 最终弱形式

整理后得到标准的变分形式：找到 $u \in \Phi$ ，使得

a(w, u) = L(w), \quad \forall w \in \Psi

其中

a(w, u) = \int_\Omega K \nabla w \cdot \nabla u \, d\Omega + \beta \int_{\Gamma_R} w u \, d\Gamma

L(w) = \int_\Omega w f \, d\Omega + \int_{\Gamma_N} w h \, d\Gamma + \beta U_\infty \int_{\Gamma_R} w \, d\Gamma

弱形式的核心价值：导数阶次降低一阶——原方程含二阶导，弱形式只要求 $u \in H^1$ ，允许分片多项式作为基函数。

Galerkin 方法：将无限维转为有限维

1 有限维子空间构造

在 IGA 中，NURBS 基函数同时构成试探空间和测试空间的离散近似：

\Phi_h = \text{span}\{N_i\}_{i=1}^{n}, \quad \Psi_h = \text{span}\{N_i\}_{i=1}^{n_{eq}}

其中 $n_{eq} < n$ ，体现迪利克雷边界条件的约束。

Galerkin 方法要求近似解 $$u_h$$ 满足：

u_h = v_h + g_h, \quad v_h \in \Psi_h, \ g_h \text{ 逼近边界条件}

将 $u_h = \sum_{i=1}^{n} c_i N_i$ 代入弱形式，利用 $$v_h$$ 的任意性，得到线性系统：

\underbrace{\mathbf{K}}_{n_{eq} \times n_{eq}} \underbrace{\mathbf{d}}_{n_{eq} \times 1} = \underbrace{\mathbf{F}}_{n_{eq} \times 1}

其中：

$\mathbf{K}$ ：刚度矩阵（stiffness matrix）， $K_{ij} = a(N_j, N_i)$

$\mathbf{d}$ ：位移/温度矢量（待求未知数）

$\mathbf{F}$ ：力矢量（load vector）， $$F_i = L(N_i)$$

2 边界条件的处理策略

由于 NURBS 基函数是高度局部化的，在迪利克雷边界 $\Gamma_D$ 上非零的基函数极少。因此存在整数 $n_{eq} < n$ ，使得前 $n_{eq}$ 个基函数在 $\Gamma_D$ 上满足齐次边界条件（ $$g_i = 0$$ ）。

对每个基函数引入 BC 函数：

\text{BC}[i] = \begin{cases} 1, & N_i \text{ 在 } \Gamma_D \text{ 上非零} \\ 0, & \text{其他} \end{cases}

装配全局矩阵时：

若 $\text{BC}[i] = 1$ ：将第 $$i$$ 行/列化为单位矩阵， $$d_i = g_i$$ （直接赋值）

若 $\text{BC}[i] = 0$ 但 $\text{BC}[j] = 1$ ：不更新 $K_{ij}$ ，将 $$F_i$$ 减去 $K_{ij} g_j$

单元刚度矩阵与力向量的构造

1 参考域到物理域的映射

将物理域单元 $\Omega^e \subset \Omega$ 通过映射 $(x^e(\xi), y^e(\eta))$ 映射到参考单元 $\hat{\Omega} = [0,1]^2$ ，引入以下关键量：

雅可比矩阵： $DF = \begin{pmatrix} x_\xi & y_\xi \\ x_\eta & y_\eta \end{pmatrix}$

雅可比行列式： $J = \det(DF) = x_\xi y_\eta - x_\eta y_\xi$

单位外法向量 $\mathbf{n}$ 、单位切向量 $\mathbf{t}$

表面行列式： $J_s = |\mathbf{n} \cdot (x_\xi, y_\xi) \times (x_\eta, y_\eta)|$

物理域积分可转化为参考域积分：

\int_{\Omega^e} \nabla w \cdot \nabla u \, d\Omega = \int_{\hat\Omega} (DF^{-T} \hat\nabla w)^T (DF^{-T} \hat\nabla u) |J| \, d\hat\Omega

2 单元矩阵与单元力向量表达式

K^e_{ab} = \int_{\hat\Omega} (B_a)^T D B_b \cdot |J| \, d\hat\xi d\hat\eta

f^e_a = \int_{\Omega^e} N_a f \, d\Omega = \int_{\hat\Omega} R_a \cdot f(x(\xi), y(\eta)) \cdot |J| \, d\hat\xi d\hat\eta

其中 $$B_a$$ 是 $$N_a$$ 的梯度在参考域中的表示， $$D$$ 是材料属性（这里 $$D = K$$ ）。

3 基函数求值与高斯积分

问题 1：如何计算基函数 $$N_i$$ 及其导数？

高斯积分点与权系数

利用 B 样条基函数与 Bernstein 基函数的关系，通过链式求导法则：

\frac{\partial N_i}{\partial \xi} = \sum_k \frac{\partial B_k}{\partial \xi} \cdot \text{weights}_k, \quad \frac{\partial N_i}{\partial x} = \frac{\partial N_i}{\partial \xi} \cdot \frac{\partial \xi}{\partial x}

问题 2：如何在参考单元上计算积分？

使用高斯-勒让德数值积分。对于 $$M$$ 阶积分精度，取 $$N$$ 个高斯点 $\{\xi_i\}$ 和权系数 $\{w_i\}$ ：

\int_{-1}^{1} f(\xi) d\xi \approx \sum_{i=1}^{N} f(\xi_i) w_i

常用配置表：

阶次 $$M$$	精度 $$N$$	高斯点位置 $\xi_i$	权系数 $$w_i$$
2	3	$\pm 0.577$	$$1$$
3	5	$0, \pm 0.775$	$$0.889, 0.556$$
4	7	$\pm 0.340, \pm 0.861$	$$0.652, 0.348$$

代码实现框架（C++ / Eigen）

IGA 求解泊松问题的整体流程：

`cpp

// 1. 初始化几何数据（控制点、节点矢量、权因子）

initial();

// 2. 构造单元连接矩阵（IEN 数组）

buildConn();

// 3. 装配全局刚度矩阵 A 和全局力向量 b

this->A = assembleA();

this->b = assembleb();

// 4. 处理迪利克雷边界条件

exeBC();

// 5. 稀疏矩阵求解（Cholesky 分解）

Eigen::SimplicialCholesky chol(A);

this->x = chol.solve(b);

assembleA 的核心：遍历所有单元，对每个高斯积分点调用 getIntegralA()，计算 $K^e_{ab}$ ，根据 IEN 数组更新全局矩阵 $\mathbf{A}$ 。

exeBC 的核心：遍历边界上的控制点，将对应行/列清零后将对角线置 1，右端项设为边界值（如 $$g = 0$$ 对应齐次迪利克雷条件）。

等几何分析中的

$h$

、

$p$

、

$r$

细化

1 $$h$$ -细化（knot insertion）

不加新控制点，不改变几何形状，仅在节点矢量中插入新的节点值。插入后基函数支撑域变小，网格局部加密，整体自由度增加。

2 $$p$$ -细化（order elevation）

不改变几何形状，将 B 样条基函数的阶数提升（一次提升一个）。阶数升高后，基函数的表达能力增强，在不加密网格的情况下提高精度。

3 $$r$$ -细化（本讲重点）

$$h$$ -细化和 $$p$$ -细化都通过增加自由度来提升精度。 $$r$$ -细化则是固定自由度数量，只优化内部控制点的位置：

$r$

-细化方法示意

给定初始控制点布局，在边界几何固定的前提下，寻找最优的内部控制点位置，使 IGA 求解误差最小。

这一思想源自形状优化（shape optimization）——区别在于形状优化调整边界控制点， $$r$$ -细化只调整内部控制点，因此不改变边界几何。

3.1 残差函数与后验误差估计

设 $$u_h$$ 为 IGA 解， $$U$$ 为精确解，误差 $$e = U - u_h$$ 。残差函数为：

r = f - \Delta u_h

后验误差指示子（K 单元上）为：

\eta_K = h_K \| r \|_{L^2(K)}

$$J$$ （残差能量模的平方）为：

J = \sum_{K \in \mathcal{T}_h} h_K^2 \int_K (f - \Delta u_h)^2 dK

3.2 $$r$$ -细化算法

在给定参数化域上求解 IGA 问题，得到初始解 $$u_h$$ 。

计算 $J = \sum_K h_K^2 \int_K (f - \Delta u_h)^2 dK$ 。

计算 $\Delta u_h$ （利用 $u_h = T(\xi, \eta)$ 和 $J = x_\xi y_\eta - x_\eta y_\xi$ 推导链式公式）。

以 $$J$$ 为目标函数，用最速下降法重排内部控制点：搜索方向 $d_k = -\nabla J$ ，沿方向做回溯线搜索找到最优步长。

返回第 1 步，迭代至收敛。

3.3 Winslow 映射方法

[Wang et al., JCAM 2019] 提出了基于 Winslow 映射与监测函数（monitor function）的 $$r$$ -自适应 IGA：

核心思想：利用解的主曲率分布作为监测函数，引导内部控制点向高曲率区域聚集，从而在解变化剧烈的区域获得更高分辨率。

G_1(\xi, \eta) = |k_1| + |k_2|

其中 $$k_1, k_2$$ 是解曲面的两个主曲率。监测函数大的区域对应解的梯度大，因此需要更密的控制点覆盖。

算法流程：

在初始参数化域上用 IGA 求解 PDE。

基于当前解计算监测函数 $$G_1$$ 。

用 Winslow 映射 $x(\xi,\eta)$ 建立新参数化，使新参数化的网格密度与 $$G_1$$ 成正比。

回到第 1 步，迭代直至 $$G_1$$ 收敛。

该方法在精度和效率上均优于固定网格 IGA，尤其当解具有尖锐特征（边界层、奇点）时优势显著。

小结

主题	核心内容
泊松方程强形式	热传导稳态方程 + 三类边界条件（迪利克雷/纽曼/罗宾）
弱形式推导	分部积分 → 变分形式 $$a(w,u)=L(w)$$
Galerkin 离散	NURBS 基函数构成有限维子空间 → 线性系统 $$Kd=F$$
单元矩阵装配	雅可比变换 + 高斯积分
$$r$$ -细化	固定自由度，优化内部控制点位置，最小化残差能量模
Winslow 方法	以解主曲率为监测函数，通过映射引导网格自适分布

IGA 的核心理念在于几何与物理的统一表示：NURBS 既是最精确的几何建模语言，也是最自然的物理场近似基底。本讲以泊松方程为载体，展示了从 PDE 弱形式推导到矩阵系统装配、再到自适应细化的完整管线，为后续更复杂的 NURBS 有限元分析和等几何拓扑优化奠定了基础。

参考论文：

Xu et al., Parameterization of computational domain in IGA: Methods and comparison, CMAME 2011

Xu et al., Efficient r-adaptive IGA with Winslow's mapping, JCAM 2019