Newton-Raphson 法¶

線形化と反復漸化式¶

外力仮想仕事と全体方程式の組立てで得た，時刻 \(t_{n+1}\) における節点変位 \(\boldsymbol{u}_{n+1}\) に関する非線形方程式を Newton-Raphson 法によって解く．時刻 \(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\) の 2 階テンソル \(\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}\) は，亜弾性構成則 \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) を Truesdell 速度ベースの構成則として接線剛性の枠組みに整合させるために必要な Cauchy 応力依存の補正行列であり，4 階テンソル成分 \(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}\) に分け，要素節点物理量のアセンブルで導入した 2 階テンソル版アセンブル集合 \(\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}\) とし，次の時間ステップへ進む．