线弹性问题求解与 GIFT 方法

2026/05/20 14:00:00·2026/05/21 11:30:00

线弹性问题的数学建模

线弹性问题是弹性力学中最经典的模型，也是理解 IGA 工作原理的最好起点。与有限元方法（FEA）相比，IGA 使用 $C^{p-1}$ 连续的 NURBS 基函数来同时表示几何和位移场，这一特性使得 IGA 在处理连续应力场时具有显著优势。

基本方程与强形式

二维线弹性问题的平衡方程强形式为：在区域 $\Omega$ 内，

\sigma_{AB,B} + f_A = 0

其中 $\sigma_{AB}$ 为 Cauchy 应力张量， $$f_A$$ 为单位体积体力。边界条件分为两类：

Dirichlet 边界 $\Gamma_D^{(A)}$ ：给定位移 $\vec{g}_A$
Neumann 边界 $\Gamma_N^{(A)}$ ：给定边界牵引力 $\bar{h}_A$

应变与位移的关系使用 工程应变（小变形假设下的线性化应变）：

\varepsilon_{AB} = \frac{1}{2}(u_{A,B} + u_{B,A})

本构方程满足胡克定理 $\sigma_{AB} = C_{ABCD} \varepsilon_{CD}$ ，其中 $$C$$ 为弹性材料系数矩阵。对于各向同性体，使用拉梅常数 $\lambda$ 和 $\mu$ （或杨氏模量 $$E$$ 和泊松系数 $\nu$ ）来表征材料性质。

线弹性问题的强形式与 Dirichlet/Neumann 边界条件

从强形式到弱形式（虚功原理）

为得到弱形式，对每个方向引入 trial 函数（位移解的候选函数） $$u_A$$ 和 weighting 函数（检验函数） $$w_A$$ 。定义函数空间：

\Phi = \{ u \mid u_A \in H^1(\Omega), \; u_A|_{\Gamma_D^{(A)}} = \bar{g}_A \}

\mathcal{V} = \{ w \mid w_A \in H^1(\Omega), \; w_A|_{\Gamma_D^{(A)}} = 0 \}

将 trial 函数和 weighting 函数代入平衡方程，乘以 $$w_A$$ 后做分步积分，得到弱形式：

\int_\Omega \varepsilon_{AB}(w) \sigma_{AB}(u) \, d\Omega = \int_\Omega w_A f_A \, d\Omega + \int_{\Gamma_N} w_A \bar{h}_A \, d\Gamma \quad \forall \vec{w} \in \mathcal{V}

为了简化表达，定义双线性形式 $$a(u, w)$$ 和线性泛函 $$L(w)$$ ：

a(u, w) = \int_\Omega \varepsilon_{AB}(w) \sigma_{AB}(u) \, d\Omega, \quad L(w) = \int_\Omega w_A f_A \, d\Omega + \int_{\Gamma_N} w_A \bar{h}_A \, d\Gamma

则弱形式可写为：Find $u \in \Phi$ s.t.

a(u, w) = L(w) \quad \forall w \in \mathcal{V}

在固体力学中，这被称为 虚功原理。 $$w$$ 被称为虚位移， $$L(w)$$ 是体力乘以虚位移加上边界牵引力做的虚功， $$a(u, w)$$ 代表应力做的功。方程描述了内力功和外力功的平衡。

Voigt 标记法

将对称张量写成向量形式可以简化后续的矩阵推导。应变向量定义为：

\vec{\varepsilon}(\vec{u}) = [u_{1,1}, \; u_{2,2}, \; u_{1,2} + u_{2,1}]^{\mathsf{T}}

虚应变向量和应力向量同理。对应的张量分量到向量分量的映射关系为：

向量索引	张量索引 A	张量索引 B
1	1	1
2	2	2
3	1	2

这样本构方程变为向量-矩阵形式 $\vec{\sigma} = \mathbf{D} \vec{\varepsilon}$ ，弹性矩阵 $\mathbf{D}$ 的维度从 $2 \times 2 \times 2 \times 2$ 降低到 $3 \times 3$ ，大幅简化了计算。

平面应变与平面应力

课程详细讨论了各向同性体的两种经典二维设定。

平面应力（薄板，忽略 $$z$$ 方向应力）： $\sigma_{33} = \sigma_{13} = \sigma_{23} = 0$ ，弹性矩阵为

\mathbf{D} = \frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix}

平面应变（长结构如大坝，约束 $$z$$ 方向应变）： $\varepsilon_{33} = \varepsilon_{13} = \varepsilon_{23} = 0$ ，弹性矩阵为

\mathbf{D} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix}

两种设定导致弹性矩阵中的系数差异，理解这一差异对于正确建立模型至关重要。对于平面应力， $$z$$ 方向法向应变 $\varepsilon_{33}$ 不为零且可由应力反推；对于平面应变， $$z$$ 方向法向应力 $\sigma_{33}$ 不为零。

Galerkin 方法：弱形式的离散化

从弱形式到代数方程组

将弱形式转化为可求解的代数方程组，使用 Galerkin 方法在有限维子空间 $\mathcal{V}^h \subset \mathcal{V}$ 和 $\mathcal{U}^h \subset \mathcal{U}$ 中计算，这些集合由等参的 NURBS 基定义。

将位移和虚位移用 NURBS 基函数展开：

u_A^h(\vec{x}) = \sum_{i \in \eta} N_i(\vec{x}) \, d_{i,A}, \quad w_A^h(\vec{x}) = \sum_{j \in \eta} N_j(\vec{x}) \, c_{j,A}

定义 $\eta = \{1, \dots, n\}$ 为包含所有 NURBS 基函数索引的集合， $\eta_D^{(A)} \subset \eta$ 为在 $\Gamma_D^{(A)}$ 上非零的基函数索引集合。将展开式代入虚功方程，利用 $$w$$ 的任意性（任意 $c_{j,A}$ ），得到：

\sum_{j} \sum_{B} a(N_i \hat{e}_A, N_j \hat{e}_B) \, d_{j,B} = L(N_i \hat{e}_A) \quad \text{for } i \in \eta - \eta_D^{(A)}, \; A=1,\dots,d

其中 $\hat{e}_1 = [1, 0]$ , $\hat{e}_2 = [0, 1]$ 。整理后得到最终矩阵形式方程：

\mathbf{K} \vec{d} = \vec{F}

其中 $\mathbf{K}$ 为 全局刚度矩阵， $\vec{d}$ 为节点位移向量， $\vec{F}$ 为外力向量。

ID 数组与自由度映射

对于每个全局基函数有 $$d$$ 个方程（对应 $$d$$ 个自由度），引入 ID 数组将自由度 $$A$$ 和全局基索引 $$i$$ 关联到方程索引 $$P$$ ：

P = \text{ID}(A, i), \quad Q = \text{ID}(B, j)

ID 数组有两种构造策略：

方法一（基于基函数索引）：将未知量按基函数组织，第 $$i$$ 个基函数的 $$d$$ 个自由度排在一起。刚度矩阵维度为 $n \times n$ ，每个元素 $k_{ij}$ 是 $d \times d$ 的子块：

K_{pq} = k_{ij}^{AB} \quad p = d(i-1)+A, \; q = d(j-1)+B

这种方法得到的刚度矩阵稀疏性更好。

方法二（基于自由度索引）：将未知量按物理自由度组织，自由度 $$A$$ 的所有基函数索引排在一起。刚度矩阵维度为 $(d \times n) \times (d \times n)$ ：

K_{PQ} = k_{ij}^{AB} \quad P = n(A-1) + i, \; Q = n(B-1) + j

这种方法使矩阵系统按物理量排序，可高效使用基于物理的线性求解器和预处理器。

Galerkin 方法推导与最终的矩阵形式方程

\mathbf{Kd} = \mathbf{F}

矩阵装配与单元构造

装配矩阵系统

根据虚功原理，构建大量局部刚度矩阵和局部载荷向量后，需要将其装配到全局系统中。在 $$d$$ 维空间和 $n_{loc}$ 个局部基函数的条件下，单元 $\Omega^e$ 的局部刚度矩阵为：

k_{ab}^{e,(AB)} = \int_{\Omega^e} (B_{a,A})^{\mathsf{T}} D \, B_{b,B} \, d\Omega

其中 $$p = d(a-1)+A$$ , $$q = d(b-1)+B$$ 。局部载荷向量为：

f_{p}^{e,(A)} = \int_{\Omega^e} N_a f_A \, d\Omega + \int_{\Gamma_N^e \cap \Gamma_N} N_a \bar{h}_A \, d\Gamma

装配过程使用以下数组：

IEN 数组： $i = \text{IEN}(a, e)$ ，局部基索引 $$a$$ 到全局基索引 $$i$$ 的映射
LM 数组： $P = \text{LM}(A, a, e) = \text{ID}(A, \text{IEN}(a, e))$ ，局部到全局方程索引的映射

Dirichlet 边界条件的处理方式是：不更新刚度矩阵中对应的行和列，而是直接修正载荷向量。若 $$BC(A,i) = 1$$ （表示自由度受约束），设定 $K_{PP} = 1$ , $F_P = (\bar{g}_A)_i$ 。

单元刚度矩阵的数值构造

使用等参思想、Bezier 提取和四边形高斯积分。二维情况下：

k_{ab}^{e,(AB)} = \int_{\hat{\Omega}^e} (\hat{B}_{a,A})^{\mathsf{T}} D(\vec{x}(\hat{\xi})) \, \hat{B}_{b,B} \, |J| \, d\hat{\xi}

数值积分采用二维高斯求积：

k_{ab}^{e,(AB)} \approx \sum_{q_1=1}^{n_g} \sum_{q_2=1}^{n_g} (\hat{B}_{a,A})^{\mathsf{T}} D \, \hat{B}_{b,B} \, |J| \, w_{q_1} \, w_{q_2}

IEN/LM 数组映射关系与全局刚度矩阵装配流程

载荷向量分为两部分构造：

1. 体积力部分（如重力、离心力）：与刚度矩阵构造方式相同，在单元内对基函数和体力做积分

2. 表面力部分（如牵引力）：需要精确获取边界 $\Gamma_3^e \cap \Gamma^e$ ，通过 Neumann 数组标识对应的 Bezier 单元边，沿边界做降维积分

单元刚度矩阵和体积力/表面力载荷向量的构造算法伪代码

悬臂梁算例对比

课程通过悬臂梁经典算例展示了 IGA 与 FEA 的数值差异。在相同单元规模下，IGA 得到的冯·米斯应力分布精度更高；在粗网格下 IGA 仍能保持极高精度，这意味着可以通过减少单元规模来加速仿真过程。

端点处有限元分析、等几何分析和解析解的位移大小对比

冯·米斯应力分布及 x/y 方向位移分布对比。相同单元下等几何比有限元数值精度更高。

GIFT：几何无关的场逼近

问题的提出

IGA 的核心思想是几何和位移场使用同一套 NURBS 样条空间表示，实现了 CAD 与 CAE 的无缝融合。然而在实际工程中，这一设计也带来了一些约束：

给定 NURBS 几何往往无法精确转换为指定的样条空间
几何表示和场逼近共用一套节点向量，降低了灵活性
解场的连续性受几何域固定，不利于处理局部强梯度问题

GIFT（Geometry-Independent Field approximaTion） 方法正是为了解决这一问题而提出。其核心思想是：几何与场使用独立的样条空间，两者通过映射关系建立联系。

IGA 与 GIFT 的核心区别：几何与解场是否共用样条空间

GIFT 的数学框架

GIFT 在几何无关的样条空间中寻求解：

u^h(\vec{x}) = \sum_{i_1,i_2,\dots,i_d} M_{i_1,i_2,\dots,i_d}(\vec{\xi}) \, c_{i_1,i_2,\dots,i_d}

其中 $M_{i_1,i_2,\dots,i_d}(\vec{\xi})$ 可以是 NURBS、B-spline、T-spline、PHT-spline 等任意样条的 tensor product。几何映射 $\vec{x} = \mathbf{F}(\vec{\xi})$ 保持不变，而解空间则可以选择更灵活的样条表示。

IGA 与 GIFT 的对比

特性	IGA	GIFT
解场样条空间	与几何相同	几何无关
计算域	$$h/p$$ 加密时变化	固定不变
积分单元	几何参数域的节点区间	解场参数域的节点区间
自由度	几何域决定	可独立选择
连续性控制	随几何域固定	灵活控制

GIFT 的关键优势在于：

灵活性：对于解场有局部强梯度的问题（如含尖锐峰值的解），可以独立地选择高连续性的场基函数而不受几何限制
连续控制：在多片区域中，相邻片之间解场的连续性可通过选择适当的场样条空间自动保证
独立加密：加密操作直接在解场的样条空间上执行，与几何无关
预计算优化：雅可比矩阵可预计算，减少重复计算开销

分片检验（Patch Test）

分片检验由 Irons 于 1965 年提出，是评估位移型单元收敛性的经典方法。虽然 Stummel 证明了通过分片检验并非收敛的充要条件，但它仍是测试单元收敛性的有效工具。

课程通过圆环域（quarter annulus）的分片检验验证了 GIFT 方法的收敛可靠性。检验几何由 1x2 次 NURBS 精确表示，控制点经过精心设计以生成三种不同形态的参数量化单元（均匀、非均匀、不规则）。

在多种组合下进行检验：

$$p_u = p_g$$ , $\Xi_u = \Xi_g$ （等价于 IGA）
$$p_u = p_g$$ , $\Xi_u \neq \Xi_g$ （同次数不同节点向量）
$$p_u < p_g$$ , $\Xi_u = \Xi_g$
$$p_u < p_g$$ , $\Xi_u \neq \Xi_g$
$$p_u > p_g$$ , $\Xi_u = \Xi_g$
$$p_u > p_g$$ , $\Xi_u \neq \Xi_g$

结果显示：当节点向量相同时（等价于 IGA），各种参数组合均通过检验（误差达机器精度 $10^{-15}$ 量级）；当节点向量不同时，GIFT 在矩形参数域上仍能通过检验，且解精度优于 IGA。

自适应 PHT-GIFT

课程介绍了基于 PHT-spline（分片多项式层次 T 样条） 的自适应 GIFT 方法。PHT-spline 由 Deng 等人于 2008 年提出，具有线性无关性和简单局部加密特性，非常适合作为 GIFT 的解场样条空间。

自适应流程为：

1. 在 NURBS 几何域上使用 GIFT 求得 PHT-spline 解场

2. 在解场上逐片计算局部误差指标

3. 通过 均值标记算法 标记需加密的参数单元

4. 将标记单元细分为 4 个子单元

5. 在细化后的 T-mesh 上构造新的 PHT-spline 基函数

6. 求解新的解场

7. 重复直到估计误差小于阈值

这一方法结合了 NURBS 的精确几何表示能力和 PHT-spline 的局部加密特性。预处理阶段需要构造几何和解场的公共参数域：对于单片 NURBS 或多片拼接域，解场的参数域由 NURBS 参数域的节点线划分构成的网格形成。

算例以含尖锐峰值的圆环域问题验证了 PHT-GIFT 的自适应能力——在解场的峰值区域，T-mesh 会自动加密以捕获局部高梯度特征，而平滑区域保持粗网格，显著提高了计算效率。

小结

第 7 讲从线弹性问题的虚功原理出发，系统推导了 IGA 方法的 Galerkin 离散化框架，建立了从连续方程到矩阵方程的完整链路。关键内容包括：

线弹性问题的强形式与弱形式（虚功原理）
Voigt 标记法与平面应力/应变设定
Galerkin 离散化与 $\mathbf{Kd} = \mathbf{F}$ 的完整推导
IEN/LM 数组的矩阵装配与边界条件处理
悬臂梁算例验证 IGA 相对于 FEA 的精度优势
GIFT 方法通过解耦几何样条空间与场逼近样条空间，突破了 IGA 的几何绑定限制
PHT-GIFT 结合 NURBS 精确几何与 PHT-spline 自适应加密

GIFT 的核心思想可以总结为一句话：几何只告诉你问题在哪里，解场告诉你该如何求解 —— 两者应该独立设计，而非绑定在一起。