一致性切线模量
一致性切线模量（也称为算法切线模量）反映了数值算法中应力对总应变的线性化增量关系，用于非线性有限元分析 Jacobian 矩阵的计算中，以显著提高全局非线性迭代的收敛速度和数值稳定性。因此，在本构积分后，计算并返回一致性切线模量是实现高效、稳健的有限元求解的关键步骤
一致性切线模量定义为
\[
\mathbf{D}_{ep}^{\text{alg}}=\frac{\partial \boldsymbol{\sigma}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}}.
\]
下简记为 \(\mathbf{D}\)
对于弹塑性材料，在已知前一步内变量 \(\boldsymbol{\alpha}_{n}\) 和当前步总应变 \(\boldsymbol{\varepsilon}_{n+1}\) 的情况下，通常通过本构积分算法来更新应力。这个过程自然定义了一个算子形式的增量本构函数 \(\hat{\boldsymbol{\sigma}}\)，即
\[
\boldsymbol{\sigma}_{n+1}=\hat{\boldsymbol{\sigma}}(\boldsymbol{\alpha}_{n},\boldsymbol{\varepsilon}_{n+1})
\]
有限应变塑性理论中，常常使用乘子应变分解，即
\[
\mathbf{F} = \mathbf{F}^{e}\cdot\mathbf{F}^{p},
\]
其中
因此在有限应变塑性中，不定义总非线性应变 \(\boldsymbol{\varepsilon}_{n+1}\)，但弹性试应变 \(\boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}}\) 仍自然出现在弹性预测-塑性校正算法中。由于
\[
\boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}} = \boldsymbol{\varepsilon}_{n}^{e} + \Delta\boldsymbol{\varepsilon}_{n+1} = \boldsymbol{\varepsilon}_{n} - \boldsymbol{\varepsilon}_{n}^{p} + \Delta\boldsymbol{\varepsilon}_{n+1} = \boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_{n}^{p},
\]
因此，将 \(\boldsymbol{\sigma}_{n+1}\) 表示为 \(\boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}}\) 的函数
\[
\boldsymbol{\sigma}_{n+1}=\bar{\boldsymbol{\sigma}}(\boldsymbol{\alpha}_{n},\boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}})=\hat{\boldsymbol{\sigma}}(\boldsymbol{\alpha}_{n},\boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}}+\boldsymbol{\varepsilon}_{n}^{p})=\hat{\boldsymbol{\sigma}}(\boldsymbol{\alpha}_{n},\boldsymbol{\varepsilon}_{n+1}),
\]
此时有
\[
\mathbf{D} = \frac{\partial \hat{\boldsymbol{\sigma}}}{\partial \boldsymbol{\varepsilon}_{n+1}} = \frac{\partial \hat{\boldsymbol{\sigma}}}{\partial \boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}}} = \frac{\partial \bar{\boldsymbol{\sigma}}}{\partial \boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}}}.
\]
为了使小应变塑性模型和有限应变塑性模型拥有形式一致的一致性切线模量，我们采用如下形式
\[
\mathbf{D} = \frac{\partial \hat{\boldsymbol{\sigma}}}{\partial \boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}}} = \frac{\partial \bar{\boldsymbol{\sigma}}}{\partial \boldsymbol{\varepsilon}_{n+1}^{e,\text{trail}}}.
\]
分段性
\(\hat{\boldsymbol{\sigma}}\) 是一个分段函数，当
\(\Phi^{\text{trial}} < 0\)，则 \(\hat{\boldsymbol{\sigma}}\) 随着弹性曲线演化，此时 \(\mathbf{D} = \mathbf{D}^{e}\)
\(\Phi^{\text{trial}} > 0\)，则 \(\hat{\boldsymbol{\sigma}}\) 随着塑性曲线演化，此时 \(\mathbf{D} = \mathbf{D}^{ep}\)
\(\Phi^{\text{trial}} = 0\)，则可能发生弹性卸载，或发生塑性演化，在切换点处，\(\hat{\boldsymbol{\sigma}}\) 不可微
各向同性硬化的 Mises 屈服准则
接下来以 Mises 屈服准则和各向同性硬化准则为例，介绍一致性切线模量 \(\mathbf{D}\) 的计算过程，根据 return mapping 算法，有
\[
\begin{equation}
\begin{aligned}
\boldsymbol{\sigma}_{n+1}=\left(\mathbf{D}^e - \frac{6G^{2}\Delta\gamma}{q^{\text{trial}}}\mathbf{I}_{d}\right) : \boldsymbol{\varepsilon}^{e, \text{trial}},
\end{aligned}
\end{equation}
\]
其中，\(q^{\text{trial}}=\sqrt{3J_{2}(\mathbf{s}^{\text{trial}})} = \sqrt{\frac{3}{2}}\|\mathbf{s}^{\text{trial}}\|\)，\(\Delta\gamma\) 是方程
(93)\[
\begin{equation}
q^{\text{trial}}-3G\Delta\gamma-\sigma_{y}(\bar{\varepsilon}_n^p + \Delta \gamma)=0
\end{equation}
\]
的解，于是
(94)\[\begin{split}
\begin{equation}
\begin{aligned}
\frac{\partial \boldsymbol{\sigma}_{n+1}}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}
=
\mathbf{D}^e
-
\frac{\Delta \gamma \, 6 G^2}{q^{\text{trial}}} \mathbf{I}_d
-
\frac{6 G^2}{q^{\text{trial}}} \boldsymbol{\varepsilon}^{e, \text{trial}}_{d}
\otimes
\frac{\partial \Delta \gamma}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}
\\
+
\frac{\Delta \gamma \, 6 G^2}{(q^{\text{trial}})^2} \boldsymbol{\varepsilon}^{e, \text{trial}}_{d}
\otimes
\frac{\partial q^{\text{trial}}}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}},
\end{aligned}
\end{equation}
\end{split}\]
其中，\(\boldsymbol{\varepsilon}^{e, \text{trial}}_{d} = \text{dev}(\boldsymbol{\varepsilon}^{e, \text{trial}}) = \mathbf{I}_{d}:\boldsymbol{\varepsilon}^{e, \text{trial}}\)
由于
(95)\[
q^{\text{trial}} = \sqrt{\frac{3}{2}}\|\mathbf{s}^{\text{trial}}\| = 2G\sqrt{\frac{3}{2}}\|\mathbf{I}_{d}:\boldsymbol{\varepsilon}^{e,\text{trial}}\|=2G\sqrt{\frac{3}{2}}\|\boldsymbol{\varepsilon}^{e, \text{trial}}_{d}\|,
\]
故
(96)\[
\begin{equation}
\frac{\partial q^{\text{trial}}}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}=2G\sqrt{\frac{3}{2}}\frac{\mathbf{s}^{\text{trial}}}{\|\mathbf{s}^{\text{trial}}\|}=2G\sqrt{\frac{3}{2}}\frac{\boldsymbol{\varepsilon}^{e, \text{trial}}_{d}}{\|\boldsymbol{\varepsilon}^{e, \text{trial}}_{d}\|}.
\end{equation}
\]
此外，根据式 (93)，有
(97)\[
\frac{\partial q^{\text{trial}}}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}-3G\frac{\partial \Delta\gamma}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}-\frac{\partial \sigma_{y}(\bar{\varepsilon}_n^p + \Delta \gamma)}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}=0,
\]
一般地，可从上式求解出 \(\frac{\partial \Delta\gamma}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}\)，对于线性硬化模型，有
\[
\sigma_{y}(\bar{\varepsilon}_n^p + \Delta \gamma) = \sigma_{0} + H\cdot(\bar{\varepsilon}_n^p + \Delta \gamma),
\]
得到
\[
\frac{\partial \sigma_{y}(\bar{\varepsilon}_n^p + \Delta \gamma)}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}} = H\frac{\partial \Delta\gamma}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}},
\]
故，代入式 (97) 中，得到
(98)\[
\begin{equation}
\frac{\partial \Delta\gamma}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}} = \frac{1}{3G+H}\frac{\partial q^{\text{trial}}}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}} = \frac{2G}{3G+H}\sqrt{\frac{3}{2}}\frac{\mathbf{s}^{\text{trial}}}{\|\mathbf{s}^{\text{trial}}\|}.
\end{equation}
\]
记
\[
\mathbf{\bar{N}} = \frac{\mathbf{s}^{\text{trial}}}{\|\mathbf{s}^{\text{trial}}\|} = \frac{\boldsymbol{\varepsilon}^{e, \text{trial}}_{d}}{\|\boldsymbol{\varepsilon}^{e, \text{trial}}_{d}\|},
\]
将式 (96) 和式 (98) 代入到 (94)，并结合式 (95) 得到
\[
\begin{equation}
\begin{aligned}
\mathbf{D}=
\frac{\partial \boldsymbol{\sigma}_{n+1}}{\partial \boldsymbol{\varepsilon}^{e, \text{trial}}}
&=
\mathbf{D}^e
-
\frac{\Delta \gamma \, 6 G^2}{q^{\text{trial}}} \mathbf{I}_d
+
6G^2\left(\frac{\Delta\gamma}{q^{\text{trial}}}-\frac{1}{3G+H}\right)\mathbf{\bar{N}}\otimes\mathbf{\bar{N}},
\end{aligned}
\end{equation}
\]
此处，\(\mathbf{D}\) 是对称的
连续切线模量
对于各向同性硬化
\[
\sigma_{y}(\bar{\varepsilon}^p) = \sigma_{0} + \kappa(\bar{\varepsilon}^p),
\]
有
\[
\boldsymbol{\alpha} = \{\bar{\varepsilon}_{n}^{p}\},\quad \mathbf{A} = \{\kappa\},
\]
故
\[
\begin{equation}
\frac{\partial^2 \psi^p}{\partial \boldsymbol{\alpha}^2}=\frac{\partial \psi^p}{\partial {\bar{\varepsilon}^{p}}^{2}}=\frac{\partial \kappa}{\partial \bar{\varepsilon}^{p}} = H(\bar{\varepsilon}^{p}),
\end{equation}
\]
对于线性硬化，有 \(H(\bar{\varepsilon}^{p}) \equiv H\). 此外，对于 Mises 屈服准则，有
\[
\mathbf{H} = -\frac{\partial \Phi}{\partial \mathbf{A}} = -\frac{\partial \Phi}{\partial \kappa} = 1
\]
且对于关联塑性流动，有
\[
\mathbf{N} = \frac{\partial \Phi}{\partial \boldsymbol{\sigma}},
\]
于是
\[\begin{split}
\begin{equation}
\begin{aligned}
\mathbf{D}^{ep}_{c} &= \mathbf{D}^{e} - \frac{(\mathbf{D}^{e}:\mathbf{N})\otimes(\mathbf{D}^{e}:\frac{\partial \Phi}{\partial \boldsymbol{\sigma}})}{\frac{\partial \Phi}{\partial \boldsymbol{\sigma}}:\mathbf{D}^{e}:\mathbf{N}-\frac{\partial \Phi}{\partial \mathbf{A}}\frac{\partial^2 \psi^p}{\partial \boldsymbol{\alpha}^2}\mathbf{H}}\\
&=\mathbf{D}^{e} - \frac{(\mathbf{D}^{e}:\mathbf{N})\otimes(\mathbf{D}^{e}:\mathbf{N})}{\mathbf{N}:\mathbf{D}^{e}:\mathbf{N}+H}
\end{aligned}
\end{equation}
\end{split}\]
由于 \(\mathbf{N}\) 是偏张量，故
\[
\mathbf{D}^{e}:\mathbf{N} = 2G\mathbf{N}\quad \Longrightarrow\quad \mathbf{N}:\mathbf{D}^{e}:\mathbf{N} = 3G,
\]
最终，连续切线模量化简为
\[
\mathbf{D}^{ep}_{c} =\mathbf{D}^{e} - \frac{6G^2}{3G+H}\bar{\mathbf{N}}\otimes\bar{\mathbf{N}}.
\]
连续切线模量与一致性切线模量关系如下
\[
\mathbf{D}=\mathbf{D}_{c}^{ep}-\frac{\Delta \gamma \, 6 G^2}{q^{\text{trial}}}\left[\mathbf{I}_d-\bar{\mathbf{N}}\otimes\bar{\mathbf{N}}\right],
\]
当 \(\Delta\gamma\rightarrow0\) 时，\(\mathbf{D}_{c}^{ep}\rightarrow\mathbf{D}\)；然而，当 \(\Delta\gamma\) 较大（例如时间步长较大时），\(\mathbf{D}{c}^{ep}\) 将逐渐偏离 \(\mathbf{D}\)。如果此时仍继续采用 \(\mathbf{D}_{c}^{ep}\)，可能会导致明显的数值不一致性，从而影响牛顿迭代的收敛性，甚至导致收敛速度显著降低
一般情形
在隐式 Return Mapping 算法中，需满足
(99)\[\begin{split}
\begin{equation}
\left\{
\begin{array}{l}
\boldsymbol{\varepsilon}^e_{n+1} - \boldsymbol{\varepsilon}^{e,\text{trial}} + \Delta\gamma\, \mathbf{N}_{n+1} \\
\boldsymbol{\alpha}_{n+1} - \boldsymbol{\alpha}_n - \Delta\gamma\, \mathbf{H}_{n+1} \\
\Phi(\boldsymbol{\sigma}_{n+1}, \mathbf{A}_{n+1})
\end{array}
\right\}
=
\left\{
\begin{array}{l}
\mathbf{0} \\
\mathbf{0} \\
0
\end{array}
\right\},
\end{equation}
\end{split}\]
求解变量为
\[
\boldsymbol{\varepsilon}^{e}_{n+1},\quad\boldsymbol{\alpha},\quad\Delta\gamma.
\]
对方程 (99) 两边微分，得到（为了方便，将下标 \(n+1\) 隐去）
\[\begin{split}
\begin{equation}
\left\{
\begin{array}{l}
\mathrm{d}\boldsymbol{\varepsilon}^e + \Delta\gamma \dfrac{\partial \mathbf{N}}{\partial \boldsymbol{\sigma}} : \mathrm{d}\boldsymbol{\sigma} + \Delta\gamma \dfrac{\partial \mathbf{N}}{\partial \mathbf{A}} \mathrm{d}\mathbf{A} + \mathrm{d}\Delta\gamma\, \mathbf{N} \\[2ex]
\mathrm{d}\boldsymbol{\alpha} - \Delta\gamma \dfrac{\partial \mathbf{H}}{\partial \boldsymbol{\sigma}} \mathrm{d}\boldsymbol{\sigma} - \Delta\gamma \dfrac{\partial \mathbf{H}}{\partial \mathbf{A}} \mathrm{d}\mathbf{A} - \mathrm{d}\Delta\gamma\, \mathbf{H} \\[2ex]
\dfrac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \mathrm{d}\boldsymbol{\sigma} + \dfrac{\partial \Phi}{\partial \mathbf{A}} \mathrm{d}\mathbf{A}
\end{array}
\right\}
=
\left\{
\begin{array}{l}
\mathrm{d}\boldsymbol{\varepsilon}^{e,\text{trial}} \\
\mathbf{0} \\
0
\end{array}
\right\},
\end{equation}
\end{split}\]
整理得到
\[\begin{split}
\begin{equation}
\begin{bmatrix}
\mathrm{d}\boldsymbol{\sigma} \\
\mathrm{d}\mathbf{A} \\
\mathrm{d}\Delta\gamma
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{D}_{11} & \mathbf{D}_{12} & \mathbf{D}_{13} \\
\mathbf{D}_{21} & \mathbf{D}_{22} & \mathbf{D}_{23} \\
\mathbf{D}_{31} & \mathbf{D}_{32} & D_{33}
\end{bmatrix}
\begin{bmatrix}
\mathrm{d}\boldsymbol{\varepsilon}^{e,\mathrm{trial}} \\
\mathbf{0} \\
0
\end{bmatrix}.
\end{equation}
\end{split}\]
于是
\[
\mathbf{D}_{11} = \frac{\mathrm{d} \boldsymbol{\sigma}_{n+1}}{\mathrm{d} \boldsymbol{\varepsilon}^{e,\mathrm{trial}}},\quad
\mathbf{D}_{21} = \frac{\mathrm{d} \mathbf{A}_{n+1}}{\mathrm{d} \boldsymbol{\varepsilon}^{e,\mathrm{trial}}},\quad
\mathbf{D}_{31} = \frac{\mathrm{d} \Delta\gamma}{\mathrm{d} \boldsymbol{\varepsilon}^{e,\mathrm{trial}}}.
\]
