完全陰解法による分割背応力モデルに基づく損傷弾塑性解析

日本機械学会論文集（A 編） 原著論文 No.2013-JAR-0608

完全陰解法による分割背応力モデルに基づく損傷弾塑性解析

山王丸 将吾

*1

_{，岡澤 重信}

*2

_{，田中 智行}

*2

Elasto-Plastic Damage Analysis Based on Divided Back Stress Model

by Full-Implicit Algorithm

Shogo SANNOMARU

*1

, Shigenobu OKAZAWA and Satoyuki TANAKA

*1

Mazda Motor Corporation, Engine Performance Development Dept. 3-1 Shinchi, Huchu-cho, Aki-gun, Hiroshima, 730-8670 Japan

This paper describes elasto-plastic damage analysis based on divided back stress by full implicit algorithm. The proposed algorithm guarantees second order convergence for geometrical non-linear iterations using consistent tangent stiffness. We validate the present method in localized necking analysis for checking convergence of residual force and cyclic loading analysis to simulate elastic softening and work-softening due to damage. The residual force is well converged in elasto-plastic damage analysis and the point of crack initiation agreed with that of experiment in necking model. The work-softening and deterioration of macroscopic elastic stiffness are well simulated in cyclic loading analysis. Key Words 1. 緒 言 近年，産業分野において機械部品の設計にCAEが活発に利用されている．その中でも強度設計問題で取り扱う 現象は大負荷を受ける繰返し負荷問題というように複雑化しており，もはや線形弾性解析では取り扱う事は困難 となる場合が多く非線形構造解析が必須となる．このような複雑な強度設計問題の例の１つとして自動車業界で は運転者の急なクラッチミートなどで発生する繰返し大負荷による自動車部品の極低サイクル疲労破壊問題が挙 げられる．このような強度問題に対してより軽量かつ高強度な部品を設計してゆくためには極低サイクル疲労に 伴う材料損傷や材料損傷も考慮した変形挙動を精度良く予測する事が工学的に非常に重要となる． この極低サイクル疲労損傷を机上で予測するためには以下の2点が大きな要となるものと考えられる． ① 弾塑性状態における適切な材料損傷モデルの導入 ② 繰返し入力負荷に対する時系列応力，ひずみといった変形履歴を表現するモデルの導入 まず①に関して材料損傷を記述するモデルは幾つかある．その 1 つのアプローチとして散逸ポテンシャルや降 伏関数に材料の微視的損傷を組込み構成則へ反映させて損傷計算を行う手法が用いられている．例えば汎用ソフ トウェアで代表的に用いられる損傷モデルとして延性破壊の微視的損傷機構であるボイドをボイド体積率として 降伏関数に組込んだ Gurson の降伏関数(1)_{による損傷モデルがあり多くの計算事例と実績} (2) (3) _{がある．しかし繰返} し負荷を受ける状態においてこの Gurson モデルは引張り負荷下で損傷が進展するが，逆に圧縮状態では損傷が回 復してしまうため実験で観測されるような損傷進展を表現する事ができない．よって Gurson モデルは極低サイク ル疲労のような繰返し負荷問題を取り扱うにはあまり適していない． ＊ * 原稿受付 2013 年 8 月 9 日 *1 _{正員，マツダ（株）} エンジン性能開発部 PT 解析グループ（〒730-8670 広島県安芸郡府中町新地 3-1） *2 _{正員，広島大学} 工学研究院 機械システム・応用力学部門（〒739-8527 広島県東広島市鏡山一丁目 4-1） E-mail: [email protected]

: Damage Mechanics, Plasticity, Finite Element Method, Divided Back Stress Model, Return-Mapping , Full-Implicit Algorithm, Consistent Tangent Modular Tensor

さらに文献(4)(5)_{によると図 1 のように弾塑性変形の進行による材料損傷の進展により見かけの巨視的弾性材料} 剛性が低下する事が報告されており，この事実は工学的には大負荷を受けた機械部品の巨視的弾性材料剛性が変 化し，その後の変形挙動に大きな影響を与え機械部品の性能低下に繋がる危険を含んでいる．よってこの巨視的 な弾性材料剛性を机上で予測する事は工学的に重要となる．しかし前述の Gurson モデルではこの巨視的弾性材料 剛性の変化挙動を表現する事ができない．そこで本研究では材料内の微視的損傷とそれによる巨視的弾性材料剛 性の変化を共に表現する事ができる連続体損傷力学(6)_{に着目した．}

(a) Tensile test with un-loading(4) (b) Hysteresis loop(5) Fig. 1 Deterioration of macroscopic elastic stiffness

一方，②について材料が弾塑性変形を伴う繰返し大負荷を受けると Bauschinger 効果(7)_{が現れる事が知られてい} る．そして Bauschinger 効果を表現する手法の 1 つに応力空間における降伏曲面の移動を表現する移動硬化則を用 いる手法が挙げられる．移動硬化モデルとして Prager 則(8)_{や Armstrong-Frederick 則}(9)_{といったモデルが提案され} ている．しかしこれらのモデルでは極低サイクル疲労問題のように数十回の多くの負荷サイクルを伴う応力-ひず み履歴を実測値と上手く追従させて表現する事が難しい．そこで本研究では Chaboche らによって提案された分割 背応力形 Armstrong-Frederick 則(10)_{に着目した．} このモデルでは背応力を複数の分割背応力の組み合わせによって表現するため，より複雑な非線形性を持つ背 応力発展挙動を表現する事ができる．そのためこのモデルはより多くの負荷サイクルに渡って実測値の応力-ひず み履歴を追従可能である事から工学的実用性が高く，現在 ADVENTURE - Cluster や ABAQUS などといった汎用 ソフトウェアにも積極的に取り入れられている． ところで一般に静的な弾塑性問題を有限要素法により解く場合，各物理量の積分手法と接線剛性の計算手法の 違いにより残差力の収束性，計算安定性，計算速度に大きな違いが現れる事が知られている．すなわち散逸ポテ ンシャルから導かれる流れ則の積分計算を行う場合，前進形 Euler 積分による陽解法を用いると残差力の収束性 が極めて悪く十分小さい時間増分を用いなければ計算が不安定となり多大な計算時間を要する．この課題は後退 形 Euler 積分による完全陰解法を用いる事で解決される事が知られている(11)_． しかし後退形 Euler 積分による有限増分アルゴリズムを用いる場合，接線剛性としてこの有限増分アルゴリズ ムと整合の取れた Consistent 接線剛性を用いなければ残差力の良好な収束性が得られない．図 2 はこの事を Simo と Taylor が円孔板引張問題の静的弾塑性計算で調べた検証事例(12)_{である．前進形 Euler 積分による陽解法に対し} て後退形 Euler 積分と Consistent 接線剛性を用いた完全陰解法で解いた場合ついて残差力収束性の違いについて示 されている．この図からも完全陰解法の方が圧倒的に優れた収束性を示している事がわかる．しかし完全陰的ア ルゴリズムは陽的アルゴリズムに対してその導出過程が非常に複雑になる事が難点となる． 以上より数十回に渡る大負荷サイクルの極低サイクル疲労損傷問題をより実用的に解いてゆくためには，損傷 状態，変形状態の表現と計算収束性,計算安定性が重要となる．本研究ではこのために以下の 3 つの条件が必要で あると考える．これら 3 条件を全て満たそうとするとその導出過程は複雑となる． A) 分割背応力形 Armstrong-Frederick 則による移動硬化モデルを用いた加工硬化則 B) A)と連続体損傷力学による損傷を強連成させた損傷弾塑性モデル

有限要素法を考慮した損傷弾塑性解析に関する過去の研究を調べると体系的には主に Lemaitre(13)_や村上(14)_らに よってまとめられており，Consistent 接線剛性について研究されたものや分割背応力形 Armstrong-Frederick 則を用 いた事例も報告(15)~(17)_{されている．しかし上記の A), B), C)の 3 条件を同時に満たした事例が見られない．} そこで本研究ではこれら 3 条件を全て満足する完全陰的アルゴリズムを提案し，簡便な計算モデルを用いてそ の妥当性の検証を行った． Explicit Full-Implicit N u m b er o f it er ati o n 0 5 10 15 20 25 1 2 3 4 5 Step Explicit Full-implicit Iteration Explicit Full-Implicit 10-8 10-6 10-4 10-2 100 102 104 M ax im u m r es id u al f o rce [N ] 0 5 10 15 20 25 Explicit Full-implicit

(a) Hole plane model (b) Number of iteration (c) Maximum residual force in step 4 Fig. 2 Convergence of residual force (12)

2. 分割背応力形 Armstrong-Frederick 則に基づく損傷弾塑性モデル 2・1 分割背応力形 Armstrong-Frederick 則 Chaboche ら(10)によって提案された分割背応力形 Armstrong-Frederick 則は式(1)で示されるように N 組の分割背応 力βkを重ね合わせる事で全体背応力を表現するものである．以下，表記に関してスカラ(0 階テンソル)量は細字， 2 階テンソル量は太字，4 階テンソル量は太字にチルド記号を付けて示す．



  N k k 1 β β （1） 応力空間における降伏曲面の移動量はこの全体背応力を用いて表現される．また各分割背応力は式(2)のように 非線形の背応力挙動を表現する Armstrong-Frederick 則によって表現される． k k ε β β a bkγ p k   _{  （2）} ここで ak , bkは各分割背応力に対する材料定数，ε p _{は塑性ひずみ} , γ は塑性乗数を示す．このように分割背応 力形 Armstrong-Frederick 則は全体背応力を複数の非線形特性の組み合わせにより表現する事ができるため，単一 の背応力を用いた場合の Armstrong-Frederick 則よりも比較的複雑な背応力発展挙動を表現する事が可能となる． 2・2 流れ則 分割背応力形 Armstrong-Frederick 則と連続体損傷力学に基づく材料損傷を考慮した散逸ポテンシャル Ψ (13)(18) は式(3)で示される．

 







                            



1 1 1 1 2 1 2 3 S D N k k k k k y d P D P r Y S D r F a b F R D F F F F Ψ β β β σ I : : ~    （3） ここで F P_{は降伏関数，F} β_{は背応力項，F} D_{は損傷項を示す．σ は Cauchy 応力，D は材料損傷値，r, S は損傷発} 展速度を決める定数である．

式(3)中の Y は損傷発展による内部エネルギ密度の変化量を示すもので損傷エネルギ解放率と呼ばれ式(4)のよ うに応力三軸度関数 Rνと関与する量である．また R は等方硬化内部変数，E はヤング率，ν はポアソン比，p は 静水圧，σeqは相当応力(Mises 応力)を示す．







 



_                    2 2 2 2 1 3 1 3 2 1 2 eq eq p R R D E Y









  （4） さらに_D~e_{を弾性材料剛性テンソル，偏差射影テンソルを} d I~ ，対称射影テンソルをIs ~_{，4 階の恒等テンソルを} I ~_， 2 階の等方テンソルを I とおくとそれぞれ式(5)で表記される．ただし I の成分はクロネッカーのデルタ δ と等し くなる．λ，μ は Lame の定数である．

 





 

                 jl ik ijkl jl il jl ik ijkl s s d s e

















I I I I I I I I I D ~ ~ ~ ~ ~ ~ 2 1 3 1 2 （5） 式(1)のように散逸ポテンシャルに降伏関数以外の項が含まれるため式(6)および式(7)のような非関連流れ則が導 かれる．               Y k k p N D a N R y                N β N ε （6） ここで

















_                                                    S Y k k k d d y d d r Y D Y Ψ D D Y N a b D D Ψ D Ψ N D D D Ψ D y 1 1 1 1 2 3 1 1 1 1 1 2 3 , , : ~ : ~ , , : ~ : ~ , ,      β β σ I β σ I β β σ N β σ I β σ I σ β σ N （7） である．ただし式(6)第 3 式には式(8)で示されるように分割背応力 βkに関する自由エネルギψ βに関して熱力学的 応力すなわち分割背応力βkとその内部変数 Xkとの関係式(18) を用いた．            k k k k k k k a a X X β X X     ： 2 （8） 連続体損傷力学ではマイクロクラックやマイクロボイドなどの微視的な材料損傷によって内力を受け持つ有効 面積の減少が材料の巨視的な弾性材料剛性の低下を引き起こすものと考えており，構成則を式(9)のように表現 する．





e e D D ε σ 1 ~ : （9） 式(6)に示される分割背応力をそのまま使用して応力積分を行うと，全ての分割背応力 βk が未知数となり応力 積分に多大の時間を要してしまう．これを避けるために時刻 n から時刻 n+1 に掛けて塑性ひずみ増分の方向が一 定であると仮定して式(6)第 3 式に対してこの時間増分間での解析解を求めると最終的に式(10)および式(11)が得 られる．ただし添字 n, n+1 はその時刻における量を示す．式(10)および式(11)のように時間増分間における解析 解を用いる事で各分割背応力は時刻 n の既知量となり全体背応力 βn+1のみが背応力の未知数となる．

 

N

 

0 E βn    F   2 3 1 （10） ただし

 

 





                                       





1 1 1 1 1 1 1 1 1 1 1 1 1 2 3 1 n n n d n n n d n n k b k k n k b kn D D D e b a F e k k β σ I β σ I N N N β E ： ： ~ ~        （11） である． 2・3 後退形 Euler 積分 本研究では損傷弾塑性問題を完全陰解法に基づく静的釣合い問題として解く事を目指しているため，式(6)およ び式(10)を後退形 Euler の差分表記で整理すると時刻 n から時刻 n+1 に掛けて各物理量が満たすべき方程式 fi は 式(12)となる．ただし時刻 n+1 において弾塑性負荷中の応力状態は後続降伏曲面上に存在するという条件を式(12) 第 5 式に加えた．

 

_

_

 

 





 

                                                       0 1 2 3 0 2 3 0 1 1 1 1 1 1 5 1 1 1 4 1 3 1 2 1 1 1 1 1 1 1 n n y n n n d n n Y n n n n n n n try e n n n e R D f D N D D f F R R f D D













  β σ I σ 0 N E β f 0 σ N ε σ D f : ~ : ~ , , （12） ただし， etry n 1

ε

は試行弾性ひずみと呼ばれ

ε

ε

ε

_



e





n try e n 1 （13） の様に前回の時刻 n の弾性ひずみに剛性方程式を解いて得られた全ひずみ増分を足し合わせたものである． さらにσ yは等方硬化による相対相当応力を示し，その加工硬化則は式(14)で示すような M 個の指数関数則を組み 合わせた形式を用いる．











     M k R B k y y k e A 1 0 1   （14）

ここでσy0は初期降伏応力, Ak , Bk は材料定数を示す．式(12)は時刻 n から時刻 n+1 の状態を求めるための Return-Mapping 方程式となり σn+1 ,Rn+1 ,βn+1, Dn+1 ,Δγ を未知数とする完全陰的な非線形連立方程式になる． 本研究では式(10)で示される非線形連立方程式を Newton-Raphson 法により求解する．この時，解の逐次修正量 を dσn+1のように d を付けて表現すると式(15)を解く事で修正解を得る事ができる．

 

 

 

 

 

 

                                                                                                                                          5 4 3 2 1 1 1 1 5 1 5 1 5 1 5 1 5 4 1 4 1 4 1 4 1 4 3 1 3 1 3 1 3 1 3 2 1 2 1 2 1 2 1 2 1 1 1 1 1 1 1 1 1 f f f d dD d dR d f D f f R f f f D f f R f f D R f D f f R f f D R n n n n n n n n n n n n n n n n n n n n n n n n f f β σ β σ β σ f f β f f σ f β σ f f β f f σ f       （15） この係数行列である Jacobi マトリックスの各成分を具体的に計算すると最終的に式(16)が導かれる．

 

 

 

 

 

 

 

 

 

_

_

_                                                                                                                                                                                              0 1 2 3 2 3 1 1 2 3 1 0 2 3 2 3 2 3 1 0 1 1 2 3 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 1 5 1 5 1 5 1 5 4 1 4 1 4 1 4 1 4 3 1 3 1 3 1 3 1 3 2 1 2 1 2 1 2 1 2 1 1 1 1 1 1 1 1 1 n n n n Y n Y n Y n n n n n n n e n n n n e n n n n n n n n n n n n n n n n n n n n D R H D N D N N F F F F D D D D f D f f R f f f D f f R f f D R f D f f R f f D R σ N N N 0 σ N E σ σ N β N I 0 σ N 0 0 N N ε β N 0 σ N D β σ β σ f f β f f σ f β σ f f β f f σ f : ~ ~                        ： （16） ただし各微分項は式(17)となる．

















 









S n n n Y e n S n S n n Y n n n n n n n n d n n n r Y D S D N r D Y S N D D D q D q D                                                         1 2 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 2 1 1 2 3 1 2 3 1 1 2 3 1 1 ε σ N σ σ N N N N I β N N N I σ N          : ~ ~

 





 





 

 

1 1 1 1 1 1 1 1 1 1 1 2 3                           





n n n y n n n n d n N k b k N k b k kn dR R d R H D q e a F e b k k      β σ I β E : ~                （17） さらに Updated-Lagrange 形の有限変形モデルに対応できるよう応力積分計算の前に Cauchy 応力と背応力に対し て式(18)で示される共回転変換(11)_{を行った．} T L σ L （18） ただし L は増分内でのスピンが一定であると仮定した時の物質点の回転を示す直交テンソルである． 2・4 Consistent 接線剛性 有限要素法における剛性方程式を Newton-Raphson 法で解くに当たり節点残差力の 2 次収束性を得るには式(12) の後退差分アルゴリズムと整合した弾塑性材料剛性を用いなければならない．これは Consistent 接線剛性と呼ば れ式(19)で定義される． try e n n n n ep

d

d

d

d

1 1 1 1    

_



ε

σ

ε

σ

D

~

（19） ここで式(19)の右辺の微分を求めるために式(12)の全微分を取って整理すると式（20）が得られる．このように式 (20)の係数行列は式(16)と同じ成分になる．

 

 

 

 

 

 

                                                                                                                                              0 0 0 1 1 1 1 5 1 5 1 5 1 5 1 5 4 1 4 1 4 1 4 1 4 3 1 3 1 3 1 3 1 3 2 1 2 1 2 1 2 1 2 1 1 1 1 1 1 1 1 1 0 ε β σ β σ β σ f f β f f σ f β σ f f β f f σ f try e n n n n n n n n n n n n n n n n n n n n n n n n n d d dD d dR d f D f f R f f f D f f R f f D R f D f f R f f D R       （20） 式(16)の左辺の係数行列の逆行列の成分を Cijとおくと式(20)から

 

                                                       0 0 0 1 55 54 53 52 51 45 44 43 42 41 35 34 33 32 31 25 24 23 22 21 15 14 13 12 11 1 1 1 0 ε C C C C C C C C C C C C C C C C β σ etry n n n n n d C C C C C C C C C d dD d dR d ~ ~ ~ ~  （21） となる．さらに 1 行目を展開すると式(22)が得られる． try e n n d dσ 1C11： ε 1 ~ _（22）

これより Consistent 接線剛性は 11 1 1 _C ε σ D~   ~   try e n n ep d d _（23） となる．以上により完全陰的アルゴリズムに基づく損傷弾塑性計算を行う事が可能となる． 本研究では以上の材料モデルのプログラムを主に静的陰解法による大規模高速構造解析を得意とする汎用ソフ トウェア ADVENTURE - Cluster のユーザ関数機能を用いて作成した． 3. 数値計算 以上のアルゴリズムの妥当性を示すための計算として以下の 2 つを検証した． ・検証 ①：局所くびれ問題に対する妥当性と残差力収束性の検証 ・検証 ②：繰返し大負荷を受ける時の応力ひずみ挙動の妥当性の検証（図 1 との比較） 3・1 局所くびれ問題 この検証のために用いた 4 面体 2 次要素の有限要素モデルを図 3 に示す．全体寸法 10 mm × 10 mm × 30 mm に 対し要素寸法 1 mm とした．局所くびれを表現するために上面を X 方向および Y 方向を固定し，Z 方向に引張強 制変位を負荷した．その他の面は単軸引張負荷を表現する対称条件を満たすように拘束し Updated-Lagrange 形の 有限変形を考慮した損傷弾塑性モデルとして計算した．ステンレス鋼のデータ(13) (19)_{を参考にした材料定数を表} １に示す．さらに巨視的き裂発生条件を決める損傷値の臨界値 Dcについて金属材料ではおおよそ 0.5 以内となる 事が知られており(13) ，本計算では Dc = 0.5 を用い全積分点のうち最大の損傷値が Dcに達した時点で計算を終了さ せた．また損傷モデルに対する比較用のモデルとして損傷定数以外の材料定数は表 1 と同じ値を用いた非損傷モ デルについても計算した．

Table 1 Material constants

Young's Modulus E[MPa] 160000.0

Poisson's Ratio ν 0.34

Initial Yield Stress σy0[MPa] 250.0

A1[MPa] 350.0 B1 1.0 A2[MPa] 160.0 B2 4.0 A3[MPa] 260.0 B3 1.0 a1[MPa] 19000.0 b1 4000.0 a2[MPa] 19000.0 b2 5000.0 a3[MPa] 16000.0 b3 200.0 S 0.500 r[MPa] 0.350 Isotropic Hardening Constant Kinematic Hardening Constant Damage Constant Fig. 3 Necking model

図 4 に非損傷モデルと損傷モデルの変形に伴う相当応力(Mises 応力)の分布を示す．ただし引張強制変位量を変

形前の長さで割った値を公称ひずみεnとして定義しεn= 1%，5%，9.5%の状態を示す．図 4(a)から非損傷モデル

では公称ひずみが増加しても局所的な変形は見られず均一に伸びており応力は最小断面部でほぼ均一となってい る．一方，図 4(b)から損傷モデルでは下面で局所くびれ変形が発生している様子が見られ，局所くびれ部の断面 内では応力が低下している事が確認される．

図 5 は横軸に公称ひずみ，縦軸にモデルの下面の断面積減少率を示し局所くびれ状態を調べたものである．こ の図から非損傷モデルでは断面がほぼ一定の割合で減少している事がわかる．一方，損傷モデルでは公称ひずみ 3%程度までは非損傷モデルと同等の断面積減少の振舞いを示すが，それ以降は加速的な局所くびれ変形により 一気に断面積が減少している．さらに図 6 は横軸に公称ひずみ，縦軸に引張面の荷重合計値を変形前の断面積で 割った公称応力を用いた公称応力−公称ひずみ線図を示す．この図からも損傷モデルでは局所くびれによる応力 積載能力の低下挙動が確認できる． 550 490 430 370 310 250 [MPa] 550 490 430 370 310 250 [MPa] εn=1% εn=5% εn= 9.5% εn=1% εn=5% εn=9.5%

(a) Non damage model (b) Damage model Fig. 4 Equivalent stress( Mises stress )

0 2 4 6 8 10 12 14 0.00 0.02 0.04 0.06 0.08 0.10 R educ ti on ra ti o[ ％ ] Nominal strain

Non damage model Damage model 0 100 200 300 400 500 600 0.00 0.02 0.04 0.06 0.08 0.10 N o m in a l st re ss [ M P a ] Nominal strain

Non damage model Damage model

Fig. 5 Reduction of cross section area Fig. 6 Stress-strain curve

以上の様に損傷モデルの有無で局所くびれが発生するか否かが明瞭に別れる結果となった．この事について以 下に考察を行う．図 7 に損傷モデルで計算した損傷値の分布を示す．この図から変形の進行と共に局所くびれ部 で損傷が集中している様子が見られ，断面中央で損傷値が最大となっている様子がわかる．また図 8 に損傷モデ ルで計算した応力三軸度の分布を示す．変形の初期段階（εn=1%）では比較的均一に応力三軸が分布しているが， 変形が進むに連れて応力三軸度の高い領域が断面中央部に集中する． 損傷発展則は式(6)で表現され損傷エネルギ解放率 Y は式(4)のように応力三軸度が高くなるに従って大きくな る事から応力三軸度の高い領域では損傷が発展し易くなる．よって応力三軸度の高い断面中央部では図 7 のよう に損傷が集中する．さらに式(3)のように損傷の増加と共に降伏関数が縮小し応力積載能力の低下による加工軟 化が表現されるため損傷が集中している断面中央部では変形が加速されさらに損傷が進展し易くなる．その結果 として損傷モデルでは局所くびれ変形が発生し，図４(b)のような局所くびれ領域での応力積載能力の低下や図 6 のような公称応力−公称ひずみ線図における応力低下挙動が表現される．

連続体損傷力学では損傷散逸ポテンシャルによる損傷値の発展を解いているが，以上のような変形挙動や損傷 挙動は延性材料を用いた実験ではボイドの発生，成長，合体による材料損傷とそれに伴う巨視的応力積載能力の 低下や局所くびれの発生として観測される．また断面中央部で損傷が最大となる事は単軸負荷下の延性破壊は周 知(20) の様に断面中央部において高い応力三軸度のために初期ボイドの成長による材料損傷が顕著に起こり図 9 の ように局所くびれ発生箇所の断面中央部において初期破壊が起きる事と対応している． 0.5 0.4 0.3 0.2 0.1 0.0 0.50 0.44 0.38 0.32 0.26 0.20 0.50 0.44 0.38 0.32 0.26 0.20 1.0 0.8 0.6 0.4 0.2 0.0 εn=1% εn=5% εn=9.5% εn=1% εn=5% εn=9.5%

Fig. 7 Damage ( damage model ) Fig. 8 Stress triaxiality ( damage model )

(a) SAE8620 (b) Cupper Fig. 9 Crack initiation(20)

次に計算の収束性について述べる．本計算は 1ステップ当たり公称ひずみ増分を0.1%ずつ与えて計算終了まで に 95 ステップ必要とした．図 10 に各計算ステップにおいて非線形剛性方程式を Newton-Raphson 法で解く際に残 差力が収束するまでに要した反復回数を示す．ただし残差力の最大値が 1.0×10-4 N 以下に達した時を収束状態と 定義した． 非損傷材料では概ね 4 回程度の反復回数を必要としたが損傷モデルでは公称ひずみの増加に従って反復回数 が増加し最大 6 回の反復回数となった．これは損傷モデルでは上述のように局所くびれ変形が生じるため１ス テップあたりの局所変形量が大きくなるためであると考えられる．いずれにせよ図 2 で示した陽解法のように数 十回の反復を必要とせず，非常に少ない反復回数で収束解を得る事が出来ている． また図 11 は最終ステップ（εn = 0.095）における反復回数に対する残差力の遷移を示したものであるが，損傷モ デルでも非損傷モデルと同様に優れた 2 次収束性を有している事が確認できる．このように良好な収束性が得ら れる理由は本計算に完全陰的アルゴリズムを用いているためであり，以上から第 2 章で示した完全陰的アルゴリ ズムが妥当であると考える．

0 1 2 3 4 5 6 7 0.00 0.02 0.04 0.06 0.08 0.10 N u m b er o f it er at io n Nominal strain

Non damage model Damage model 1 2 3 4 5 M ax im u m r es idu al f o rce [N ] Iteration

Non damage model Damage model 10-5 10-4 10-3 10-2 10-1 100 101 102

Fig. 10 Number of iteration Fig. 11 Residual force

3・2 繰返し大負荷を受ける問題 計算時間を省くため図 12 のように要素数の少ないモデルを用い，上面に対して Z 方向の両振りの繰返し強制 変位を与えた．またその他の面については単軸負荷の対称条件を満たす様に拘束し損傷弾塑性モデルとして計算 した． 比較用の実測データの選定として一般的な鋼材やアルミ材でも同様の傾向は見られるため差し支えないが，図 1(b)に示したように損傷に伴う変形挙動がより顕著に見られるため事前に腐食させた S500s 鋼のデータ(5)を参考 にした．ヤング率が一般的な鋼に対して低めとなっている理由は腐食によるものと考えられる． またここで参考にした実材料の詳細な損傷定数が不明なため各材料定数を表 2 のように現実的な金属材料の範 囲で便宜的な値(13) _{を入力し定性的な変形挙動を表現できるかのみ検証した．}

Table 2 Material constants

Young's Modulus E[MPa] 65000.0 Poisson's Ratio ν 0.34 Initial Yield Stress σy0[MPa] 520.0

A1[MPa] 10.0 B1 2.0 A2[MPa] 20.0 B2 4.0 A3[MPa] 10.0 B3 1.0 a1[MPa] 19000.0 b1 4000.0 a2[MPa] 19000.0 b2 5000.0 a3[MPa] 16000.0 b3 200.0 S 0.500 r[MPa] 15.000 Dc 0.500 Damage Constant Isotropic Hardening Constant Kinematic Hardening Constant Fig. 12 Cyclic loading model

図 13 にこのモデルに公称ひずみ振幅 4%の両振負荷を 10 c 掛けた時の応力ひずみ線図を非損傷モデルと損傷モ デルについて示す．横軸は強制変位量を元の長さで割った公称ひずみを示し，縦軸は引張面の荷重合計値を変形 前の断面積で割った公称応力を示す．実測結果では図 1(b)のようにサイクル数の増加に伴う明瞭な加工軟化挙動 が見られるが，非損傷モデルでは図 13(a)のようにヒステリシスループが塑性シェイクダウンを起こし実測のよ うな加工軟化挙動が見られない．一方，損傷モデルでは図 13(b)のサイクル数の増加に伴う加工軟化挙動が確認 できる．図 14 はサイクル数に伴う各ループの引張時のピーク応力を示したものである．この図からも損傷モデ

ルでは図 1(b)の実測と同様にサイクル数に伴う応力低下が見られ加工軟化挙動を定性的に表現できている事がわ かる． また図 15 にサイクル数に伴う損傷値の変化を示す．図のように損傷モデルではサイクル数の増加に伴い損傷 値が増加してゆく傾向が見られる．すなわち損傷モデルにおいて加工軟化傾向が表現されるのは，この損傷値の 増加に対して式(3)のように降伏関数の収縮による応力積載能力の低下が要因であると考えられる． さらに実測結果では図 1(b)のようにサイクル数の増加に伴い弾性領域のヤング率が低下している様子が確認で きる．図 16 は本計算で得られたサイクル数の増加に伴うヤング率の変化を示したものである．非損傷モデルで は弾塑性変形の履歴に関わらずヤング率は一定値となっているが，損傷モデルでは実測結果と同様にサイクル数 の増加に対してヤング率が低下している．これは図 15 のようにサイクル数の増加に伴って大きくなる損傷値に 対して式(9)で示されるように巨視的弾性材料剛性の低下が表現されるためである． 以上のように本計算手法を用いる事で繰返し大負荷を掛けた時の応力，ひずみ挙動を定性的に良く表現でき ている事がわかる． -800 -600 -400 -200 0 200 400 600 800 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 N o m in al s tr ess[ M P a] Nominal strain -800 -600 -400 -200 0 200 400 600 800 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 N o m in a l st re ss [M P a] Nominal strain (a) Non damage model (b) Damage model

Fig. 13 Stress-Strain curve(FEM )

0 100 200 300 400 500 600 700 1 2 3 4 5 6 7 8 T en si le p e ak s tr e ss [M P a] Cycle

Non damage model Damage model 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 1 2 3 4 5 6 7 8 D am a g e Cycle

0 10 20 30 40 50 60 70 1 2 3 4 5 6 7 8 Y oung 's m odul us [M P a ] Cycle Non damage model Damage model

Fig. 16 Deterioration of macroscopic elastic stiffness(FEM) 4. 結 言 得られた知見を以下にまとめる． 1) 分割背応力形 Armstrong-Frederick 則に基づく損傷弾塑性モデルを完全陰解法で解くため応力積分アルゴリ ズムおよび Consistent 接線剛性を示し，残差力が 2 次収束性を持つ安定で良好な計算手法を提案した． 2) 局所くびれ問題を有限要素法で計算し，非損傷モデルに対して損傷モデルでは損傷に伴う応力積載能力の 低下により局所くびれ変形が加速され最終的に断面中央部で巨視的き裂が生じるという実測結果と良く 一致する計算結果が得られた． 3) 変位一定の繰返し大変形問題を有限要素法で計算し非損傷モデルに対して損傷モデルでは損傷発展に伴う 加工軟化挙動や巨視的ヤング率の低下挙動といった実測に見られる定性的な振舞いを表現する事ができ た． 謝 辞 本研究を進めるに当たり，ADVENTURE- Clusterについて助言を戴いた株式会社アライドエンジニアリングに 謝意を表す． 文 献

(1) Gurson, L., “Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I - Yield Criteria and Flow Rules for Porous Ductile Media”, Transactions of the ASME, Journal of Engineering Materials and Technology, Vol. 99 (1977), pp. 2-15.

(2) 菊池正紀, 山王丸将吾, “シャリップ破壊も含めた延性破壊過程の研究”, 日本機械学会論文集 A 編, Vol.73, No. 732 (2007), pp. 934-941.

(3) 菊池正紀, 山王丸将吾, “混合モード荷重下における延性破壊機構の研究”, 日本機械学会論文集 A 編, Vol. 74, No.

745 (2008), pp. 1235-1242.

(4) Lemaitre, J., “A Continuous Damage Mechanics Model for Ductile Fracture”, Journal of Engineering

Materials and Technology, Trans. ASME, Vol.107 (1985), pp.83-89.

(5) Apostolopoulos, C. A. and Rodopoulos, C. A., “Inelastic Cyclic Behavior of As-Received and Pre-Corroded S500s Tempcore Steel Reinforcement”, International Journal of Structural Integrity, Vol. 1 (2010), pp.52-62. (6) Lemaitre, J., A Course on Damage Mechanics (1996), Springer.

(7) 八高隆雄, 長谷川正, “バウシンガー効果に関する研究の歴史と現状”, 日本鉄鋼協会論文集 (1984), pp.1551-1558.

(8) Prager, “Recent Developments in the Mathematical Theory of Plasticity”, Journal of Applied Phisics, Vol. 20 (1949), pp.235-241.

(9) Armstrong, P. J., Frederick, C.O., “A mathematical Representation of the Multiaxial Bauschinger Effect”, CEGB Report RD/BN/731, Berkeley Nuclear Laboratories (1966).

(10) Chaboche, J.L., “Time-Independent Constitutive Theories for Cyclic Plasticity”, International Journal of Plasticity, Vol. 2, No. 2 (1986), pp. 149-188.

(11) 久田俊明，野口裕久，非線形有限要素法の基礎と応用 (1996), 丸善株式会社.

(12) Simo, J.C. and Taylor, R.L., “Consistent Tangent Operatiors for Rate-Independent Elastoplasticity”, Computer Methods in Applied Mechanics and Engineering, Vol.49 (1985), pp.101-118.

(13) Lemaitre, J., Desmorat,R. , Engineering Damage Mechanics (2005), Springer. (14) 村上澄男, 連続体損傷力学 (2008), 森北出版株式会社.

(15) Hayakawa, K., Nakamura, T. and Tanaka, S., “Elastic-Plastic Behavior of WC-Co Cemented Carbide Used for Forging Tool Considering Anisotropic Damage and Stress Unilaterality”, International Journal of Damage Mechanics (2010), pp. 421-439.

(16) Sawyer, J.P.G., Wang, C.H. and Jones, R., “An Implicit Algorithm Using Explicit Correctors for the Kinematic Hardening Model with Multiple Back Stress”, International Journal of for Numerical Methods in Engineering, Vol. 50 (2001), pp.2093-2107.

(17) Doghri, I., “Numerical Implementation and Analysis of a Class of Metal Plasticity Models Coupled with Ductile Damage”, International Journal for Numerical Methods in Engineering, Vol. 38, Issue 20 (1995), pp. 3403-3431.

(18) 寺田賢二郎 監訳, 非線形有限要素法 (2012), 森北出版株式会社.

(19) 櫻庭健一郎, “ステンレス鋼における最適疲労設計基準の確立”, 産業技術研究センター研究報告, Vol. 4 (2009). (20) Bluhm, J. I. and Morrissey, R. J., “Fracture in a Tensile Specimen”, Proceedings of the First International