微小ひずみの可塑性¶

GetFEM における可塑性モデルの近似のフレームワークについて説明します．ブリック要素を実装し新しい可塑性モデルを拡張する場合は interface/src/gf_model_set.cc と interface/src/gf_model_set.cc を参照してください．

理論的背景¶

微小ひずみの可塑性について簡単に紹介します．より詳細な説明については，主に [SI-HU1998] および [SO-PE-OW2008] を参照してください．

微小ひずみテンソルの加算分解¶

\(\Omega \subset \rm I\hspace{-0.15em}R^3\) は変形可能な物体の基準配置であり， \(u : \Omega \rightarrow \rm I\hspace{-0.15em}R^3\) は変形フィールドです．微小ひずみの可塑性は，微小ひずみテンソル \(\varepsilon(u) = \dfrac{\nabla u + \nabla u^T}{2}\) の加法分解に基づいています．

\[\varepsilon(u) = \varepsilon^e + \varepsilon^p\]

ここで， \(\varepsilon^e\) はひずみテンソルの弾性部分であり， \(\varepsilon^p\) は弾性部分です．

内部変数，自由エネルギーポテンシャル，弾性法則¶

以下を考慮します

\[\alpha : \Omega \rightarrow \rm I\hspace{-0.15em}R^{d_{\alpha}},\]

これはひずみ型の内部変数 \(d_{\alpha}\) のベクトルフィールドです．(内部変数が考慮されていない場合は，\(d_{\alpha} = 0\) となります．)自由エネルギーポテンシャルも考慮します．

\[\psi(\varepsilon^e, \alpha),\]

対応する応力型変数は以下で決定されます．

\[\sigma = \dfrac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^e, \alpha), ~~~~ A = \dfrac{\partial \psi}{\partial \alpha}(\varepsilon^e, \alpha),\]

ここで， \(\sigma\) はCauchyの応力テンソルであり， \(A\) は応力型の内部変数です．降伏条件は

\[\sigma:\dot{\varepsilon}^p - A.\dot{\alpha} \ge 0.\]

通常， \(\psi(\varepsilon^e, \alpha)\) は次のように分解されます．

\[\psi(\varepsilon^e, \alpha) = \psi^e(\varepsilon^e) + \psi^p(\alpha).\]

線形化弾性体の場合， \(\psi^e(\varepsilon^e) = \frac{1}{2} ({\cal A}\varepsilon^e) :\varepsilon^e\) は4次の弾性テンソルです．等方性線形弾性の場合，この式は \(\psi^e(\varepsilon^e) = \mu \mbox{dev}(\varepsilon^e) : \mbox{dev}(\varepsilon^e) + \frac{1}{2} K (\mbox{tr}(\varepsilon^e))^2\) のようになります．ここで， \(\mu\) はせん断弾性率で \(K = \lambda + 2\mu/3\) はBulk弾性率です．

塑性ポテンシャル，降伏関数，塑性流れ則¶

応力が臨界値に達すると，塑性降伏が起こると考えられます．この限界は，降伏関数 \(f(\sigma, A)\) と条件によって決定されます．

\[f(\sigma, A) \le 0.\]

表面 \(f(\sigma, A) = 0\) は，塑性変形が起こる降伏面です．

また，(両方の変数に関して凸である)塑性ポテンシャル \(\Psi(\sigma, A)\) を考慮してみましょう．塑性の流れ則と方向は次のように定義されます．

\[\dot{\varepsilon}^p = \gamma \dfrac{\partial \Psi}{\partial \sigma}(\sigma, A), ~~~~~~ \dot{\alpha} = -\gamma \dfrac{\partial \Psi}{\partial A}(\sigma, A),\]

付加的な相補性条件として以下が成り立ちます．

\[f(\sigma, A) \le 0, ~~~ \gamma \ge 0, ~~~ f(\sigma, A) \gamma = 0.\]

変数 \(\gamma\) は塑性乗数と呼ばれます．以下の場合に注意してください． \(\psi(\varepsilon^e, \alpha), f(\sigma, A) \mbox{ と } \Psi(\sigma, A)\) が微分可能でないとき，副差を使用しなければなりません．関連する可塑性は \(\Psi(\sigma, A) = f(\sigma, A)\) に対応します．

初期境界値問題¶

動的弾塑性問題の弱定式化は，任意の運動学的に許容される次のような試験関数 \(v\) で書くことができます．

\[\begin{split}\left| \begin{array}{l} \int_{\Omega} \rho \ddot{u}\cdot v + \sigma : \nabla v dx = \int_{\Omega} f\dot v dx + \int_{\Gamma_N} g\dot v dx, \\ u(0,x) = u_0(x), ~~~\dot{u}(0) = \mathrm{v}_0(x), \\ \varepsilon^p(0,x) = \varepsilon^p_0, ~~~ \alpha(0,x) = \alpha_0, \end{array} \right.\end{split}\]

\(u_0, \mathrm{v}_0, \varepsilon^p_0, \alpha_0\) は初期値で，\(f\) と \(g\) は領域の内部で， \(\Omega\) と境界の部分には \(\Gamma_N\) があります．

可塑性モデルは，しばしば， \(\rho \ddot{u}\) 項が無視される準静的問題に適用されることに注意してください．

時間ステップ \(\Delta t = t_{n+1} -t_n\) が与えられたとします，続いて時刻 \(t_n\) から \(t_{n+1}\) までの， \(u_n, \varepsilon^p_n \mbox{ と } \alpha_n\) 時刻 \(t_n\) の \(u(t_n), \varepsilon^p_n(t_n) \mbox{ と } \alpha(t_n)\) の近似について示します．選択された時間積分スキームの近似は（例えば，提案されたスキームの1つ過渡問題の積分のためのモデルツール）流れ則の積分に使用される時間積分スキームとは異なっても構いません（以下を参照）．

流れ則の積分¶

塑性流れ則は，独自の時間積分法で積分する必要があります．標準法の中で，後方Euler法では， \(\theta\) 法（または汎用された台形ルール）と汎用された中間点法が，この文脈で最も一般的に使用されます．ここでは， \(\theta\)- 法（\(\theta = 1\) は特殊なケースとして後方Euler法に対応します）を選択します．

\(u_{n+1}\) は，考慮する時間ステップの変位であり，前のステップの変位は \(u_{n}\) です．

塑性流れ則の積分のための \(\theta\)- 法は

(1)¶\[\varepsilon^p_{n+1} - \varepsilon^p_{n} = (1-\theta)\Delta t \gamma_n \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n}, A_{n}) + \theta \Delta t \gamma_{n+1} \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n+1}, A_{n+1}),\]

(2)¶\[\alpha_{n+1} - \alpha_n = -(1-\theta)\Delta t \gamma_n \dfrac{\partial \Psi}{\partial A}(\sigma_{n}, A_{n}) - \theta\Delta t \gamma_{n+1} \dfrac{\partial \Psi}{\partial A}(\sigma_{n+1}, A_{n+1}),\]

相補性条件は次の通りです．

\[f(\sigma_{n+1}, A_{n+1}) \le 0, ~~~ \gamma_{n+1} \ge 0, ~~~ f(\sigma_{n+1}, A_{n+1}) \gamma_{n+1} = 0.\]

ここで， \(0 < \theta \le 1\) は \(\theta\) - 法のパラメータです．\(\theta = 0\) は，可塑性の陽な積分を考慮しないため除外します．\(\theta = 1\) は後方Euler法に対応し，\(\theta = 1/2\) は2次整合法であるCrank-Nicolson（または台形法）法に対応します．時間ステップ \(n\) における量の相補性条件は，前の時間ステップにより処理されています（ \(\sigma_{n}, \alpha_n, \mbox{と } \gamma_n\) は，すでに決定しています）．

A solution would be to solve the whole problem with all the unknowns, that is \(u_{n+1}, \gamma_{n+1}, \varepsilon^p_{n+1} \mbox{ and } A_{n+1}\). This is of course possible but would be a rather expensive strategy because of the resulting high number of degrees of freedom. A classical strategy (the return mapping one for instance, see [SO-PE-OW2008] or the closest point projection one) consists in integrating locally the plastic flow on each Gauss point of the considered integration method separately, or more precisely to consider on each Gauss point the maps

\[ \begin{align}\begin{aligned}{\mathscr E}^p : (u_{n+1}, \zeta_n, \eta_n) \mapsto \varepsilon^p_{n+1}\\{\mathscr A} : (u_{n+1}, \zeta_{n}, \eta_n) \mapsto \alpha_{n+1}\end{aligned}\end{align} \]

\(\eta_n, \zeta_{n}\) により方程式 (1), (2) の右辺は次のようになります．

\[ \begin{align}\begin{aligned}\zeta_n = \varepsilon^p_{n} + (1-\theta)\Delta t \gamma_n \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n}, A_{n}) ,\\\eta_n = \alpha_n - (1-\theta)\Delta t \gamma_n \dfrac{\partial \Psi}{\partial A}(\sigma_{n}, A_{n})\end{aligned}\end{align} \]

これは，特に次のことを意味します \((\varepsilon^p_{n+1}, \alpha_{n+1}) = ({\mathscr E}^p(u_{n+1}, \zeta_n, \eta_n), {\mathscr A}(u_{n+1}, \zeta_{n}, \eta_n))\) は方程式 (1) と (2) の解です．これらの写像とそれらの接線モジュラス（通常，一貫した接線モジュラスと呼ばれる）は，Newton法を用いた問題のグローバルな解法で使用され，\(u_{n+1}\) は一意な残りの変数です．Return Mapping法の利点は，グローバル解の一意の変数が変位 \(u_{n+1}\) であることです．各Gauss点上の非線形解法がしばしば必要であり，通常はローカルNewton法で実行されます．

GetFEM では我々はReturn Mapping法と，主に [PO-NI2016] ， [SE-PO-WO2015] および [HA-WO2009] に触発された，より簡単な接線モジュラスを可能にする以下に展開される代替戦略の両方を提案します．それは \(u_{n+1}\) に関して，さらに未知数なものとして， \(\gamma_{n+1}\) （倍の数）の値を保存する必要があります．それからわかるように，これは，降伏関数のより汎用的な処理を可能にします．これは，この追加の未知数のスカラー場の単純さと引き換えとなっています．

まず，後ほど現れる（単純なローカル逆変換が可能になる） \(\alpha(\sigma_{n+1}, A_{n+1}) > 0\) という追加の（そしてオプションの）関数と新しい未知のスカラー場を考えます

\[\xi_{n+1} = \dfrac{\gamma_{n+1}}{\alpha(\sigma_{n+1}, A_{n+1})} ,\]

私たちの2つの主要な未知数は \(u_{n+1} \mbox{ と } \xi_{n+1}\) になりました．離散化された塑性流れ則の積分は，次のようになります．

(3)¶\[\varepsilon^p_{n+1} - \varepsilon^p_{n} = (1-\theta)\alpha(\sigma_n,A_n)\Delta t \xi_n \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n}, A_{n}) + \theta \alpha(\sigma_{n+1},A_{n+1}) \Delta t \xi_{n+1} \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n+1}, A_{n+1}),\]

(4)¶\[\alpha_{n+1} - \alpha_n = (1-\theta) \alpha(\sigma_n,A_n)\Delta t \xi_n \dfrac{\partial \Psi}{\partial A}(\sigma_{n}, A_{n}) + \theta \alpha(\sigma_{n+1},A_{n+1}) \Delta t \xi_{n+1} \dfrac{\partial \Psi}{\partial A}(\sigma_{n+1}, A_{n+1}),\]

(5)¶\[f(\sigma_{n+1}, A_{n+1}) \le 0, ~~~ \xi_{n+1} \ge 0, ~~~ f(\sigma_{n+1}, A_{n+1}) \xi_{n+1} = 0.\]

\(u_{n+1} \mbox{ と } \xi_{n+1}\) が与えられたら，以下の2つの写像を定義します

\[ \begin{align}\begin{aligned}\tilde{\mathscr E}^p : (u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}, \eta_n) \mapsto \varepsilon^p_{n+1}\\\tilde{\mathscr A} : (u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}, \eta_n) \mapsto \alpha_{n+1}\end{aligned}\end{align} \]

ここでペア \((\varepsilon^p_{n+1}, \alpha_{n+1}) = (\tilde{\mathscr E}^p(u_{n+1}, \theta \xi_{n+1}, \zeta_{n}, \eta_n), \tilde{\mathscr A}(u_{n+1}, \theta \xi_{n+1}, \zeta_{n}, \eta_n))\) は方程式 (3), (4) の解です ( (5) は考慮しません．)．後ほど，少なくとも単純な等方性の塑性流れの則では，これらの写像は単純な表現を持ち，場合によっては \(u_{n+1}\) に対し線形になることを示しています．

その場合でも \(u_{n+1} \mbox{ と } \xi_{n+1}\) は応力 \(\sigma_{n+1}\) を与えます．

\[\sigma_{n+1} = \dfrac{\partial \psi^e}{\partial \varepsilon^e}(\varepsilon(u_{n+1}) -\varepsilon^p_{n+1}).\]

\[A_{n+1} = \dfrac{\partial \psi^p}{\partial \alpha}(\alpha_{n+1}).\]

相補性方程式 (5) は， [HA-WO2009] のように， \(r > 0\) に対し，適切に選択された相補性関数を使用して規定されます．

\[\int_{\Omega} (\xi_{n+1} - (\xi_{n+1} + r f(\sigma_{n+1}, A_{n+1}))_+) \lambda dx = 0, \forall \lambda\]

または

\[\int_{\Omega} (f(\sigma_{n+1} + (-f(\sigma_{n+1}, A_{n+1}) - \xi_{n+1}/r)_+ , A_{n+1}) ) \lambda dx = 0, \forall \lambda\]

注記：表記法 \(\Delta \xi_{n+1} = \Delta t \xi_{n+1}\) は書体でよく使われます．ここでは，2つの量の間の差異を保存するためこれを選択しています．これは主に，適応可能な時間ステップの使用のためです．時間ステップが変化するとき，値 \(\xi_n\) は， \(\theta\) -法で， \(\Delta \xi_{n}\) の代わりに \(\xi_n\) を使用することが望ましいです．

平面ひずみ近似¶

平面ひずみ近似は，長手方向の歪み（ \(z\) 軸に沿っていると仮定します）が他の方向の歪みに比べて小さく無視できると考えられる長い円筒形の物体の変形に対応する2次元問題です．これは，平面ひずみテンソルの形をとります

\[\begin{split}\varepsilon(u) = \left(\hspace{-0.5em}\begin{array}{ccc} \varepsilon_{1,1} & \varepsilon_{1,2} & 0 \\ \varepsilon_{1,2} & \varepsilon_{2,2} & 0 \\ 0 & 0 & 0 \end{array}\hspace{-0.5em}\right).\end{split}\]

次の式を示します．

\[\begin{split}\bar{\varepsilon}(u) = \left(\hspace{-0.5em}\begin{array}{cc} \varepsilon_{1,1} & \varepsilon_{1,2} \\ \varepsilon_{1,2} & \varepsilon_{2,2} \end{array}\hspace{-0.5em}\right)\end{split}\]

この式はひずみテンソルの成分を示します．ひずみテンソルの塑性と弾性部分の分解では，

\[\varepsilon^p_{1,3} = \varepsilon^p_{2,3} = \varepsilon^e_{1,3} = \varepsilon^e_{2,3} = 0\]

と次式となります．

\[\varepsilon^e_{3,3} + \varepsilon^p_{3,3} = \varepsilon_{3,3} = 0.\]

塑性モデルへの平面ひずみ近似の適応は，ほとんどの場合，簡単な作業です．等方性線形弾性応答は，

\[\sigma = \lambda \mbox{tr}(\varepsilon(u)) I + 2\mu(\varepsilon(u) - \varepsilon^p),\]

であり，したがって

\[\bar{\sigma} = \lambda \mbox{tr}(\bar{\varepsilon}(u)) \bar{I} + 2\mu(\bar{\varepsilon}(u) -\bar{\varepsilon}^p),\]

応力テンソルの非ゼロの \(\sigma_{3,3}\) 成分は

\[\sigma_{3,3} = \lambda \mbox{tr}(\bar{\varepsilon}(u)) - 2\mu \varepsilon^p_{3,3}\]

均等塑性ひずみが想定される汎用的なケースでは，以下の通りになることに留意してください．

\[\mbox{ tr}(\varepsilon^p) = 0 ~~~~ \rm I\hspace{-0.15em}Rightarrow ~~~ \varepsilon^p_{3,3} = - (\varepsilon^p_{1,1} + \varepsilon^p_{2,2}).\]

平面応力近似¶

平面応力近似は，一般に，薄板の2次元膜変形を記述します．この近似は，面内非ゼロ成分のみを有するように応力テンソルを規定します．

\[\begin{split}\sigma = \left(\hspace{-0.5em}\begin{array}{ccc} \sigma_{1,1} & \sigma_{1,2} & 0 \\ \sigma_{1,2} & \sigma_{2,2} & 0 \\ 0 & 0 & 0 \end{array}\hspace{-0.5em}\right).\end{split}\]

さらに，

\[\begin{split}\bar{\sigma} = \left(\hspace{-0.5em}\begin{array}{cc} \sigma_{1,1} & \sigma_{1,2} \\ \sigma_{1,2} & \sigma_{2,2} \end{array}\hspace{-0.5em}\right)\end{split}\]

のように応力テンソルの面内成分を与えます．弾塑性の場合，一般に，付加的な未知変数 \(\varepsilon^e_{1,3}\), \(\varepsilon^e_{2,3}\), \(\varepsilon^e_{3,3}\) により，応力テンソルの面外成分をゼロにするように2次の塑性流れ則を適用することにより成り立ちます．（例えば， [SO-PE-OW2008] を参照）．

等方性の線形弾性応答 \(\sigma = \lambda \mbox{tr}(\varepsilon^e) + 2\mu\varepsilon^e\) の場合，次のようになります．

\[\begin{split}\varepsilon^e = \left(\hspace{-0.5em}\begin{array}{ccc} \varepsilon^e_{1,1} & \varepsilon^e_{1,2} & 0 \\ \varepsilon^e_{1,2} & \varepsilon^e_{2,2} & 0 \\ 0 & 0 & \varepsilon^e_{3,3} \end{array}\hspace{-0.5em}\right).\end{split}\]

また

\[\varepsilon^e_{3,3} = -\dfrac{\lambda}{\lambda+2\mu}(\varepsilon^e_{1,1} + \varepsilon^e_{2,2})\]

そのため，

(6)¶\[\bar{\sigma} = \lambda^* \mbox{tr}(\bar{\varepsilon}^e) + 2\mu\bar{\varepsilon}^e ~~~\mbox{ with } \lambda^* = \dfrac{2\mu\lambda}{\lambda+2\mu}\]

さらに，

(7)¶\[\|\mbox{Dev}(\sigma)\| = \left(\|\bar{\sigma}\|^2 - \dfrac{1}{3}(\mbox{tr}(\bar{\sigma}))^2\right)^{1/2}.\]

均等塑性ひずみが仮定されている場合，以下が成り立ちます．

\[\mbox{ tr}(\varepsilon^p) = 0 ~~~~ \rm I\hspace{-0.15em}Rightarrow ~~~ \varepsilon^p_{3,3} = - (\varepsilon^p_{1,1} + \varepsilon^p_{2,2}).\]

いくつかの古典的な法則¶

Tresca : \(\rho(\sigma) \le \sigma_y\) ．ここで， \(\rho(\sigma)\) はCauchy応力テンソルのスペクトル半径， \(\sigma_y\) は一軸降伏応力です．(これはいくつかの硬化内部変数に依存する可能性があります）．

Von Mises : \(\|\mbox{Dev}(\sigma)\| \le \sqrt{\frac{2}{3}}\sigma_y\) ．ここで \(\mbox{Dev}(\sigma) = \sigma - \frac{1}{3}\mbox{tr}(\sigma)I\) は \(\sigma\) と \(\|\sigma\| = \sqrt{\sigma:\sigma}\) の偏差部分です．

Von-Mises基準（Prandl-Reussモデル）を伴う完全な等方性弾性熱可塑性¶

内部変数はなく，等方性弾性応答を考慮します．流動化則は次の式になります．

\[\dot{\varepsilon}^p = \gamma \dfrac{\mbox{Dev}(\sigma)}{\|\mbox{Dev}(\sigma)\|}\]

これは式 \(\Psi(\sigma) = f(\sigma) = \|\mbox{Dev}(\sigma)\| - \sqrt{\frac{2}{3}}\sigma_y\) に対応します．

塑性流れ則の積分のための \(\theta\)- 法は次式となります．

\[\varepsilon^p_{n+1} - \varepsilon^p_{n} = (1-\theta)\alpha(\sigma_{n}) \Delta t \xi_n \dfrac{\mbox{Dev}(\sigma_{n})}{\|\mbox{Dev}(\sigma_{n})\|} + \theta\alpha(\sigma_{n+1}) \Delta t \xi_{n+1} \dfrac{\mbox{Dev}(\sigma_{n+1})}{\|\mbox{Dev}(\sigma_{n+1})\|}.\]

因子を選択すると，次のようになります． \(\alpha(\sigma_{n}) = \|\mbox{Dev}(\sigma_{n})\|\) と \(\xi_n = \dfrac{\gamma_n}{\alpha(\sigma_{n})}\) ．ゆえに次式となります．

\[\varepsilon^p_{n+1} - \varepsilon^p_{n} = (1-\theta)\Delta t \xi_n \mbox{Dev}(\sigma_{n}) + \theta \Delta t \xi_{n+1} \mbox{Dev}(\sigma_{n+1}).\]

\(\mbox{Dev}(\sigma_{n+1}) = 2\mu\mbox{Dev}(\varepsilon(u_{n+1})) - 2\mu\varepsilon^p_{n+1}\) であるため次式となります．

\[\tilde{\mathscr E}^p(u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}) = \zeta_n + \left(1-\dfrac{1}{1+2\mu\theta\Delta t\xi_{n+1}}\right)(\mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n),\]

これは， \(u_{n+1}\) （ただし， \(\xi_{n+1}\) とは関係ありません）に関する線形表現です．

さらに， \(\zeta_n\) は次のように定義されます．

\[\zeta_n = \varepsilon^p_n+(1-\theta)\Delta t \xi_n (\mbox{Dev}(\sigma_n)) = \varepsilon^p_n+(1-\theta)\Delta t \xi_n 2\mu \left(\mbox{Dev}(\varepsilon(u_{n}))-\varepsilon^p_n\right).\]

乗数の除去（Return Mapping法の場合）

次式が成り立ちます

\[\|\mbox{Dev}(\sigma_{n+1})\| = 2\mu\|\mbox{Dev}(\varepsilon(u_{n+1})) -\varepsilon^p_{n+1}\| = \dfrac{2\mu}{1+2\mu\theta\Delta t \xi_{n+1}}\|\mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n\|,\]

したがって， \(B = \mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n\) であり，

\[2\mu\|B\| \le \sqrt{\frac{2}{3}}\sigma_y,\]

そして弾性範囲の場合は， \(\xi_{n+1} = 0\) ，または \(\|\mbox{Dev}(\sigma_{n+1})\| = \sqrt{\frac{2}{3}}\) となり，

\[1+2\mu\theta\Delta t \xi_{n+1} = \dfrac{2\mu\|B\|}{\sqrt{\frac{2}{3}}\sigma_y},\]

であり，したがって

\[\varepsilon^p_{n+1} = \zeta_n + \left( 1 - \sqrt{\frac{2}{3}}\dfrac{\sigma_y}{2\mu\|B\|}\right) B.\]

2つの選択肢は次のように表されれます．

\[\varepsilon^p_{n+1} = {\mathscr E}^p(u_{n+1}, \zeta_{n}) = \zeta_n + \left( 1 - \sqrt{\frac{2}{3}}\dfrac{\sigma_y}{2\mu\|B\|}\right)_+ B.\]

乗算器 \(\xi_{n+1}\) ( \(\theta \ne 1\) の \(\theta\) -法である必要があります）は次のように与えられます．

\[\xi_{n+1} = \dfrac{1}{\theta\Delta t}\left(\sqrt{\frac{3}{2}}\dfrac{\|B\|}{\sigma_y} - \dfrac{1}{2\mu}\right)_+.\]

平面ひずみ近似

平面ひずみ近似 \(\bar{\varepsilon}^p\) と \(\bar{\varepsilon}(u_{n+1})\) は3次元歪みテンソルを面内歪みで置き換えるのと同じ表現をしています．

\[\bar{\tilde{\mathscr E}}^p(u_{n+1}, \theta \Delta t \xi_{n+1}, \bar{\zeta}_{n}) = \bar{\zeta}_n + \left(1-\dfrac{1}{1+2\mu\theta\Delta t\xi_{n+1}}\right)(\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - \bar{\zeta}_n),\]

ここで， \(\overline{\mbox{Dev}}(\bar{\varepsilon}) = \bar{\varepsilon} - \dfrac{\mbox{tr}(\bar{\varepsilon})}{3} \bar{I}\) は3次元偏差器の2次元制限です．

さらに，降伏条件については，

\[\|\mbox{Dev}(\sigma)\|^2 = 4\mu^2\left(\|\overline{\mbox{Dev}}\bar{\varepsilon}(u) - \bar{\varepsilon}^p\|^2 + \left(\dfrac{\mbox{tr}(\bar{\varepsilon}(u))}{3} -\mbox{tr}(\bar{\varepsilon}^p) \right)^2\right).\]

そして乗数の除去のために次式とします．

\[\bar{\mathscr E}^p(\bar{u}_{n+1}, \bar{\varepsilon}^p_{n}) = \bar{\zeta}^p_{n} + \left( 1 - \sqrt{\frac{2}{3}}\dfrac{\sigma_y}{2\mu\|B\|}\right)_+ \bar{B}\]

ここで \(\bar{B} = \overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1}))-\bar{\zeta}_{n}\) and \(\|B\|^2 = \|\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - \bar{\zeta}_n\|^2 + \left(\dfrac{\mbox{tr}(\bar{\varepsilon}(u_{n+1}))}{3} -\mbox{tr}(\bar{\zeta}_n) \right)^2\)．

面応力近似

平面応力近似の場合， (6) を使用して，3次元の場合の式から導きます．

\[\bar{\varepsilon}^p_{n+1} = \dfrac{1}{1+2\mu\theta\Delta \xi}\left(\bar{\zeta}_{n} +2\mu\theta\Delta \xi\left(\bar{\varepsilon}(u_{n+1}) - \dfrac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+1})) - \mbox{tr}(\bar{\varepsilon}_{n+1}^p))\bar{I}\right) \right),\]

since \(\mbox{Dev}(\varepsilon(u)) = \varepsilon(u) - \dfrac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u)) - \mbox{tr}(\bar{\varepsilon}^p))\). Of course, this relation still has to be inverted. Denoting \(\alpha = 1+2\mu\theta\Delta \xi\), \(\beta = \dfrac{4\mu^2\theta\Delta \xi}{3\lambda+6\mu}\) and \(C = \bar{\zeta}_{n} +2\mu\theta\Delta \xi\left(\bar{\varepsilon}(u_{n+1}) - \dfrac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+1}))))\bar{I}\right)\) one obtains

\[\bar{\varepsilon}^p_{n+1} = \dfrac{\beta \mbox{tr}(C)}{\alpha(\alpha-2\beta)}\bar{I} + \dfrac{1}{\alpha}C.\]

さらに，降伏条件の場合，式 (7) を使用できます．

線形等方性および移動硬化とVon-Mises基準を用いた等方性弾性率¶

等方性弾性応答と内部変数 \(\alpha : \Omega \rightarrow \rm I\hspace{-0.15em}R\) を以下の条件を満たす累積塑性歪として考察します．

\[\dot{\alpha} = \sqrt{\dfrac{2}{3}}\gamma\]

線形硬化は等方的硬化係数 \(H_i\) であり，次式となります．

\[\psi^p(\alpha) = \frac{1}{2}H_i\alpha^2\]

すなわち \(A = H_i\alpha\) および一軸降伏応力は次のように定義されます．

\[\sigma_y(a) = \sigma_{y0} + A = \sigma_{y0} + H_i\alpha,\]

\(\sigma_{y0}\) は最初の一軸降伏応力です．降伏関数（およびこれは塑性モデル関連しているため，塑性ポテンシャル）は次のように定義されます．

\[\Psi(\sigma, A) = f(\sigma, A) = \|\mbox{Dev}(\sigma - \frac{2}{3}H_k\varepsilon^p)\| - \sqrt{\frac{2}{3}}(\sigma_{y0} + A),\]

ここで， \(H_k\) は動力学的硬化係数です．前節と同じ計算により次式となります．

\[\tilde{\mathscr E}^p(u_{n+1}, \theta\Delta t \xi_{n+1}, \zeta_n) = \zeta_n + \dfrac{1}{2(\mu+H_k/3)}\left(1 - \dfrac{1}{1+2(\mu+H_k/3)\theta\Delta t\xi_{n+1}}\right)(2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n)\]

\[\begin{split}\begin{array}{rcl} \tilde{\mathscr A}(u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}, \eta_n) &=& \eta_n + \sqrt{\dfrac{2}{3}} \theta \Delta t \xi_{n+1}\|\mbox{Dev}(\sigma_{n+1} - \frac{2}{3}H_k\varepsilon^p_{n+1})\| \\ &=& \eta_n + \sqrt{\dfrac{2}{3}} \theta \Delta t \xi_{n+1}\|2\mu\mbox{Dev}(\varepsilon(u_{n+1})) - 2(\mu+H_k/3)\varepsilon^p_{n+1}\| \\ &=& \eta_n + \sqrt{\dfrac{2}{3}} \dfrac{\theta \Delta t \xi_{n+1}}{1+2(\mu+H_k/3)\theta\Delta t\xi_{n+1}}\|2\mu\mbox{Dev}(\varepsilon(u_{n+1})) - 2(\mu+H_k/3)\zeta_{n}\| \\ &=& \eta_n + \sqrt{\dfrac{2}{3}}\dfrac{1}{2(\mu+H_k/3)}\left(1 - \dfrac{1}{1+2(\mu+H_k/3)\theta\Delta t\xi_{n+1}}\right) \|2\mu\mbox{Dev}(\varepsilon(u_{n+1})) - 2(\mu+H_k/3)\zeta_{n}\|\end{array}\end{split}\]

ここで， \(\zeta_n\) と \(\eta_n\) は次のように定義されます．

\[\zeta_n = \varepsilon^p_n+(1-\theta)\Delta t \xi_n (\mbox{Dev}(\sigma_n)-\frac{2}{3}H_k\varepsilon^n_p) = \varepsilon^p_n+(1-\theta)\Delta t \xi_n \left(2\mu\mbox{Dev}(\varepsilon(u_{n}))-2(\mu+H_k/3)\varepsilon^n_p\right),\]

\[\eta_n = \alpha_n+(1-\theta)\sqrt{\dfrac{2}{3}}\Delta t \xi_n \|\mbox{Dev}(\sigma_n)-\frac{2}{3}H_k\varepsilon^n_p\| = \alpha_n+(1-\theta)\sqrt{\dfrac{2}{3}}\Delta t \xi_n \|2\mu\mbox{Dev}(\varepsilon(u_{n}))-2(\mu+H_k/3)\varepsilon^n_p\|.\]

等方性硬化係数は，\(f(\sigma, A)\) 内のみで \(\tilde{\mathscr E}^p(u_{n+1}, \theta \Delta \xi, \varepsilon^p_{n})\) の式には影響しないことに注意してください．

乗数の除去（Return Mapping法の場合）

式 \(\delta = \dfrac{1}{1+2(\mu+H_k/3)\theta\Delta t\xi_{n+1}}\), \(\beta = \dfrac{1-\delta}{2(\mu+H_k/3)}\) と， \(B = 2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n\) によって \(\varepsilon^p_{n+1}\) と \(\alpha_{n+1}\) の式は次のようになります．

(8)¶\[\varepsilon^p_{n+1} = \zeta_n+\beta B, ~~~ \alpha_{n+1} = \eta_n + \sqrt{\dfrac{2}{3}}\beta \|B\|,\]

そして，塑性の拘束は次の通りです．

\[\delta \|B\| \le \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1}).\]

従って，弾性の場合， \(\xi_{n+1} = 0, \delta = 1\) であり次式となります．

\[\|B\| \le \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \eta_n),\]

塑性領域の場合 \(\xi_{n+1} > 0, \delta < 1\) ，\(\delta \|B\| = \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1})\) と \((1-\delta)\) を使用して方程式を解きます．

\[\|B\| - (1-\delta)\|B\| = \sqrt{\dfrac{2}{3}}\left(\sigma_{y0}+H_i \eta_n + \sqrt{\dfrac{2}{3}} \dfrac{H_i}{2(\mu+H_k/3)}(1-\delta)\|B\|\right),\]

それにより，次式となります．

\[1-\delta = \dfrac{2(\mu+H_k/3)}{\|B\|(2\mu+\frac{2}{3}(H_k+H_i))}\left(\|B\|-\sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \eta_n) \right)\]

この2つのケースをまとめると，

\[\beta = \dfrac{1}{\|B\|(2\mu+\frac{2}{3}(H_k+H_i))}\left(\|B\|-\sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \eta_n) \right)_+\]

(8) により \({\mathscr E}^p(u_{n+1}, \zeta_n, \eta_n)\) と \({\mathscr A}(u_{n+1}, \zeta_n, \eta_n)\) が与えられます．乗算器 \(\xi_{n+1}\) は次のように与えられます．

\[\xi_{n+1} = \dfrac{1}{(2(\mu+H_k/3))\theta\Delta t}(\dfrac{1}{\delta}-1) = \dfrac{1}{\theta\Delta t}~\dfrac{\beta}{1-2(\mu+H_k/3)\beta}.\]

平面ひずみ近似

\(\delta = \dfrac{1}{1+2(\mu+H_k/3)\theta\Delta t\xi_{n+1}}\), \(\beta = \dfrac{1-\delta}{2(\mu+H_k/3)}\), \(B = 2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n\) と \(\overline{B} = 2\mu\overline{Dev}(\bar{\varepsilon}(u_{n+1}))-2(\mu+H_k/3)\bar{\zeta}_n\) により面内部では

\[\bar{\tilde{\mathscr E}}^p(u_{n+1}, \theta\Delta t \xi_{n+1}, \bar{\zeta}_n) = \bar{\zeta}_n + \beta \overline{B},\]

\[\tilde{\mathscr A}(u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}, \eta_n) = \eta_n + \sqrt{\dfrac{2}{3}}\beta\|B\|,\]

また

\[\|B\|^2 = \|2\mu\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - 2(\mu+H_k/3)\bar{\zeta}_n\|^2 + \left(2\mu\dfrac{\mbox{tr}(\bar{\varepsilon}(u_{n+1}))}{3} -2(\mu+H_k/3)\mbox{tr}(\bar{\zeta}_n) \right)^2.\]

降伏条件は次の通りです．

\[\delta \|B\| \le \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1}).\]

\(\beta\) は，前の節で述べた式と同じ式 \(\|B\|\) を持ちます． \(\bar{\zeta}_n\) と \(\eta_n\) の表現は，条件に応じて適用する必要があります．

Souza-Auricchio弾塑性則（形状記憶合金用）¶

この流動化則の構成の正当性については，たとえば [GR-ST2015] を参照してください．Von-Mises応力基準と等方性弾性応答，内部変数および特殊タイプの移動硬化は，次のような制約を考慮して検討されます \(\|\varepsilon^p\| \le c_3\) ．塑性ポテンシャルと降伏関数は，

\[\Psi(\sigma) = f(\sigma) = \left\|\mbox{Dev}\left(\sigma - c_1\dfrac{\varepsilon^p}{\|\varepsilon^p\|} - c_2\varepsilon^p - \delta \dfrac{\varepsilon^p}{\|\varepsilon^p\|}\right)\right\| - \sqrt{\frac{2}{3}}\sigma_{y},\]

相補性条件で

\[\delta \ge 0, \|\varepsilon^p\| \le c_3, \delta (\|\varepsilon^p\|-c_3) = 0,\]

ここで， \(c_1, c_2 \mbox{ と } c_3\) はいくつかの物理的なパラメータです． \(\dfrac{\varepsilon^p}{\|\varepsilon^p\|}\) は \(\varepsilon^p = 0\) の単位球全体として処理されることに注意してください．

準備中 ...

弾塑性ブリック要素¶

テストプログラム tests/plasticity.cc, interface/tests/matlab/demo_plasticity.m, interface/tests/matlab/demo_plasticity.py と contrib/test_plasticity 内を参照してください．

汎用的なブリック要素¶

汎用的なブリック要素には2つのバージョンがあります．1つめは塑性乗数が，問題の変数として保持されているときのもので，追加される項は次の形式です．

\[\int_{\Omega} \sigma_{n+1} : \nabla \delta u dx + \int_{\Omega} (\xi_{n+1} - (\xi_{n+1} + r f(\sigma_{n+1}, A_{n+1}))_+) \delta\xi dx = 0,\]

\(r > 0\) がブリック要素によって（弾性係数の点で）選択された特定の値を有し，Return Mapping法が選択されたとき（塑性乗数は単なるデータです）は，次の項を追加するだけです．

\[\int_{\Omega} \sigma_{n+1} : \nabla v dx.\]

ブリック要素をモデル md に追加する関数は次の通りです．

getfem::add_small_strain_elastoplasticity_brick
  (md, mim, lawname, unknowns_type,
   const std::vector<std::string> &varnames,
   const std::vector<std::string> &params, region = size_type(-1));

ここで lawname は実装された塑性法則の名前であり， unknowns_type は塑性乗数が問題もしくは次のイテレーションのために保存されたモデルのデータが(写像手法を返す)未知数である離散化の選択を示します．いずれにしても乗数が格納されることに注意してください． varnames は，塑性法則（少なくとも変位，塑性乗数，塑性ひずみ）に依存する長さの変数とデータ名の集合です． params はパラメータ（少なくとも弾性係数と降伏応力）表現のリストです．これらの表現は，モデルのいくつかのデータ名（または変数名）でもかまいませんが，GWFL汎用弱形式言語の有効なスカラー式（ "1/2", "2+sin(X[0])", "1+Norm(v)" ... など）も可能です． params にオプションで指定される最後の2つのパラメータは，塑性ひずみ積分と時間ステップ dt に使用される theta -法（汎用された台形ルール）の theta パラメータです．省略された場合の theta のデフォルト値は1であり，これは1次整合性がある古典的な後方 Euler 方式に対応します． theta=1/2 は，2次整合であるCrank-Nicolson法（台形則）に対応します． 1/2と1の間の値は有効な値です． dt のデフォルト値は 'timestep' であり，モデルで定義された時間ステップを単に示すだけです（ md.set_time_step(dt) による）．あるいは，任意の式（データ名，定数値など ...）でもかまいません．時間ステップは，1つの反復から次のイテレーションまで変更することができます． region はメッシュ領域です．

利用可能な塑性法則は次の通りです．

"Prandtl Reuss"（または "等方性完全可塑性" ）．硬化なしの等方性弾塑性．変数は変位，塑性乗数，塑性ひずみです．変位は変数でなければならず，前の時間ステップでの変位（通常は "u" と "Previous_u" ）に対応する "Previous_" の接頭辞に対応するデータがあります．塑性乗数には2つのバージョン（通常は "xi" と "Previous_xi" ）が必要です． unknowns_type = DISPLACEMENT_ONLY の場合はデータとして定義され， unknowns_type = DISPLACEMENT_AND_PLASTIC_MULTIPLIER の場合は変数として定義されます．塑性ひずみは，mesh_femまたは（好ましくは）im_data（ mim に対応する）に格納された n x n のデータテンソルフィールドを表す必要があります．データは，第1のパラメータ(Lame係数)，第2のパラメータ（せん断弾性率）および一軸降伏応力です．重要：この規則は3次元表現を実装していることに注意してください． 2次元で使用されている場合，式は単に2次元に転置されます．平面ひずみ近似については，以下を参照してください．
"平面ひずみ Prandtl Reuss" （または "平面ひずみ等方性完全塑性" ）平面ひずみ近似に適合したものと同じ規則． 2次元でのみ使用できます．
"Prandtl Reuss線形硬化" （または "等方性塑性線形硬化" ）．線形等方性および移動硬化による等方性弾性塑性． "Prandtl Reuss" 則と比較した追加変数：累積塑性ひずみ．塑性ひずみと同様に，時間ステップの終了時にのみ記憶されるので，単純なデータが必要である（ im_data 上にあるのが好ましい）．2つの追加パラメータ：移動硬化係数および等方性パラメータ．3次元表現のみです．
" 平面ひずみPrandtl Reuss線形硬化" （または "平面ひずみ等方性可塑性線形硬化" ）．前のものと同じ規則ですが，平面ひずみ近似に適合しています．2次元でのみ使用できます．

重要： small_strain_elastoplasticity_next_iter は，各時間ステップの終わり，次のステップの前（および後処理の前に，これが塑性ひずみおよび塑性乗数の値を設定する）で呼び出されなければならないことを忘れないでください．

さらに，以下の関数は，微小塑性ひずみブリック要素の時間ステップをすすめることを許可します．

getfem::small_strain_elastoplasticity_next_iter
  (md, mim, lawname, unknowns_type,
   const std::vector<std::string> &varnames,
   const std::vector<std::string> &params, region = size_type(-1));

パラメータは add_small_strain_elastoplasticity_brick のパラメータとまったく同じでなければなりません．そのため，この関数の説明を参照してください．基本的には，このブリック要素は塑性ひずみと塑性積演算子を計算し，次のステップのために保存します．さらに，計算された変位を前の時間ステップの変位を格納するデータに（通常は "u" から "Previous_u" に）コピーします．この関数は compute_small_strain_elastoplasticity_Von_Mises を使う前に呼び出さなければなりません．

次の関数は

getfem::compute_small_strain_elastoplasticity_Von_Mises
  (md, mim, lawname, unknowns_type,
   const std::vector<std::string> &varnames,
   const std::vector<std::string> &params,
   const mesh_fem &mf_vm, model_real_plain_vector &VM,
   region = size_type(-1));

mf_vm で近似された微小ひずみ弾塑性項に対するVon Mises応力場を計算し，その結果を`VM`に格納します．他のすべてのパラメータは add_small_strain_elastoplasticity_brick とまったく同じです． small_strain_elastoplasticity_next_iter は，この関数を呼び出す前に呼び出さなければならないことに注意してください．

完全な可塑性のための低水準汎用構築に基づく特定のブリック要素¶

This is a previous version of an elastoplasticity brick which is restricted to isotropic perfect plasticity and is based on the low-level generic assembly. Its specificity which could be interesting for testing is that the flow rule is integrated on finite element nodes (not on Gauss points).

このブリック要素をモデルに追加する関数は以下の通りです．

getfem::add_elastoplasticity_brick
    (md, mim, ACP, varname, previous_varname, datalambda, datamu, datathreshold, datasigma, region);

ここで

varname はブリック要素が追加されたメインの未知変位(u)を表します．
previous_varname は直前の時間ステップでの変位です．
datalambda と datamu はLame係数に対応するデータです．
datathreshold は，調査した材料の塑性閾値を表します．
datasigma は，材料によって支持される応力制約値を表します．使用されるNewton法に必要な時間スキームの2つの反復で構成される必要があります． datasigma が定義されている有限要素法は varname の導関数を表すことができることに注意してください．
ACP は，使用される投影のタイプに対応します．これは abstract_constraints_projection 型を持ち，現時点では，Von Misesのものに対応する VM_projection しか存在しません．

注意: datalambda, datamu と datathreshold は定数でも，同じ有限要素法で記述してもかまいません．

この関数は，Newton法を使用して解かれる接ベクトルと右辺ベクトルを構築します．

さらに，次の関数は

getfem::elastoplasticity_next_iter
    (md, mim, varname, previous_varname, ACP, datalambda, datamu, datathreshold, datasigma);

載荷または除荷後に材料がサポートしている新しい応力制約値を計算し，（一度解を求めたら）次のように変数 varname および datasigma をアップロードします．

\[u^{n+1} \rm I\hspace{-0.15em}Rightarrow u^n \ \ \ \ \ \textrm{ and } \ \ \ \ \ \sigma^{n+1} \rm I\hspace{-0.15em}Rightarrow \sigma^n\]

そして， \(u^n\) と \(\sigma^n\) は計算された新しい値を含み，プロセスを再起動することができます．

次の関数は

getfem::compute_elastoplasticity_Von_Mises_or_Tresca
    (md, datasigma, mf_vm, VM, tresca=false);

datasigma に格納されているVon Mises（または tresca = true の場合Tresca）基準の応力テンソルのを計算します．応力は mesh_fem mf_vm で評価され，ベクトル VM に格納されます．もちろん，この関数は以前の関数 elastoplasticity_next_iter がこの関数の前に呼び出された場合にのみ使用できます．

次の関数は

getfem::compute_plastic_part
    (md, mim, mf_pl, varname, previous_varname, ACP, datalambda, datamu, datathreshold, datasigma, Plast);

載荷と除荷後に出現する材料の可塑性部分を Plast に mf_pl で計算します．

datasigma は，新しい応力制約値を含むベクトル，すなわち，材料の載荷または除荷後のベクトルでなければならないことに留意してください．

微小ひずみの可塑性¶

理論的背景¶

微小ひずみテンソルの加算分解¶

内部変数，自由エネルギーポテンシャル，弾性法則¶

塑性ポテンシャル，降伏関数，塑性流れ則¶

初期境界値問題¶

流れ則の積分¶

平面ひずみ近似¶

平面応力近似¶

いくつかの古典的な法則¶

Von-Mises基準（Prandl-Reussモデル）を伴う完全な等方性弾性熱可塑性¶

線形等方性および移動硬化とVon-Mises基準を用いた等方性弾性率¶

Souza-Auricchio弾塑性則（形状記憶合金用）¶

弾塑性ブリック要素¶

汎用的なブリック要素¶

完全な可塑性のための低水準汎用構築に基づく特定のブリック要素¶

GetFEM

ナビゲーション

Related Topics