Newton-Raphson 法¶

线性化与迭代递推式¶

用 Newton-Raphson 法求解外力虚功与整体方程组装中得到的、关于时刻 \(t_{n+1}\) 节点位移 \(\boldsymbol{u}_{n+1}\) 的非线性方程。设时刻 \(t_n\) 的节点位移 \(\boldsymbol{u}_n\) 已知，以位移增量 \(\Delta\boldsymbol{u}\) 为未知量，求

\[ \boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u} \]

。以下忽略外力向量对节点位移的依赖，令 \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\)，求解

\[ \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

。

在当前解 \(\Delta\boldsymbol{u}\) 处的切线刚度

\[ \boldsymbol{K} = \left. \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \right|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}} \]

用于线性化非线性方程后，得到

\[ \boldsymbol{K}\, d\boldsymbol{u} + \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

。将第 \(i\) 次迭代的修正量记为 \(d\boldsymbol{u}_i\)，将迭代开始时的残差向量记为

\[ \boldsymbol{R}_{i-1} = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) \]

则迭代递推式为

\[ \boldsymbol{K}_i\, d\boldsymbol{u}_i = \boldsymbol{R}_{i-1}, \qquad \Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i \]

。残差 \(\boldsymbol{R}_i\) 表示相对于平衡状态的力不平衡。

切线刚度矩阵的构造¶

切线刚度 \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) 通过对内力虚功的离散化中得到的单元内力向量关于节点位移求偏导，并将各单元被积项在单元区域上积分后汇集构成。将单元层面的被积项记为 \(\boldsymbol{K}^e_X\)（参考构形表示，TL 法）或 \(\boldsymbol{K}^e_x\)（当前构形表示，UL 法），则单元切线刚度为

\[ \boldsymbol{K}^e = \int_{\Omega^e_0} \boldsymbol{K}^e_X\, dV \quad (\text{TL}), \qquad \boldsymbol{K}^e = \int_{\Omega^e} \boldsymbol{K}^e_x\, dv \quad (\text{UL}) \]

。下面给出 TL/UL 被积项的最终形式。两者均可分解为材料刚度项（初始位移项）与几何刚度项（初始应力项）之和。

Total Lagrange 法¶

Total Lagrange 法假定第 2 Piola-Kirchhoff 应力率 \(\dot{\boldsymbol{S}}\) 与 Green-Lagrange 应变率 \(\dot{\boldsymbol{E}}\) 之间满足线性关系 \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\)。这对应线性弹性体（St.Venant-Kirchhoff 体）和超弹性体的本构关系；FrontISTR 对这些材料采用 Total Lagrange 法。此时单元切线刚度的被积项以张量形式写为

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_X\, \dot{\boldsymbol{u}}^e = \dot{\boldsymbol{S}}:\delta\boldsymbol{E} + \boldsymbol{S}:(\delta\boldsymbol{F}^T \dot{\boldsymbol{F}}) \]

。右端第 1 项为材料刚度项（初始位移项），第 2 项为几何刚度项（初始应力项）。

在 FrontISTR 的实现中，该被积项使用 Voigt 表示按矩阵形式

\[ \boldsymbol{K}^e_X = (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \tilde{\boldsymbol{C}}\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) + \boldsymbol{F}_9^T\, \boldsymbol{S}_9\, \boldsymbol{F}_9 \]

进行计算。各矩阵如下。\(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) 是内力虚功的离散化中引入的 B 矩阵，\(\tilde{\boldsymbol{C}}\) 是本构张量 \(\boldsymbol{\mathsf{C}}\) 的 Voigt 表示所对应的 \(6\times 6\) 材料刚度矩阵（见张量记法与数学基础）。\(\boldsymbol{S}_9, \boldsymbol{F}_9\) 是用于将几何刚度项表示为矩阵乘积的重排矩阵。首先定义将 \(3\times 3\) 二阶张量 \(\boldsymbol{A}\) 重排成 9 分量向量的记号 \([\,\cdot\,]\)：

\[ [\boldsymbol{A}] = (A_{11}, A_{21}, A_{31}, A_{12}, A_{22}, A_{32}, A_{13}, A_{23}, A_{33})^T \]

。由此，\(\boldsymbol{F}_9\) 是将变形梯度的变分写成 \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) 形式的 \(9\times d n_e\) 矩阵；与单元节点 \(\alpha = 1, \ldots, n_e\) 对应的 \(9\times d\) 块为

\[ [\boldsymbol{F}_9]_\alpha = \begin{bmatrix} (\partial N_\alpha^e/\partial X_1)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_2)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_3)\, \boldsymbol{I} \end{bmatrix} \qquad (\boldsymbol{I} \text{ 为 } 3\times 3 \text{ 单位矩阵}) \]

，按单元节点顺序横向排列，得到 \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\)。\(\boldsymbol{S}_9\) 是与其组合后，使几何刚度项可表示为 \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\) 的 \(9\times 9\) 矩阵

\[ \boldsymbol{S}_9 = \begin{bmatrix} S_{11} \boldsymbol{I} & S_{12} \boldsymbol{I} & S_{13} \boldsymbol{I} \\ S_{21} \boldsymbol{I} & S_{22} \boldsymbol{I} & S_{23} \boldsymbol{I} \\ S_{31} \boldsymbol{I} & S_{32} \boldsymbol{I} & S_{33} \boldsymbol{I} \end{bmatrix} \]

。

Updated Lagrange 法¶

Updated Lagrange 法假定相对 Kirchhoff 应力张量的 Jaumann 率 \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) 与变形速率张量 \(\boldsymbol{D}\) 满足线性关系 \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\)。这是线性弹性体、弹塑性体和蠕变材料共同采用的次弹性本构形式，FrontISTR 对这些材料使用 Updated Lagrange 法。此时以当前构形表示的单元切线刚度被积项，以张量形式写为

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_x\, \dot{\boldsymbol{u}}^e = \boldsymbol{\sigma}^{\nabla T}:\delta\boldsymbol{A}_{(L)} + \boldsymbol{\sigma}:(\delta\boldsymbol{F}_t^T\, \boldsymbol{L}) \]

（\(\boldsymbol{\sigma}^{\nabla T}\) 为 Truesdell 率，\(\boldsymbol{A}_{(L)}\) 为 Almansi 应变的线性部分，\(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) 为相对于当前构形的位移梯度，\(\boldsymbol{L}\) 为速度梯度张量）。右端第 1 项为材料刚度项，第 2 项为几何刚度项。

在 FrontISTR 的实现中，该被积项使用 Voigt 表示按矩阵形式

\[ \boldsymbol{K}^e_x = \boldsymbol{b}^T\, (\tilde{\boldsymbol{C}} - \boldsymbol{G})\, \boldsymbol{b} + \boldsymbol{f}_9^T\, \boldsymbol{\sigma}_9\, \boldsymbol{f}_9 \]

进行计算。\(\boldsymbol{b}\) 是在当前构形中构造的 B 矩阵（内力虚功的离散化）。\(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) 是在 TL 法定义的 \(\boldsymbol{S}_9, \boldsymbol{F}_9\) 中，将第 2 PK 应力 \(\boldsymbol{S}\) 替换为 Cauchy 应力 \(\boldsymbol{\sigma}\)，将参考构形梯度 \(\partial N_\alpha^e/\partial X_i\) 替换为当前构形梯度 \(\partial N_\alpha^e/\partial x_i\) 所得。

\(\boldsymbol{G}\) 是一个依赖 Cauchy 应力的修正矩阵，用于使次弹性本构关系 \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) 作为基于 Truesdell 率的本构关系与切线刚度框架一致。它将四阶张量分量 \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) 排列为 \(6\times 6\) Voigt 形式：

\[ \boldsymbol{G} = \begin{bmatrix} 2\sigma_{11} & 0 & 0 & \sigma_{12} & 0 & \sigma_{31} \\ 0 & 2\sigma_{22} & 0 & \sigma_{12} & \sigma_{23} & 0 \\ 0 & 0 & 2\sigma_{33} & 0 & \sigma_{23} & \sigma_{31} \\ \sigma_{12} & \sigma_{12} & 0 & \tfrac{\sigma_{11}+\sigma_{22}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{23}}{2} \\ 0 & \sigma_{23} & \sigma_{23} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{22}+\sigma_{33}}{2} & \tfrac{\sigma_{12}}{2} \\ \sigma_{31} & 0 & \sigma_{31} & \tfrac{\sigma_{23}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{33}+\sigma_{11}}{2} \end{bmatrix} \]

。

整体刚度矩阵的组装¶

整体切线刚度 \(\boldsymbol{K}\) 通过将单元刚度 \(\boldsymbol{K}^e\) 分成节点对对应的 \(d\times d\) 块 \(\boldsymbol{K}^e_{\alpha\beta}\)，并使用单元节点物理量的组装中引入的二阶张量组装集合 \(\mathcal{E}^2(i_g, i_h)\)，按下式构造：

\[ \boldsymbol{K}_{i_gi_h} = \sum_{(e,\alpha,\beta) \in \mathcal{E}^2(i_g, i_h)} \boldsymbol{K}^e_{\alpha\beta} \]

所得值排列在第 \(i_g\) 行、第 \(i_h\) 列形成矩阵。在实现中不显式构造集合 \(\mathcal{E}^2\)，而是在单元循环中直接累加到相应块。矩阵是阶数为“每节点自由度 \(\times\) 整体节点数 \(n_g\)”的方阵，但由于除通过单元相连的节点之间外，其余分量均为 \(0\)，因此采用稀疏矩阵形式存储。

TL 法与 UL 法的单元刚度除参考构形的切换（节点坐标与 B 矩阵的构造来源）以及是否包含 \(\boldsymbol{G}\) 矩阵外具有相同形式，因此 FrontISTR 使用公共子程序实现两者。

迭代算法¶

综上，在迭代开始时令 \(\Delta\boldsymbol{u} = \boldsymbol{0}\)，计算初始残差 \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\)，然后在第 \(i\) 次迭代中执行以下步骤。

按切线刚度矩阵的构造中的过程，计算当前位移 \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) 处的切线刚度 \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\)。
为施加几何边界条件，对受位移约束的自由度修改切线刚度矩阵和残差向量，得到 \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\)（见几何边界条件的处理）。
求解线性方程 \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\)，得到修正量 \(d\boldsymbol{u}_i\)。该过程通常占迭代计算的大部分计算开销。
将位移增量更新为 \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\)，并据此计算内力向量 \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) 及残差 \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\)。
进行收敛判定，若已达到收敛则结束迭代。残差 \(\boldsymbol{R}_i\) 中，施加几何边界条件的自由度会出现相当于约束反力的分量，因此使用去除这些分量后的 \(\tilde{\boldsymbol{R}}_i\) 构造判定指标。具体判定指标和阈值见收敛判定。若未收敛并达到迭代上限，则视为迭代失败。

迭代收敛后，将此时的 \(\Delta\boldsymbol{u}\) 加到 \(\boldsymbol{u}_n\)，得到时刻 \(t_{n+1}\) 的累计位移 \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\)，并进入下一个时间步。