onigiri_note’s diary

概要

弾塑性FEMの陰的アルゴリズム（俗にいうリターンマッピング）を実装するにあたり最低限必要と思われる事柄について説明。
リターンマッピングは弾塑性FEMの計算におけるどの部分の処理にあたるかについても説明。
弾性FEMまでに必要な事項は簡単に振り返る程度に留める。
今回の一連の記事は等方性材料、微小変形、微小回転を仮定している。

前提となる知識の準備
　構造解析(静解析)における有限要素法(以下、単にFEM)では、式(1)のような方程式を $u$ について解いている。 $K_{all}$ は全体剛性マトリクス、 $u$ は変位ベクトル、 $f_{all}$ は全体荷重ベクトルと呼ばれる。 $f_{all}$ は境界条件から決定される。極々、簡単に言えば、ばね定数 $K_{all}$ と外力 $f_{all}$ から物体の変形 $u$ を求めるということをFEMでは行っている。
$\displaystyle K_{all} \cdot u = f_{all} \tag{1}$
$K_{all}$ は式(2)のように要素剛性マトリクス $K_{e}$ の足し合わせで与えられ、接線係数 $D$ とBマトリクス $B$ から決まる。ただし、式(2)における $\sum$ 記号は要素アッセンブリー作用素の意味で用いている。接線係数 $D$ は構成則（応力‐ひずみ関係）から決まる物性値でDマトリクスとも呼ばれ、１軸問題における応力‐ひずみ線図の接線(Fig.1)に対応する。Bマトリクスは形状関数の偏微分で表現され、多項式になる。Bマトリクスの具体的な形については、ググったらすぐ出てくるので気になる人は調べてみてください。tex形式で数式を書くのが大変だったので省略します。ごめんなさい。
$\displaystyle K_{all} = \sum_{elem}{K_e} = \sum_{elem}{ \int_{V_e} B^T D B {\rm d}V} \tag{2}$

f:id:onigiri_note:20181230141755j:plain — Fig.1 応力‐ひずみ線図

弾性変形範囲では、Dマトリクスはひずみ（変位）に依存しない定数となり、式(2)から計算される全体剛性マトリクス $K_{all}$ は変位 $u$ に依存しない定数となる。よって、弾性変形範囲で式(1)は $u$ についての連立１次方程式となり、連立１次方程式を解くプログラムを持っていれば解くことができる。一方、塑性変形域ではDマトリクスがひずみ（変位）に依存するため、 $K_{all}$ も変位 $u$ に依存する。つまり、塑性変形域では、式(1)は $u$ についての非線形方程式となり、あからさまに書くと式(3)のようになる。
$\displaystyle K_{all}(u) \cdot u = f_{all} \tag{3}$
ここで、式(3)を式(4)~(6)のように離散化し、未知数 $\Delta u_{n+1}$ を未知数として反復法で解くことを考えてみる。nステップまでの値と境界条件から決まる $\Delta f_{n+1}$ は既知の値とする。例えば、Fig.2のような反復法のアルゴリズムで $\Delta u_{n+1}$ が求められる。
$\displaystyle K_{all}(u_{n+1}) \cdot u_{n+1} = {f_{all}}_{n+1} \tag{4}$

$\displaystyle u_{n+1} = u_n + \Delta u_{n+1} \tag{5}$

$\displaystyle \Delta f_{n+1} = {f_{all}}_{n+1} - {f_{all}}_{n} \tag{6}$

f:id:onigiri_note:20181230223529j:plain — Fig.2 変位増分を求める反復法アルゴリズム

リターンマッピングの概要と陽解法との対比
　弾塑性FEMにおける陰的アルゴリズム（以下、リターンマッピング）はFig.2のどの箇所に該当し、どのような処理を行っているかを陽解法のアルゴリズムと対比しながら簡単に説明する。
　リターンマッピングではFig.2における②の行程を陰解法的に処理する。つまり、リターンマッピングとは接線係数（Dマトリクス）を陰解法的に求める処理に相当する。通常、定式化の過程が複雑になりすぎないように、式(7)のようなオイラーの後退差分が利用される（n+1ステップ目の接線係数をn+1ステップ目の情報も利用して求める）。また、「降伏曲面上に応力が乗る」という条件に相当する式(8)も合わせて利用される。ここで、 $\bar{\sigma}_{n+1}$ はMisesの相当応力、 $F(p_{n+1})$ は降伏応力、 $p_{n+1}$ は累積相当塑性ひずみをあらわす。「降伏曲面上に応力が乗る」という静力学の基本的な条件をリターンマッピングの処理は満たすので、リターンマッピングで求めた接線係数Dはコンシステント接線係数（整合接線係数）と呼ばれる。式(7)と式(8)を連成させて接線係数Dを求めるわけだが、ごく簡単な構成式を除いて、通常は非線形方程式となっているので、反復法によりコンシステント接線係数を求める。１軸問題におけるリターンマッピングの処理のイメージをFig.3に図示する。
※式(7)のような表現は本来誤りであるが、リターンマッピングの処理のイメージを掴んで頂くためにこのような表現をした。この辺りの細かいことについては続きの記事で解説する。
$\displaystyle D = \left. \frac{\partial \sigma}{\partial \epsilon} \right|_{u=u_{n+1}} \tag{7}$
$\displaystyle \bar{\sigma}_{n+1} - F(p_{n+1}) = 0 \tag{8}$

f:id:onigiri_note:20181230235449j:plain — Fig.3 リターンマッピングのイメージ図

　次に、陽解法で接線係数（Dマトリクス）を求める場合はどのような処理をしているかを説明する。陽解法ではn+1ステップ目の接線係数をnステップ目の情報だけから求める。陰解法（リターンマッピング）との対比のために、陽解法における接線係数の求め方のイメージ図をFig.4に図示する。ここで、式(8)に相当する「降伏局面に応力が乗る」という条件は利用されないため、構成式から決まる応力‐ひずみ関係（力‐変位関係）を満たさない接線係数が求められてしまう。このことから、Fig.2における③の行程の全体剛性マトリクスも本来設定していた力‐変位関係と矛盾した値になってしまう。これが原因でFig.2で示した反復計算の収束性が陰解法（リターンマッピング）に比べて著しく低下してしまう。条件によっては陰解法に比べて陽解法は100倍以上の計算時間がかかってしまう場合もある。計算時間短縮のために、定式化はやや複雑になるが汎用のFEMコードでは陰解法（リターンマッピング）が採用されている。
$\displaystyle D = \left. \frac{\partial \sigma}{\partial \epsilon} \right|_{u=u_{n}} \tag{9}$

f:id:onigiri_note:20181231142156j:plain — Fig.4 陽解法における接線係数の求め方のイメージ図

まとめ

リターンマッピングとはコンシステント接線係数を求める処理のこと。
コンシステント接線係数を用いると変位ベクトルを求める際の収束性が著しく向上する。

今後の記事の方針

等方硬化と移動硬化に従う場合のリターンマッピングの定式化を具体的な構成式を用いて行う。
余裕があれば、フリーのFEMソフトのユーザー関数を利用して、リターンマッピングを自分で実装する際のプログラムも公開する。