1. 商业有限元真正每一步在做什么
超弹性、return mapping、几何刚度、Newton-Raphson、line search、arc-length......实际上，它们分别回答不同的问题。
假设一个大变形隐式分析的第 \(n\) 个增量已经收敛：
\[\mathbf x_n,\qquad \boldsymbol\sigma_n,\qquad \mathrm{state}_n.
\]
现在要求 \(t_{n+1}\)。一次 Newton 迭代可以浓缩为
\[\boxed{\mathbf u\rightarrow\mathrm{kinematics}\rightarrow\mathrm{constitutive\ update}\rightarrow\mathbf f_{\mathrm{int}}\rightarrow\mathbf K\rightarrow\Delta\mathbf u}.
\]
反复迭代，直到
\[\mathbf R=\mathbf f_{\mathrm{ext}}-\mathbf f_{\mathrm{int}}\approx\mathbf0.
\]
真正复杂的地方不是 Newton 公式本身，而是每一个箭头背后都有一套理论。
2. Nonlinear FEM 的八个理论层次
可以按照下面的顺序理解整个 nonlinear FEM：
\[\begin{gathered}
\boxed{\text{1. Configuration / description}}\\
\text{TL、UL、ALE、Eulerian}
\end{gathered}
\]
↓
\[\begin{gathered}
\boxed{\text{2. Kinematics}}\\
\mathbf F,\mathbf C,\mathbf E,\mathbf D,\ \text{有限转动}
\end{gathered}
\]
↓
\[\begin{gathered}
\boxed{\text{3. Stress/strain measures + work conjugacy}}\\
\boldsymbol\sigma,\boldsymbol\tau,\mathbf P,\mathbf S
\end{gathered}
\]
↓
\[\begin{gathered}
\boxed{\text{4. Constitutive model}}\\
\text{hyperelastic / hypoelastic / plasticity / viscoelasticity}
\end{gathered}
\]
↓
\[\begin{gathered}
\boxed{\text{5. Objectivity + constitutive integration}}\\
\text{objective rate、corotational frame、return mapping、history update}
\end{gathered}
\]
↓
\[\begin{gathered}
\boxed{\text{6. Weak form + consistent linearization}}\\
\mathbf f_{\mathrm{int}},\quad
\mathbf K_{\mathrm{mat}},\quad
\mathbf K_{\mathrm{geo}}
\end{gathered}
\]
↓
\[\begin{gathered}
\boxed{\text{7. Element discretization}}\\
\text{shape function、Gauss integration、locking、hourglass、mixed formulation}
\end{gathered}
\]
↓
\[\begin{gathered}
\boxed{\text{8. Nonlinear solution}}\\
\text{Newton、line search、arc-length、explicit dynamics}
\end{gathered}
\]
有了这个框架，商业有限元软件中大量看似零散的选项，都可以重新放回它们所属的理论层。比如，TL/UL 首先属于运动和平衡方程的描述方式；Full Newton、line search、arc-length 则属于全局非线性方程的求解策略。两者不是同一层次的概念。
2.1 第一件事：到底在哪个构形上写平衡方程？
有限变形首先区分
\[\mathbf X\longrightarrow\mathbf x_n\longrightarrow\mathbf x_{n+1}.
\]
其中 \(\mathbf X\) 是初始构形，\(\mathbf x_n\)
是已经收敛的当前构形，\(\mathbf x_{n+1}\) 是待求构形。
2.1.1 Total Lagrangian
TL 始终以初始构形 \(\Omega_0\) 为参考：
\[\mathbf F=\frac{\partial\mathbf x}{\partial\mathbf X},\qquad \mathbf C=\mathbf F^T\mathbf F,
\]
\[\mathbf E=\frac12(\mathbf C-\mathbf I).
\]
与 \(\mathbf E\) 功共轭的第二类 Piola-Kirchhoff 应力为 \(\mathbf S\)：
\[\delta W_{\mathrm{int}}=\int_{\Omega_0}\mathbf S:\delta\mathbf E\,dV_0.
\]
所以经典 TL 中经常看到
\[\boxed{\Omega_0+\mathbf E\leftrightarrow\mathbf S}.
\]
但"选择某个构形就唯一决定某种应力和应变"并不严格。不同度量可以通过 push-forward / pull-back 转换。真正需要保证的是 kinematics、stress measure、virtual work 和 tangent 在数学上保持一致。
2.1.2 Updated Lagrangian
UL 以最近已经收敛的构形作为下一增量的参考：
\[\mathbf x_n\rightarrow\mathbf x_{n+1}.
\]
空间描述中常使用
\[\mathbf L=\nabla_{\mathbf x}\mathbf v=\mathbf D+\mathbf W,
\]
其中
\[\mathbf D=\frac12(\mathbf L+\mathbf L^T),\qquad \mathbf W=\frac12(\mathbf L-\mathbf L^T).
\]
对于 rate-form 本构关系，经常出现
\[\overset{\circ}{\boldsymbol\sigma}=\mathbb C:\mathbf D.
\]
这里就会引出 objective stress rate。但需要特别强调：UL 并不等价于必须使用 Jaumann rate。 构形描述与本构模型如何保 objectivity 是两个相关但不同的问题。
2.1.3 Eulerian
Eulerian 描述固定空间网格，材料可以穿过网格：
\[\frac{D(\cdot)}{Dt}=\frac{\partial(\cdot)}{\partial t}+\mathbf v\cdot\nabla(\cdot).
\]
因此除了材料响应，还必须处理 advection / transport。
ALE 中，网格也可以运动，但是材料怎么运动和计算网格怎么运动可以分开指定。
2.2 Kinematics
Newton solver 得到的是节点自由度，但 constitutive model 真正需要的是 Gauss point 处的变形信息。
对于实体单元，
\[\mathbf x(\boldsymbol\xi)=\sum_aN_a(\boldsymbol\xi)\mathbf x_a.
\]
由此得到 \(\mathbf F\)，再计算 \(\mathbf C\)、\(\mathbf b\)、\(\mathbf E\)，或者率形式中的 \(\mathbf L\)、\(\mathbf D\)、\(\mathbf W\)。
所以程序中的真实链条是
\[\boxed{\mathrm{nodal\ DOFs}\rightarrow\mathrm{element\ interpolation}\rightarrow\mathrm{Gauss\ point\ kinematics}}.
\]
2.3 Stress measure 与 work conjugacy
有限变形下不存在唯一的"应力"。常见量包括 Cauchy stress \(\boldsymbol\sigma\)、Kirchhoff stress \(\boldsymbol\tau\)、第一类Piola-Kirchhoff 应力 \(\mathbf P\) 和第二类 Piola-Kirchhoff 应力 \(\mathbf S\)。他们之前可以相互转换：
\[J=\det\mathbf F,\qquad \boldsymbol\tau=J\boldsymbol\sigma,
\]
\[\mathbf P=\mathbf F\mathbf S,\qquad \boldsymbol\sigma=\frac1J\mathbf F\mathbf S\mathbf F^T.
\]
关键是 work conjugacy。例如 \(\mathbf S:\dot{\mathbf E}\) 和 \(\boldsymbol\tau:\mathbf D\) 都是典型的功共轭组合。有限变形代码中非常危险的错误，就是把不同构形上的 stress、strain 和 tangent 随意混合。
2.4 商业软件的材料点算法
第 \(n\) 个增量已经保存
\[\mathrm{state}_n=\{\boldsymbol\sigma_n,\boldsymbol\varepsilon_n^p,\alpha_n,\ldots\}.
\]
Newton 第 \(i\) 次迭代给材料点新的运动学输入。材料积分器完成
\[\boxed{\mathrm{kinematic\ input}+\mathrm{state}_n\rightarrow\mathrm{stress}_{n+1}^{(i)}+\mathrm{state}_{n+1}^{\mathrm{trial}}+\mathbb C_{\mathrm{alg}}^{(i)}}.
\]
三个关键输出是 trial stress、trial history variables 和 algorithmic consistent tangent。
2.5 Objectivity：为什么刚体转动不能产生假材料响应？
如果一个物体只发生刚体转动而没有真实 deformation，材料不应该凭空产生额外的材料响应，这个特性叫做 objectivity。
对于 rate-form constitutive model，普通的 \(\dot{\boldsymbol\sigma}\) 通常不是 objective 的，因此需要适当的 objective stress rate，或者在 corotational frame 中更新应力。
2.6 虚功方程和它的导数
全局离散平衡为
\[\mathbf R(\mathbf u)=\mathbf f_{\mathrm{ext}}(\mathbf u)-\mathbf f_{\mathrm{int}}(\mathbf u)=\mathbf0.
\]
内部力在一种常见表示下可以概念性写成
\[\mathbf f_{\mathrm{int}}=\int\mathbf B^T\boldsymbol\sigma\,dv.
\]
Newton 需要线性化：
\[\mathbf R(\mathbf u+\Delta\mathbf u)\approx\mathbf R(\mathbf u)+\frac{\partial\mathbf R}{\partial\mathbf u}\Delta\mathbf u.
\]
若采用上面的 residual 定义，
\[\mathbf K=-\frac{\partial\mathbf R}{\partial\mathbf u}.
\]
大变形下常理解为
\[\boxed{\mathbf K=\mathbf K_{\mathrm{mat}}+\mathbf K_{\mathrm{geo}}}.
\]
\(\mathbf K_{\mathrm{mat}}\) 来自材料响应，\(\mathbf K_{\mathrm{geo}}\) 来自当前应力状态和 nonlinear kinematics。因此即使材料完全线弹性，只要考虑大变形，有限元问题仍然可以是非线性的。对于 follower load 等随构形变化的外载，还可能出现额外的 load stiffness。
对于增量型材料，Newton 真正需要的是材料切线刚度：
\[\mathbb C_{\mathrm{alg}}=\frac{\partial\boldsymbol\sigma_{n+1}}{\partial\boldsymbol\varepsilon_{n+1}}.
\]
"材料模型能给出正确 stress"和"材料模型能给出与离散更新一致的 tangent"是两个不同的问题。后者直接影响 global Newton 的收敛性质。
2.7 Element discretization：理论最终必须落到单元上
连续体理论最终必须通过单元离散进入程序，于是出现 shape function、isoparametric mapping、Gauss integration、reduced/full
integration、volumetric locking、shear locking、hourglass、mixed \(u-p\) formulation、enhanced strain 等问题。例如近不可压材料中，即使 constitutive model 完全正确，一个普通低阶 displacement element 仍可能因为 volumetric locking 得到很差的结果。
2.8 全局非线性求解
2.8.1 Full Newton
\[\mathbf K_T(\mathbf u_i)\Delta\mathbf u_i=\mathbf R_i,
\]
\[\mathbf u_{i+1}=\mathbf u_i+\Delta\mathbf u_i.
\]
每次 iteration 重新形成当前 tangent。在 tangent一致并且解足够接近时，可以获得很快的局部收敛。
2.8.2 Modified Newton
一个 increment 内冻结某次 tangent：
\[\mathbf K_0\Delta\mathbf u_i=\mathbf R_i.
\]
残差继续更新，但 tangent 不再每次重新形成。这样可以减少 assembly/factorization 成本，但通常牺牲收敛速度。
2.8.3 Quasi-Newton
Quasi-Newton 不重新计算精确 tangent，而构造
\[\widetilde{\mathbf K}_i\approx\mathbf K_T(\mathbf u_i).
\]
利用
\[\mathbf s_i=\mathbf u_{i+1}-\mathbf u_i,\qquad \mathbf y_i=\mathbf R_{i+1}-\mathbf R_i
\]
更新 Jacobian 近似，使其满足类似 secant condition：
\[\widetilde{\mathbf K}_{i+1}\mathbf s_i\approx\mathbf y_i.
\]
因此三者可以概括为：
- Full Newton：重新计算真实 tangent；
- Modified Newton：冻结旧 tangent；
- Quasi-Newton：不重算真实 tangent，但根据迭代信息更新其近似。
2.8.4 Line search
标准 Newton 使用完整 correction，而 line search 引入
\[\mathbf u_{i+1}=\mathbf u_i+\alpha_i\Delta\mathbf u_i.
\]
它不改变 Newton direction，而是控制沿这个方向走多远。
2.8.5 Automatic increment 与 cutback
成熟商业软件不会在一次 Newton 失败后立即结束，而是执行
\[\boxed{\mathrm{failed\ increment}\rightarrow\mathrm{rollback}\rightarrow\mathrm{cutback}\rightarrow\mathrm{retry}}.
\]
如果连续多个 increment 很容易收敛，还可以自动增大 increment size。
2.8.6 Arc-length / Riks
普通 load control 把 load factor 当作已知量。Arc-length 则同时求 \(\mathbf u\) 和 \(\lambda\)：
\[\mathbf R(\mathbf u,\lambda)=\mathbf0,
\]
再增加一个 path constraint：
\[g(\Delta\mathbf u,\Delta\lambda)=0.
\]
这样可以跨越 limit point，用于 snap-through、snap-back 和 post-buckling 等问题。
3. 总结
理解了上面所说的大框架，再看商业软件里的各种选项，就不应该只问"这个算法是什么"，而应该先问：它属于 nonlinear FEM的哪一层？它修改的是运动学、本构积分、单元离散，还是全局非线性求解？
浙公网安备 33010602011771号