Multiscale Finite Element Analysis of Yield-point Phenomenon in Ferrite–Pearlite Duplex Steels

Abstract

Yield-point phenomena in Ferrite–Pearlite duplex steels were investigated using multi-scale computational simulations. In this multi-scale simulations, stress–strain relationship of Ferrite phase was characterized by an elastoplastic constitutive model considering yield-drop behavior and its material constants were determined by minimizing residual error between a computational simulation and experiment of tensile test, where yield-point phenomenon in a tensile test of Ferrite steel was reproduced. Using the determined material response of Ferrite phase, finite element analyses of Ferrite–Pearlite duplex microstructure were performed to examine its macroscopic material response and its microscopic deformation mechanism. Besides, finite element analyses of tensile test based on the numerical results of microscopic finite element analysis were conducted to reproduce yield-point phenomena in Ferrite–Pearlite duplex steels.

1. はじめに

降伏点現象は明確な降伏強度と塑性変形開始後の応力降下を示す材料挙動であり，様々な材料で確認される¹⁾。特に，炭素鋼では顕著に確認され，炭素などの侵入型固溶元素によって転位運動が阻害されることがメカニズムとして考えられている^2,3,4)。一般に，引張試験では試験片の並行部全体が均一な単軸応力状態となることを想定しているが，降伏点現象のように応力降下を伴う場合，試験片内に局所変形が発生し，伝播することで応力–ひずみ関係において特有の降伏棚を形成する⁵⁾。この降伏棚は試験片内の不均一変形に起因しているので，この引張試験で評価された応力–ひずみ関係は試験片形状やサイズの影響を含み，真の材料応答と言えない。そこで，有限要素法のような離散化数値解法を用いることで，応力降下による局所変形およびその際の降伏棚を含む応力–ひずみ関係を再現することができ，真の材料特性を検討できる^6,7,8,9)。

炭素鋼であっても，パーライト鋼のように炭化物が十分に析出した鉄鋼では，降伏点現象は発現しない。Setoyamaら^10,11)はフェライト–パーライト複合組織鋼において高強度組織であるパーライト組織の連結性に注目して，実験と数値シミュレーションを比較した。ここで，フェライト–パーライト複合組織鋼の引張試験においても，降伏点現象が発現していることが確認されるが，フェライト相の応力–ひずみ関係にはフェライト単相鋼の応力–ひずみ関係から応力降下を除いた応答を使用し，無視されている。また，Watanabeら¹²⁾はセメンタイトの第一原理計算からマルチスケールにフェライト–パーライト複合組織の力学挙動を評価した。ここでも，フェライト相の応力–ひずみ関係には応力降下を除いて定式化された転位密度に基づく結晶塑性構成モデル¹³⁾が用いられた。複合組織において，フェライト相の応力降下が発生した場合，引張試験と同様に局所変形の発生・伝播が現れるはずであるが，これまでの数値シミュレーションに関する研究ではほとんど議論されていない。降伏点現象は一種の不安定現象であり，計算安定性を損なうため，材料不均一性を考慮した数値解析のような複雑な数値シミュレーションでは避けられてきた。

Watanabe and Iwata¹⁴⁾は応力降下を表現する降伏点現象の弾塑性構成モデルを提案した。この構成モデルでは一般的な応力–塑性ひずみ関係を記述する関数に対して，降伏開始直後の応力降下を表現する漸近関数を掛け合わせた数理構造を持ち，既存の研究では議論されていない応力降下の急峻さを議論できる。このような応力降下を伴う構成モデルを有限要素解析へ適用した場合，局所変形を伴う不安定現象が再現されるが，安定的に収束解を得るために配慮すれば，静的陰解法を用いて数値シミュレーションを実施することは可能である。

以上を踏まえて，本研究では，Watanabe and Iwata¹⁴⁾が提案した弾塑性構成モデルにおいて，計算安定性を向上させるための簡略化を施した上で，降伏点現象を含む引張試験の数値シミュレーションを実施し，実験結果と比較することでフェライト鋼の真の応力–ひずみ関係を抽出する。そのフェライト相の応力–ひずみ関係を用いて，フェライト–パーライト複合組織の変形シミュレーションを実施し，金属材料微視組織の変形状態とそれに対応するマクロ応力–ひずみ関係を評価し，その変形メカニズムを議論する。さらに，金属材料微視組織の数値シミュレーションで得られたマクロ応答を用いて引張試験の数値シミュレーションを実行するマルチスケール解析を行う。

2. 降伏点応力降下を表現する弾塑性構成モデル

降伏点現象を有限要素解析で再現するために，フェライト相において塑性変形開始後に応力降下が生じる弾塑性構成モデルを適用する必要がある。本研究では，Watanabe and Iwata¹⁴⁾が提案した構成モデルを簡略化したモデルを採用する。応力積分アルゴリズムや弾塑性接線係数の導出はWatanabeら¹⁵⁾の定式化に従う。降伏点現象では局所変形とともに急激な塑性ひずみの変化を生じるため，ひずみ速度依存性の影響もあると考えられるが，本研究では議論を簡略化するために，ひずみ速度非依存型塑性モデルを採用する。有限ひずみ弾塑性論に基づき弾塑性乗算分解を採用する。

  

$$F=F^{\text{e}}{F}^{\text{p}}$$(1)

ここで，Fは全変形勾配，F^e，F^pはそれぞれ変形勾配の弾性および塑性部分である。式（1）によって，基準配置，中間配置，現配置という三つの基底が定義される。弾性構成モデルとして次式のSt. Venant弾性構成モデルを採用する。

  

$$\widehat{S}={ℂ}^{\text{e}}\frac{1}{2}\left(C^{\text{e}}-1\right)$$(2)

ここで，$\widehat{S}$は中間配置における第2Piola-Kirchhoff応力，${ℂ}^{\text{e}}$は中間配置における弾性接線係数，C^e=F^eTF^eは弾性右Cauchy–Green変形テンソルである。ひずみ速度非依存型塑性構成モデルにおいて，降伏関数は次式のように定義される。

  

$$\varphi =\Sigma -\Gamma \le 0$$(3)

ここで，Σ, Γはそれぞれ応力ノルムと降伏応力である。本研究では，応力ノルムをvon Mises応力に対応する次式で定義する。

  

$$\Sigma :=\sqrt{\frac{3}{2}\text{dev}\left[\widehat{T}\right]:\text{dev}{\left[\widehat{T}\right]}^{\text{T}}}$$(4)

ここで，$\widehat{T}=C^e\widehat{S}$は中間配置におけるMandel応力である。

Watanabe and Iwata¹⁴⁾の提案した現象論的構成モデルでは，降伏応力Γを次式で定義する。

  

$$\Gamma :=q\left(1+r\right)$$(5)

ここで，qは塑性変形履歴の内部変数ξの関数であり，塑性硬化を含む塑性負荷状態の応力–ひずみ関係の概形を定義する。rは降伏点現象を表現するために導入した内部変数ζの関数である。本研究では，次式のVoce型非線形硬化則¹⁶⁾を用いてqを次式で定義する。

  

$$q:={\tau }_{\text{Y}}+\left({\tau }_{\text{sat}}-{\tau }_{\text{Y}}\right)\left(1-\exp \left[-h\xi \right]\right)+H\xi$$(6)

ここで，τ_Yはマスターカーブの初期降伏応力，τ_satは塑性履歴変数が十分に大きくなった際の収束降伏応力，hは非線形加工硬化の感度，Hは線形硬化係数である。rはΓの初期値が上降伏応力τ_upとなり，ゼロへ漸近するような関数として次式で定義する。

  

$$r:=\left(\frac{{\tau }_{\text{up}}}{{\tau }_{\text{Y}}}-1\right)\exp \left[-k_0\zeta \right]$$(7)

ここで，k₀は漸近関数の感度，ζは関数の変化を記述する内部変数である。

Watanabe and Iwata¹⁴⁾は熱力学に基づく定式化によって，二つの内部変数の発展方程式を次のように導出した。

  

$$\dot{\xi }=-\gamma \frac{\partial \varphi }{\partial q}=\gamma \left(1+r\right)$$(8)

  

$$\dot{\zeta }=\gamma \frac{\partial \varphi }{\partial r}=\gamma q$$(9)

式（7）から式（8）においてrは降伏直後は正値を持つが，塑性変形発生とともに速やかにゼロへ漸近する。また，式（9）において両辺を初期降伏応力τ_Yで正規化すると下記のように考えることができる。

  

$$\begin{array}{cc}\dot{\xi }\approx \gamma ,& \frac{\dot{\zeta }}{{\tau }_{\text{Y}}}\end{array}\approx \gamma$$(10)

すなわち，両変数の初期値はともにゼロと考えると，ζ ≈ τ_Yξと近似でき，rもξの関数と考えることができる。また，漸近関数の感度をk=τ_Yk₀と定義しなおせば，降伏応力Γは次式のように簡略化できる。

  

$$\Gamma =\left\{{\tau }_{\text{Y}}+\left({\tau }_{\text{sat}}-{\tau }_{\text{Y}}\right)\left(1-\exp \left[-h\xi \right]\right)+H\xi \right\}\left\{1+\left(\frac{{\tau }_{\text{up}}}{{\tau }_{\text{Y}}}-1\right)\exp \left[-k\xi \right]\right\}$$(11)

以上の簡略化によって構成モデルの内部変数をひとつにでき，有限要素解析への実装が容易になるとともに数値解析をロバストに実行できるようになる。

3. フェライト鋼の降伏点現象の数値シミュレーション

降伏点現象では応力降下によって，局所変形の伝播が生じるため，応力–ひずみ関係において引張試験の形状・サイズの影響が現れる。本研究では，実験に対応した降伏点現象を含む引張試験の数値シミュレーションを実行し，フェライト相の応力降下を含む塑性特性を同定する。

3・1　引張試験片の有限要素モデル

Setoyamaら¹⁰⁾の引張試験を基に引張試験片の有限要素モデルを作成した（Fig.1）。静的陰解法を用いて，Fig.1の引張試験片端面に変位を付与する数値シミュレーションを実行する。

Fig. 1.

Finite element model of tensile specimen.

本研究では，弾性等方性を仮定し，弾性定数はWatanabeら¹²⁾より，ヤング率E_F=200 GPa，ポアソン比ν_F=0.296と定義する。また，この引張試験で得られた応力–ひずみ関係の後半部よりマスターカーブとなるτ_Y, τ_sat, h, Hを下記のように設定した。

  

$$\begin{array}{cccc}{\tau }_{\text{Y}}=0.095\text{GPa},& {\tau }_{\text{sat}}=0.305\text{GPa},& h=40.,& H=0.33\text{GPa}\end{array}$$

また，実験値より上降伏応力τ_up=0.354 GPaとした。残りの応力降下に関する材料定数kは次項にて数値シミュレーションと実験結果が合うように同定する。

3・2　数値解析結果と材料定数の同定

降伏点現象では局所変形が試験片全体へ伝播するため，応力–ひずみ関係における降伏棚発生のひずみ範囲が，材料定数だけでなく，試験片のサイズにも依存する。よって，実験に合わせた有限要素モデルに対して，引張試験の数値シミュレーションを実行し，材料定数を同定することで，フェライト鋼の真の力学特性を評価できる。ここで，数値シミュレーションの応力–ひずみ関係は実験と同様に，荷重–変位関係と試験片寸法から公称応力–公称ひずみ関係を求めて，真応力–真ひずみ関係へ変換する。

材料定数kを変えた際の構成モデル¹¹⁾の応答をFig.2（a）に示す。また，同様に引張試験の数値シミュレーションから計算された応力–ひずみ関係を対応する実験結果とともにFig.2（b）に示す。Fig.2（a）より，材料定数kが大きくなるほどに応力降下が急峻となり，降伏直後の材料応答を制御できる。また，この構成モデルを用いた引張試験の数値シミュレーションの結果から評価される応力–ひずみ関係では，明確な降伏棚が表現される。すなわち，引張試験の数値シミュレーションから計算した応力–ひずみ関係Fig.2（b）は構成式（11）とは全く異なる応答となる。

Fig. 2.

Sensistivity of parameter k on stress–strain curves.

引張試験の実験結果に対応する適切な材料定数kを同定するために，実験と構成モデルの応答の残差ノルムを次式で定義し，最小化する。

  

$$\Delta =\frac{1}{{\tau }_{\text{Y}}}{\displaystyle {\int }_{\epsilon }\left|{\sum }_{\exp }^{\epsilon }-{\sum }_{\text{sim}}^{\epsilon }\right|}d\epsilon$$(12)

ここで，Σ^ε_expとΣ^ε_expは軸ひずみεにおける引張試験の実験と数値シミュレーションにおける応力ノルム値（軸応力値に相当）である。本研究では，降伏棚が現れている軸ひずみ範囲ε∈[0.01, 0.05]にて，数値積分によって残差（12）を評価した。評価した材料定数kと残差の関係をFig.3に示す。Fig.3より，k=300と決定した。

Fig. 3.

Residual error in finite element analyses changing parameter k.

Fig.2（b）より，同定した材料定数を用いた数値シミュレーションによって，実験の応力–ひずみ関係を再現できていることがわかる。ひずみ範囲ε∈[0.01, 0.05]における相当塑性ひずみおよび相当応力の分布をFig.4に示す。上降伏点後の応力降下によって，局所変形が発生し，その領域が伝搬している様子が確認できる。本研究では，初期降伏の発生点を制御しておらず，引張試験片中で計算誤差レベルで応力が微妙に高い点にて降伏が発生する。降伏棚が発生しているひずみ範囲では試験片内が不均一な変形となっているが，試験片平行部の全領域が応力降下を経験した塑性変形状態となると，試験片全体が均一変形状態となっていく。

Fig. 4.

Numerical results: Finite element analyses of tensile test.

降伏点現象のように局所変形が現れる問題では，有限要素メッシュ分割の依存性が顕著となる。上記で同定した材料定数において，有限要素分割数を変化させたときの応力–ひずみ関係の変化をFig.5に示す。メッシュが粗い場合，降伏棚後の材料応答が異なる結果となる。メッシュサイズ1. mmと1.5 mmのモデルでは，ほぼ同じ応答が得られているため，本研究では，メッシュサイズ1.5 mmの有限要素モデルを採用した（Fig.1）。

Fig. 5.

Mesh sensitivity of finite element analyses in tensile test.

以上のように，実験と数値シミュレーションの差異が最小となるように構成モデル（11）の材料定数を決めることで，降伏棚を発現させる要因となっている応力降下を含む材料応答を抽出することができる。

4. フェライト–パーライト複合組織のマルチスケール・シミュレーション

前節で同定した応力降下を含むフェライト相の応力–ひずみ関係を用いて，フェライト–パーライト複合組織の有限要素解析を実施し，マクロ応力–ひずみ関係に対するフェライト相の応力降下の影響を議論する。次に，得られたマクロ応力–ひずみ関係を用いて，引張試験の数値シミュレーションを実行し，実験結果と比較する。ここでは，Setoyamaら¹⁰⁾の実験結果の中でパーライト組織の体積分率21.5%および40.3%の試料を解析対象とする。

4・1　複合組織の有限要素モデル

実験に対応したパーライト体積分率となるようにフェライト–パーライト複合組織の有限要素モデルをFig.6 のように作成した。ここでは，対称性を考慮して同形状の48個の微小構造で分割された切頂八面体で充填して全体領域を構成する。そして，微小領域に対してフェライト相とパーライト組織の特性をランダムに与える。Fig.6の全体領域は54個の切頂八面体を含むため，全体領域は2,592個の微小領域に分割される。各微小領域は6 つの六面体一次要素で分割する。Setoyamaら¹⁰⁾の有限要素モデルよりも小さな領域に対してフェライト相とパーライト組織を分散させることができ，想定のパーライト組織の体積分率となるように有限要素モデルを作成できる。Fig.6の有限要素モデルに対して周期境界条件を与え，マクロ変位勾配を制御して，変形解析を実行する。

Fig. 6.

Finite element model of Ferrite–Pearlite duplex microstructure.

フェライト相の力学特性として前節にて評価した応力–ひずみ関係を用いる。また，パーライト組織の力学特性について，弾性特性（ヤング率E_P，ポアソン比ν_P）はWatanabeら¹²⁾においてマルチスケール解析によって評価された次の値を用いる。

  

$$\begin{array}{cc}{E}_{\text{p}}=216\text{GPa,}& {\nu }_{\text{P}}\end{array}=0.288$$(13)

パーライト組織の塑性特性はSetoyamaら¹⁰⁾で用いられた実験データを基とした応力–塑性ひずみ関係を用いる。

4・2　数値解析結果

はじめに，有限要素モデル（Fig.6）の代表体積要素としての妥当性を確認するために，6つのひずみ成分に対応したマクロ微小ひずみモードを付与し，弾性接線係数を評価する。数値解析結果をTable 1にまとめる。どちらの有限要素モデルにおいても等方な弾性接線係数が得られており，有限要素モデル（Fig.6）が十分にランダムな不均一性を有し，代表体積要素として適切であることが示された。

Table 1. Elastic constants estimated with finite element analyses of Ferrite–Pearlite duplex microstructure.

Pearlite volume fraction	Ceiiii	Ceiijj = λ	Ceijij = μ	E	ν
21.5%	269.49	112.38	78.56	203.4	0.2943
40.3%	272.34	112.75	79.80	206.3	0.2928

次に，Y₁, Y₂, Y₃の各軸方向へマクロ単軸応力状態となるようにマクロ変位勾配を制御し，高ひずみまで付与した。数値解析結果としてマクロ相当応力–相当ひずみ関係をFig.7（a）に示す。相当ひずみは次式で定義した。

  

$${\epsilon }^{*}=\sqrt{\frac{2}{3}\text{dev}\left[\epsilon \right]:\text{dev}\left[\epsilon \right]}$$(14)

Fig. 7.

Numerical results: Stress–strain curves of Ferrite–Pearlite duplex steels.

ここで，$\epsilon =\frac{1}{2}\text{ln}\left[F{F}^{\text{T}}\right]$である。弾性応答と同様に各三軸方向のマクロ応力–ひずみ関係はほぼ同一であったため，Y₁軸方向の結果のみを示した。Setoyamaら¹⁰⁾の評価と同じく，降伏強度を含む応力レベルはパーライト組織の体積分率とともに増加する。パーライト組織の体積分率が異なる2つの有限要素モデルの双方のマクロ応力–ひずみ関係において，降伏後に応力降下が確認されるが，Table 2に示すように複合化による上降伏応力の増分はフェライト相とパーライト組織の強度差から考えると相対的に小さい。これはフェライト相の降伏後，塑性変形が局所変形を伴って急速に進展するため，複合組織全体の降伏強度もフェライト相の上降伏応力で決まってしまうためである。すなわち，弾性域ではほぼ均一変形であるため，フェライト相とパーライト組織の弾性ひずみが一致すると仮定し，フェライト相が上降伏点の応力値τ_upとなる際のフェライト–パーライト複合組織の降伏応力τ_Yは，平均場理論を用いて次式で表される。

  

$${\overline{\tau }}_{\text{Y}}=\left(1-f\right){\tau }_{\text{up}}+f{E}_{\text{P}}^{*}{\epsilon }_{\text{up}}^{\text{e}}=\left(1-f+f\frac{E_{\text{P}}^{*}}{E_{\text{F}}^{*}}\right){\tau }_{\text{up}}$$(15)

Table 2. Yield strength calculated with finite element analyses of Ferrite–Pearlite duplex microstructure and estimated with equation (17).

Pearlite volume fraction	Numerical result	Estimation
0. %	0.354 GPa	0.354 GPa
21.5%	0.357 GPa	0.361 GPa
40.3%	0.364 GPa	0.366 GPa

ここで，fはパーライト組織の体積分率，ε^e_upはフェライト相が上降伏応力に至った際の相当ひずみ値である。E*_F，E*_Pはフェライト相およびパーライト組織の相当弾性剛性であり，

  

$$\begin{array}{cc}{E}_{\text{F}}^{*}=\frac{3E_{\text{F}}}{2\left(1+{\nu }_{\text{F}}\right)}=232\text{GPa},& E_{\text{P}}^{*}=\end{array}252\text{GPa}$$(16)

となる。したがって，フェライト–パーライト複合組織の降伏応力は次式で推定できる。

  

$${\tau }_{\text{Y}}=\left(1+0.0867f\right){\tau }_{\text{up}}$$(17)

Table 2にFig.6の有限要素モデルに対応する推定値を示す。また，応力–ひずみ関係における降伏後の応力降下はパーライト組織の体積分率の増加とともに緩やかになり，パーライト組織の体積分率が50%を超える領域では応力降下が確認できなくなる。

フェライト–パーライト複合組織の応力降下時の変形状態を確認するために，マクロ相当ひずみ範囲ε∈[0.003, 0.012]における相当塑性ひずみ分布をFig.8に示す。パーライト体積分率21.5%の有限要素モデルでは，局所変形がフェライト相において伝搬し，材料組織を横断する形で発生した後，領域全体の変形へ移行している。また，パーライト体積分率40.3%の有限要素モデルにおいても，局所変形は生じているものの，パーライト組織に局所変形の伝搬が妨げられている。その結果，パーライト体積分率が高くなるほどにマクロ応力–ひずみ関係における応力降下の影響は小さくなる。

Fig. 8.

Numerical results: Distribution of equivalent plastic strain in Ferrite–Pearlite duplex steels.

数値解析によって得られたフェライト–パーライト複合組織のマクロ応力–ひずみ関係を用いて，Fig.1 の有限要素モデルに対して引張試験の有限要素解析を実行した。引張試験の数値シミュレーションの結果を対応する実験結果とともに応力–ひずみ関係をFig.7（b）に示す。フェライト–パーライト複合組織から評価されたマクロ応力–ひずみ関係を基に降伏点現象が再現されている。本研究では，パーライト組織の体積分率の増加に伴い，材料微視組織においてフェライト相における局所変形を拘束することで，引張試験において降伏点現象が現れにくくなるメカニズムをマルチスケール数値シミュレーションで再現した。

実験の応力–ひずみ関係は数値シミュレーションよりも高い応力値を示しており，定量的に現象を再現できていない。Setoyamaら¹⁰⁾の研究においても，パーライト体積分率が50%未満の領域では，同様の結果となっており，数値シミュレーションで材料特性を定義するために使用したフェライト単相鋼およびパーライト鋼と対象とするフェライト–パーライト複合組織の各相・組織の特性が整合していないと考えられる。式（17）から考えると実験に対応する降伏応力となるためには，フェライト相の上降伏応力τ_upの上昇が必要となる。パーライト体積分率40.3%の複合組織鋼において，参照の降伏応力が得られるためには，0.5 GPaほどに上昇すると推定され，複合組織鋼のフェライト相の力学特性は本研究で使用したフェライト単相鋼の力学特性とは大きく異なることになる。したがって，フェライト–パーライト複合組織のモデリングのためには，微小材料試験^17,18,19)による局所力学特性の評価が必要である。例えば，微小押し込み試験と数値シミュレーションを連携して局所力学特性を評価する手法^20,21)が提案されており，応力降下を伴う応力–ひずみ関係へ拡張することで，上記の課題を解決できると期待される。

5. おわりに

本研究では，フェライト–パーライト複合組織鋼の引張試験における降伏点現象に着目し，試験片および材料微視構造の不均一変形を考慮したマルチスケール数値シミュレーションを実行した。ここでは，引張試験における降伏点現象を数値シミュレーションで再現し，フェライト相の真の応力–ひずみ関係を抽出した。このフェライト相の力学特性を用いて，降伏点現象を示すフェライト–パーライト複合組織の変形機構を材料微視組織の数値シミュレーションを用いて議論するとともに，複合組織の数値解析から得られたマクロ応力–ひずみ関係を用いて引張試験における降伏点現象を再現した。

降伏点現象のように局所変形を伴う不安定現象は数値シミュレーションで扱うことは容易ではないが，数値解析手法の発達によって可能となってきた。本研究で扱った実験と数値シミュレーションを併用した評価手法は，材料不安定現象のメカニズム解明に寄与できる。

謝辞

本研究の一部は科学研究費補助金21H01220の支援を受けた。（株）豊田中央研究所瀬戸山大吾博士に，実験データを提供いただいた。物質・材料研究機構山本由希氏，谷口祥子博士には数値モデルの作成において多大な技術的支援を頂いた。

文献

• 1)  E.O.Hall: Yield Point Phenomena in Metals and Alloys, Springer, New York, NY, (1970). https://doi.org/10.1007/978-1-4684-1860-6
• 2)  A.H.Cottrell and B.A.Bilby: Proc. Phys. Soc., A, 62(1949), 49. https://doi.org/10.1088/0370-1298/62/1/308
• 3)  W.M.Lomer: J. Mech. Phys. Solids, 1(1952), 64. https://doi.org/10.1016/0022-5096(52)90007-0
• 4)  W.G.Johnston and J.J.Gilman: J. Appl. Phys., 30(1959), 129. https://doi.org/10.1063/1.1735121
• 5)  G.T.Hahn: Acta Metall., 10(1962), 727. https://doi.org/10.1016/0001-6160(62)90041-X
• 6)  J.A.Shaw and S.Kyriakides: Int. J. Plast., 13(1997), 837. https://doi.org/10.1016/S0749-6419(97)00062-4
• 7)  H.Tsukahara and T.Iung: Mater. Sci. Eng. A, 248(1998), 304. https://doi.org/10.1016/S0921-5093(97)00857-5
• 8)  F.Yoshida, Y.Kaneda and S.Yamamoto: Int. J. Plast., 24(2008), 1792. https://doi.org/10.1016/j.ijplas.2008.05.004
• 9)  R.Schwab and V.Ruff: Acta Mater., 61(2013), 1798. https://doi.org/10.1016/j.actamat.2012.12.003
• 10)  D.Setoyama, I.Watanabe and N.Iwata: Tetsu-to-Hagané, 98(2012), 290 (in Japanese). https://doi.org/10.2355/tetsutohagane.98.290
• 11)  I.Watanabe, G.Nakamura, K.Yuge, D.Setoyama and N.Iwata: From Creep Damage Mechanics to Homogenization Methods, Advanced Structured Materials Vol.64, Springer, Cham, (2015). https://doi.org/10.1007/978-3-319-19440-0_24
• 12)  I.Watanabe, D.Setoyama, N.Nagasako, N.Iwata and K.Nakanishi: Int. J. Numer. Methods Eng., 89(2012), 829. https://doi.org/10.1002/nme.3264
• 13)  I.Watanabe, D.Setoyama, N.Iwata and K.Nakanishi: Int. J. Plast., 26(2010), 570. https://doi.org/10.1016/j.ijplas.2009.09.005
• 14)  I.Watanabe and N.Iwata: Tetsu-to-Hagané, 100(2014), 1232 (in Japanese). https://doi.org/10.2355/tetsutohagane.100.1232
• 15)  I.Watanabe, N.Iwata and K.Nakanishi: Trans. Jpn. Soc. Comput. Eng. Sci., 2010(2010), 20100005 (in Japanese). https://doi.org/10.11421/jsces.2010.20100005
• 16)  E.Voce: J. Inst. Met., 74(1948), 537.
• 17)  R.Ando, T.Matsuno, T.Matsuda, N.Yamashita, H.Yokota, K.Goto and I.Watanabe: ISIJ Int., 61(2021), 473. https://doi.org/10.2355/isijinternational.ISIJINT-2020-546
• 18)  Y.Tampa, K.Takagi, S.Ueki, M.Ohta, Y.Mine and K.Takashima: ISIJ Int., 62(2022), 1741. https://doi.org/10.2355/isijinternational.ISIJINT-2022-088
• 19)  V.Jayaram: Annu. Rev. Mater. Res., 52(2022), 473. https://doi.org/10.1146/annurev-matsci-080819-123640
• 20)  K.Goto, I.Watanabe and T.Ohmura: Int. J. Plast., 116(2019), 81. https://doi.org/10.1016/j.ijplas.2018.12.007
• 21)  T.-T.Chen and I.Watanabe: Sci. Technol. Adv. Mater.-Methods, 2(2022), 416. https://doi.org/10.1080/27660400.2022.2129508