连续介质力学基础
变形梯度、应变与应力张量、非线性FEM公式
连续介质力学是所有结构和流体CAE仿真背后严格的数学框架。线性结构力学课程往往跳过大变形的微妙之处——但一旦仿真橡胶密封圈、金属成形过程、碰撞事故或生物组织，就会遇到多种应力度量、共轭应变对，以及拉格朗日与欧拉描述的选择问题。掌握这些概念，才能真正区分能够正确建立非线性求解器的工程师，与那些得到错误答案却浑然不知的使用者。
1. 运动学：运动描述与变形梯度
教授，Abaqus理论手册里到处都提到变形梯度 $\mathbf{F}$，它到底是什么？为什么比直接用"应变"更基本？
$\mathbf{F}$ 是对局部变形最完整的描述——它在一个张量里同时包含了拉伸、旋转和剪切。简单理解：它回答了"一个无限小的材料纤维在物体变形后如何变化？"这个问题。如果一段材料纤维在参考（未变形）构型中方向为 $d\mathbf{X}$，变形后变为 $d\mathbf{x} = \mathbf{F}\,d\mathbf{X}$。就这一个方程，包含了局部变形的全部信息。Green-Lagrange等应变度量是从 $\mathbf{F}$ 中剥离旋转部分后推导出来的——因为纯刚体旋转在材料中不产生应力，只有拉伸才会。
1.1 参考构型与当前构型
设 $\mathbf{X}$ 为材料点在 $t=0$ 时参考（拉格朗日）构型中的位置，$\mathbf{x} = \boldsymbol{\chi}(\mathbf{X},t)$ 为其在当前（欧拉）构型中的位置。位移向量为：
1.2 变形梯度
$\mathbf{F}$ 的关键性质：
- $J = \det(\mathbf{F}) > 0$——雅可比行列式等于体积比 $dv/dV$
- 不可压缩材料（橡胶、软组织、塑性变形）满足 $J = 1$
- 一般情况下 $\mathbf{F}$ 不对称
- 参考构型下：$\mathbf{F} = \mathbf{I}$
1.3 极分解
$\mathbf{F}$ 可唯一分解为：
其中 $\mathbf{R}$ 为正交旋转张量（$\mathbf{R}^T\mathbf{R} = \mathbf{I}$，$\det\mathbf{R} = 1$），$\mathbf{U}$ 为右拉伸张量（对称正定），$\mathbf{V}$ 为左拉伸张量。$\mathbf{U}$ 的特征值即为主拉伸比 $\lambda_1, \lambda_2, \lambda_3$。
$\mathbf{F} = \mathbf{R}\mathbf{U}$ 这个分解对求解器有什么实际用处？
极分解用于定义客观（坐标系无关）的应力率，确保本构模型不会因刚体旋转而虚假产生应力。在协转坐标系公式中——LS-DYNA碰撞仿真广泛使用——应力在旋转坐标系 $\mathbf{R}^T\boldsymbol{\sigma}\mathbf{R}$ 中计算，旋转已被剥离出去。否则，对于经历大刚体转动的材料单元，本构律的时间积分会累积出虚假剪切应力。Jaumann率、Green-Naghdi率和Truesdell率都以不同方式涉及 $\mathbf{R}$。
2. 应变度量
2.1 Green-Lagrange应变（物质坐标系）
Green-Lagrange应变张量定义在参考构型中，刚体运动时为零：
其中 $\mathbf{C} = \mathbf{F}^T\mathbf{F}$ 为右Cauchy-Green变形张量。用位移梯度 $H_{iJ} = \partial u_i/\partial X_J$ 展开：
非线性项 $H_{kI}H_{kJ}$ 是Green-Lagrange应变与工程小应变张量的区别所在。小变形时，$E_{IJ} \approx \varepsilon_{IJ} = \frac{1}{2}(H_{IJ} + H_{JI})$。
2.2 Almansi应变（空间坐标系）
Almansi（欧拉）应变张量定义在当前构型中：
其中 $\mathbf{B} = \mathbf{F}\mathbf{F}^T$ 为左Cauchy-Green（Finger）张量。Almansi应变用于空间描述的公式，与Cauchy应力在空间描述中形成共轭对。
2.3 对数（Hencky）应变
对数应变是大弹塑性变形最自然的度量：
一维情况：$\varepsilon_H = \ln(L/L_0)$。Hencky应变对连续变形具有可加性，是金属塑性的首选度量（Abaqus、Ansys、LS-DYNA均以"对数应变"输出）。
| 应变度量 | 参考系 | 公式 | 适用场景 |
|---|---|---|---|
| 工程应变 $\varepsilon$ | 参考构型 | $(L-L_0)/L_0$ | 小变形 |
| Green-Lagrange $\mathbf{E}$ | 参考构型（物质） | $\frac{1}{2}(\mathbf{C}-\mathbf{I})$ | 全拉格朗日、超弹性 |
| Almansi $\mathbf{e}$ | 当前构型（空间） | $\frac{1}{2}(\mathbf{I}-\mathbf{B}^{-1})$ | 更新拉格朗日 |
| Hencky $\boldsymbol{\varepsilon}_H$ | 当前构型 | $\ln\mathbf{V}$ | 金属塑性、橡胶大变形 |
3. 应力度量：Cauchy、PK1、PK2
Abaqus输出"Cauchy应力"，理论手册的本构律部分又提到"PK2应力"。我们为什么需要不止一种应力？
每种应力度量描述的是同一个物理事实，只是参考坐标系不同，各有其数学便利性。Cauchy应力是你实际感受到的——当前面积上的力。用压力表测到的就是这个。但对于写在参考构型中的本构律（比如橡胶的超弹性模型），你需要将应力参考到原始未变形面积，这就是PK2应力存在的原因——它是对称的，活在参考坐标系里，是Green-Lagrange应变的能量共轭量，写本构律最方便。PK1是中间过渡——当前构型的力除以参考构型面积——它不对称，直接用于本构律不方便。
3.1 Cauchy应力（真应力，$\boldsymbol{\sigma}$）
Cauchy应力张量定义在当前构型中。当前构型中法向量为 $\mathbf{n}$ 的面元上的面力：
由力矩平衡，$\boldsymbol{\sigma}$ 对称：$\sigma_{ij} = \sigma_{ji}$。这是Abaqus输出的S11、S22、S33、S12、S13、S23。
3.2 第一Piola-Kirchhoff应力（PK1，$\mathbf{P}$）
PK1将参考构型面积法向量映射到当前构型力：
PK1不对称：$P_{iJ} \neq P_{Ji}$。在全拉格朗日公式的弱形式中自然出现，但很少直接用于本构律。
3.3 第二Piola-Kirchhoff应力（PK2，$\mathbf{S}$）
PK2完全在物质坐标系中，对称，且是Green-Lagrange应变的能量共轭量：
单位参考体积的内虚功：
对于超弹性材料（橡胶、生物组织），本构律最简洁地表示为：
其中 $W(\mathbf{C})$ 为应变能密度函数（例如Neo-Hookean：$W = \frac{\mu}{2}(I_1 - 3) + \frac{\kappa}{2}(J-1)^2$）。
3.4 应力转换关系
4. 本构框架
4.1 超弹性
超弹性材料储存应变能并完全恢复，应变能密度 $W$ 仅取决于变形状态，而非加载路径。常见模型：
| 模型 | $W(I_1, I_2, J)$ | 适用场景 |
|---|---|---|
| Neo-Hookean | $\frac{\mu}{2}(I_1-3) + \frac{\kappa}{2}(J-1)^2$ | 中等应变、简单橡胶 |
| Mooney-Rivlin | $C_{10}(I_1-3) + C_{01}(I_2-3)$ | 橡胶（应变至约150%） |
| Ogden | $\sum_p \frac{\mu_p}{\alpha_p}(\lambda_1^{\alpha_p}+\lambda_2^{\alpha_p}+\lambda_3^{\alpha_p}-3)$ | 大应变、生物组织 |
| Arruda-Boyce | 8链网络模型 | 含锁死效应的聚合物网络 |
4.2 有限变形弹塑性
对于经历大塑性变形的金属（板料成形、碰撞），采用乘法分解：
弹性部分 $\mathbf{F}^e$ 储存弹性能；塑性部分 $\mathbf{F}^p$ 不可逆且体积守恒（金属满足 $\det\mathbf{F}^p = 1$）。Mandel应力 $\boldsymbol{\Sigma} = \mathbf{C}^e\mathbf{S}^e$ 在中间构型中驱动塑性流动。
4.3 框架客观性
本构律必须是坐标系无关的（客观的）：材料响应不能依赖于观察者的刚体运动。当前构型经历旋转 $\mathbf{Q}$ 时，Cauchy应力必须满足：
这个要求排除了朴素的速率方程 $\dot{\boldsymbol{\sigma}} = \mathbb{C}:\mathbf{d}$（不满足客观性），需要使用客观应力率，如Jaumann率：$\overset{\nabla}{\boldsymbol{\sigma}} = \dot{\boldsymbol{\sigma}} - \mathbf{W}\boldsymbol{\sigma} + \boldsymbol{\sigma}\mathbf{W}$。
5. 全拉格朗日与更新拉格朗日公式
Abaqus文档里有"全拉格朗日"和"更新拉格朗日"。做橡胶密封圈压缩仿真时，这两种方法在实践中有什么不同？
两种方法对同一问题给出完全相同的结果——正确实现时它们在数学上等价。区别纯粹在于每步使用的参考构型。全拉格朗日始终将量参考回 $t=0$ 时的原始未变形构型。更新拉格朗日在每个增量步开始时将参考更新为当前变形构型。对于橡胶密封圈，Abaqus内部使用更新拉格朗日以提高效率——刚度矩阵和残差在当前构型中计算，避免了在 $\mathbf{F}$ 中积累大的位移梯度。但本构律通常以物质坐标（全拉格朗日）框架写出，使用PK2/Green-Lagrange共轭对，求解器自动处理坐标转换。
5.1 全拉格朗日（TL）公式
所有积分对参考体积 $V_0$ 进行。虚功原理：
优点：本构律始终使用相同参考；Green-Lagrange应变计算直接；适用于超弹性材料。缺点：超大变形时，从参考到当前的映射计算代价高。
5.2 更新拉格朗日（UL）公式
参考构型在每个增量步 $n$ 开始时更新为当前构型 $V_n$。虚功原理变为：
其中 $\boldsymbol{\tau} = J\boldsymbol{\sigma}$ 为Kirchhoff应力，$\mathbf{d}$ 为变形速率张量。LS-DYNA、Abaqus隐式和显式的大多数实体单元均采用此公式。
| 方面 | 全拉格朗日 | 更新拉格朗日 |
|---|---|---|
| 参考构型 | 固定：$t=0$ 状态 | 每步更新：$t_n$ 状态 |
| 应力度量 | PK2（$\mathbf{S}$） | Cauchy（$\boldsymbol{\sigma}$）或Kirchhoff（$\boldsymbol{\tau}$） |
| 应变度量 | Green-Lagrange（$\mathbf{E}$） | Almansi（$\mathbf{e}$）或变形速率（$\mathbf{d}$） |
| 适用场景 | 橡胶、超弹性、大弹性应变 | 金属塑性、碰撞、大塑性应变 |
6. 非线性FEM实现
6.1 非线性的来源
结构FEM中存在三类不同的非线性：
- 材料非线性： 塑性屈服、超弹性、粘弹性、损伤
- 几何非线性： 大位移/旋转、预应力（薄膜/缆索刚化）、屈曲
- 接触非线性： 边界条件变化、间隙单元、摩擦
6.2 Newton-Raphson迭代
非线性平衡残差 $\mathbf{R}(\mathbf{u}) = \mathbf{f}_{\text{外}} - \mathbf{f}_{\text{内}}(\mathbf{u}) = \mathbf{0}$ 迭代求解。第 $k$ 次迭代：
更新拉格朗日公式的切线刚度矩阵包含三个贡献：
其中 $\mathbf{K}_{\text{几何}}$ 为几何（应力）刚度矩阵——负值时对应屈曲，正值时对应预应力刚化。
6.3 收敛准则
标准收敛检验（Abaqus默认值）：
- 力残差：$\|\mathbf{R}^{(k)}\| / \|\mathbf{f}_{\text{参考}}\| < 5 \times 10^{-3}$
- 位移修正：$\|\Delta\mathbf{u}^{(k)}\| / \|\mathbf{u}^{(k)}\| < 10^{-2}$
- 能量准则：$\Delta\mathbf{u}^{(k)} \cdot \mathbf{R}^{(k)} / \Delta\mathbf{u}^{(0)} \cdot \mathbf{R}^{(0)} < 10^{-5}$
我的Abaqus仿真在模拟橡胶垫片受压时总是不收敛，应该先检查什么？
橡胶受压不收敛，按顺序检查这几点。第一，超弹性模型是否由真实试验数据标定？Neo-Hookean模型如果只用单轴拉伸数据拟合，在压缩下经常失稳，因为双轴/体积行为没有被约束。第二，检查网格——橡胶被压扁时单元畸变严重，雅可比 $J \to 0$ 会导致问题。尝试使用杂交单元（Abaqus中的C3D8H），处理近不可压缩性效果更好。第三，减小初始增量步长，开启自动步长控制。第四，检查接触定义——如果密封圈同时受压和滑移，粘着接触与无摩擦接触对收敛行为影响很大。