超弹性问题求解与等几何配点法

在进入超弹性问题的详细讨论之前，有必要梳理非线性现象的各类来源。课程首先从物理角度对非线性问题做了系统分类：

• 几何非线性：大位移和大旋转，但应变仍可假设为小量（框架、膜、壳体结构）

• 有限形变：位移和应变同时大（金属成型、轮胎力学）

• 材料非线性：本构关系非线性（粘弹聚合物、弹塑性材料）

• 稳定性问题：平衡方程的不稳定性（壳体屈曲）

• 耦合问题：多物理场交互（流固耦合等）

超弹性（Hyperelasticity） 属于材料非线性中最简单的一类，其核心特征是应力-应变关系由 应变能密度函数 定义。本构方程不再满足胡克定律的线性关系，而是通过能量密度函数对应变张量的导数来计算应力。

---

1 本构方程

应变能密度函数 $\psi(\mathbf{F})$ 是超弹性材料的核心，其中 $\mathbf{F}$ 为 形变梯度，定义为物体当前坐标相对于初始坐标的偏导数 $\mathbf{F} = \partial \vec{x} / \partial \vec{X}$。由能量共轭关系，应力通过应变能密度函数对格林应变的导数给出：

$\mathbf{S} = \frac{\partial \psi(\mathbf{G})}{\partial \mathbf{G}} = 2 \frac{\partial \psi(\mathbf{C})}{\partial \mathbf{C}}$

其中 $\mathbf{G}$ 为 格林应变（Green Strain），$\mathbf{C} = \mathbf{F}^{\mathsf{T}} \mathbf{F}$ 为 Cauchy 应变张量。

常见的超弹性材料模型包括：

St. Venant-Kirchhoff 模型（格林应变的简单多项式形式）：

$\psi = \frac{\lambda}{2}(\text{tr}\,\mathbf{G})^2 + \mu \, \text{tr}(\mathbf{G}^2)$

Neo-Hookean 模型（不可压缩橡胶材料的首选模型）：

$\psi = \frac{\mu}{2}(\text{tr}(\mathbf{C}) - 3) + \frac{\lambda}{2}(\det\mathbf{F} - 1)^2$

Mooney-Rivlin 模型 和 Ogden 模型 则是针对橡胶类材料的更精细模型，后者使用主拉伸 $\lambda_1, \lambda_2, \lambda_3$ 作为自变量，形式更为灵活但参数也更多。

2 几何方程与格林应变

在线弹性中使用的工程应变 $\varepsilon = \frac{1}{2}(\mathbf{F}^{\mathsf{T}}\mathbf{F} - \mathbf{I})$ 是格林应变的线性部分 $\mathbf{G}_L$，忽略了非线性部分 $\mathbf{G}_N$。当处理大变形时，这两部分都必须保留：

$\mathbf{G} = \mathbf{G}_L + \mathbf{G}_N$

这使得几何方程从线性变为非线性，是超弹性问题区别于线弹性问题的根本所在。

3 平衡方程与离散化

在拉格朗日构型下，利用线动量守恒建立动力学平衡方程：

$\mathbf{M} \ddot{\vec{u}} + \mathbf{C} dot{\vec{u}} + \mathbf{R}(\vec{u}) = \vec{F}$

其中 $\mathbf{M}$ 为质量矩阵，$\mathbf{C}$ 为阻尼矩阵，$\mathbf{R}$ 为内力向量（由第二类皮奥拉-克希荷夫应力 $\mathbf{S}$ 的积分给出），$\vec{F}$ 为外力向量。

在等几何分析框架下进行离散化：

• 位移：$u(\vec{\xi}) = \sum_{i} N_i(\vec{\xi}) \vec{u}_i$

• 形变梯度：$\mathbf{F} = \mathbf{I} + \sum_i u_i \frac{\partial N_i}{\partial \vec{\xi}} \frac{\partial \vec{\xi}}{\partial \vec{x}}$

• 格林应变：$\mathbf{G} = \mathbf{B}_L \vec{u}_e + \mathbf{B}_N(\vec{u}_e)$

离散化后的平衡方程记为 $\mathbf{r}(\vec{u}) = \mathbf{K}\vec{u} - \vec{F} = 0$，其中 $\mathbf{K}$ 为刚度矩阵。注意这里 $\mathbf{B}_L$ 和 $\mathbf{B}_N$ 的形式不同，导致方程非线性且刚度矩阵 非对称。

---

超弹性问题的离散平衡方程是非线性的，必须通过迭代求解。课程采用 Newton-Raphson 迭代法，其核心是对残差函数做泰勒展开并取前两项：

$\vec{r}(\vec{u} + \Delta\vec{u}) \approx \vec{r}(\vec{u}) + \mathbf{J}(\vec{u}) \Delta\vec{u} = 0$

其中 $\mathbf{J} = \partial \vec{r} / \partial \vec{u}$ 为 切线刚度矩阵。

切线刚度矩阵分为两部分：

• $\mathbf{K}_L = \mathbf{B}_L^{\mathsf{T}} \mathbf{D} \mathbf{B}_L$：小位移刚度矩阵，即线弹性问题中的刚度矩阵

• $\mathbf{K}_N = \mathbf{B}_L^{\mathsf{T}} \mathbf{D} \mathbf{B}_N + \mathbf{B}_N^{\mathsf{T}} \mathbf{D} \mathbf{B}_L + \mathbf{B}_N^{\mathsf{T}} \mathbf{D} \mathbf{B}_N$：大变形关联刚度矩阵，由几何非线性贡献

牛顿迭代的步骤为：设定初始位移 $\vec{u}^0$ → 组装 $\mathbf{K}^T(\vec{u}^n)$ → 求解 $\Delta\vec{u}^n$ → 更新 $\vec{u}^{n+1} = \vec{u}^n + \Delta\vec{u}^n$ → 重复直到收敛。

初始值通常设为零，但若设得过远可能导致不收敛或收敛缓慢，需根据具体问题调整。

1 静力学算例：厚壁圆柱

课程以厚壁圆柱为例，对比线弹性与超弹性模型在大变形下的精度。使用三次细分后的超弹性仿真结果作为精确解 $u^*$，计算不同多项式次数和单元密度下的 L2 和 H1 范数相对误差。

悬臂梁大变形算例结果

结果显示：相同单元步长下超弹性模型的精度高于线弹性模型，且随着单元密度增加，超弹性模型误差收敛得更快。这说明对于大变形的精确模拟，使用超弹性本构模型是必要的。

---

1 运动方程

处理模型运动变形过程时，需考虑状态变量随时间的变化。运动方程的一般形式为：

$\mathbf{M} \ddot{\vec{u}} + \mathbf{C} \dot{\vec{u}} + \mathbf{R}(\vec{u}) = \vec{F}(t)$

2 显式与隐式方法的对比

显式时间积分（如中心差分法）仅需上一时刻的物理量，方法简单清晰，在质量矩阵具有对角结构时计算效率高。但其收敛性受步长限制，通常需要极小的步长。

隐式时间积分（如 Newmark 方法）依赖前一时刻和未知时刻的物理量，每个时间步需求解非线性方程组，通常与牛顿迭代法结合使用。其无条件稳定的优点使得步长可以选得更大，适合超弹性仿真对稳定性的高要求。

3 Newmark 方法

Newmark 方法中，下一时刻的位移和速度通过以下公式计算：

$\dot{\vec{u}}_{n+1} = \dot{\vec{u}}_n + \Delta t[(1-\gamma)\ddot{\vec{u}}_n + \gamma \ddot{\vec{u}}_{n+1}]$

$\vec{u}_{n+1} = \vec{u}_n + \Delta t \dot{\vec{u}}_n + \Delta t^2 \left[\left(\frac{1}{2} - \beta\right)\ddot{\vec{u}}_n + \beta \ddot{\vec{u}}_{n+1}\right]$

参数 $\gamma \geq 0.5$、$\beta \geq (\gamma + 0.5)^2/4$。

将位移表达式代入运动方程，得到关于下一时刻位移的非线性代数方程组。在每个时间步内，通过牛顿迭代求解该方程组。最终的迭代格式为：

$\Delta\vec{u}_{n+1}^{(i+1)} = \left[\alpha_1 \mathbf{M} + \alpha_2 \mathbf{C} + \mathbf{K}^T(\vec{u}_{n+1}^{(i)})\right]^{-1} \left[\vec{F}_{n+1} - \mathbf{R}(\vec{u}_{n+1}^{(i)})\right]$

其中 $\alpha_1, \alpha_2$ 为 Newmark 参数变换后的常数项。整个求解过程形成外层（时间步）和内层（牛顿迭代）的双层迭代结构。

4 仿真对比

课程通过数值实验对比了 FEA 和 IGA 在处理大变形问题时的表现：有限元方法在大变形下容易出现网格扭曲变形，需频繁重新划分网格；而等几何方法利用 NURBS 基函数的高阶连续性，即使在大变形下仍能保持光顺的几何和应力表示，无需重新划分网格，在精度和效率上均优于 FEA。

---

1 为什么需要配点法

Galerkin 方法使用方程的弱形式，需要在计算区域内进行 数值积分（高斯积分等），这在几何形状复杂时计算代价高昂。等几何配点法（Isogeometric Collocation） 直接在配点集上求解偏微分方程的强形式，避免了区域内部的数值积分，因此能显著降低计算开销。

配点法与 IGA 的契合之处在于：强形式方程中出现二阶导数，而 IGA 选用的 样条基函数具有高阶连续性 $C^{p-1}$，天然满足强形式对导数的需求。

2 配点的选择策略

配点的选取是配点法的关键步骤。课程介绍了两种主要方法：

Greville 坐标方法：给定节点向量 $\{\xi_0, \dots, \xi_{m+p}\}$，Greville 坐标为 $\bar{\xi}_i = \frac{1}{p}\sum_{j=1}^{i+p} \xi_j$，$i=1,\dots,n=m-p$。该方法被广泛使用，计算简便。

Demko 坐标方法：Chebyshev 样条的极值点坐标，通过迭代算法获得。

二维参数域内的配点为 $\tau_{ij} = (\bar{\xi}_i, \bar{\eta}_j)$。

3 配点处的装配

将配点代入偏微分方程的强形式，得到线性方程组 $\mathbf{K} \vec{U} = \vec{F}$。由于计算域可能由多个样条面片组成，配点分为三类：面片内部配点、面片间公共配点、边界配点，各类配点分别对应不同的装配公式。

4 计算效率对比

课程给出了 Galerkin 方法与配点法的计算时间对比实验（泊松方程，精确解 $x^2+y^2$）：

随着自由度的增大，Galerkin 方法与配点法的时间之比呈增大趋势，证明在大规模问题中配点法能带来非常可观的效率提升（接近一个数量级）。

5 应用：浮雕模型生成

课程还展示了配点法在书法浮雕生成中的应用。通过二值化图像、参数化、静态/动态浮雕面生成等流程，利用配点法高效求解泊松方程，实现书法的三维浅浮雕表达。

---

第8讲从超弹性本构模型出发，构建了从应变能密度函数到离散平衡方程的完整推导链。面对非线性带来的挑战，牛顿迭代法与 Newmark 隐式时间积分构成了动静力学的核心求解框架。等几何配点法则从另一条路径出发——通过直接在强形式上配点求解，避开了数值积分的计算瓶颈，在大规模问题中展现出接近一个数量级的效率优势。两部分内容共同构成了IGA处理非线性问题的方法论基础。