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}}) \]

。右側第一項為材料剛性項（初始位移項），第二項為幾何剛性項（初始應力項）。

在 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}\) 為速度梯度張量。右側第一項為材料剛性項，第二項為幾何剛性項。

在 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}\)，並進入下一個時間步。