疲労き裂伝播における遭遇荷重の履歴影響について

角 洋一 大川鉄平

１．緒言

船体構造設計における疲労強度評価はS-N線図に基づいた手法が一般的であり、S-N線で はき裂が構造部材の板厚を貫通するまでを疲労寿命として考える。しかしながら、現実に は作用荷重や工作精度のバラツキ等の多くの不確定要素が存在するため、板厚貫通までき 裂が成長する可能性は十分考えられ、実際に検査において発見されるき裂も板厚を貫通し ているものが多い¹⁾。一方、それらの板厚貫通き裂は船体構造が有する構造不静定性や圧縮 残留応力等の影響によって比較的緩やかな速度で進展することもあり、その伝播経路次第 では構造に特に害を及ぼさないことも考えられる。そこで、使用中の構造で疲労き裂が発 見された場合の対応として、その伝播経路も含めた疲労き裂伝播挙動のシミュレーション によりき裂を診断し、今後の補修計画を立てることが重要といえる。そのような観点から 著者らのグループは船舶等の溶接板骨構造における疲労き裂の伝播寿命、経路推定のため の有限要素法の繰り返し計算による自動き裂伝播解析システム(以後CP-Systemと呼ぶ)を 開発しており、特にき裂伝播経路については高精度で予測できることが確認されている^2)-10)。 著者らの次の目標は疲労き裂伝播寿命推定の高精度化であるが、現状の CP-System では基 本的に一定振幅荷重のみへの対応であり、簡易的に応力比の影響が考慮できる加藤の式 ¹¹⁾ によってき裂伝播寿命を計算している。一方、船体構造に作用する荷重は、主に波浪に起 因して生じる不規則な変動振幅荷重である。変動振幅荷重を受ける構造の疲労き裂伝播挙 動は、一定振幅荷重の場合に比べ明らかに複雑な挙動を呈し、荷重順序が重要な因子とな ることがこれまでの研究により明らかにされている ¹²⁾。これは過去の繰り返し荷重によっ てき裂先端に形成された塑性域中をき裂が進展することで、き裂縁に残留引張変形層が形 成されることが主要因であると考えられている。そのため、疲労き裂伝播寿命の精度良い 推定のためには、残留引張変形層厚さ及びき裂面の接触を適切に考慮できる手法を確立す ることが重要と考えられる。

一方、豊貞は、疲労き裂進展はき裂先端部で塑性仕事がなされる間に起きると考え、疲労 き裂伝播速度を律するパラメータとして ΔKRP(き裂先端に引張塑性域が形成される区間に 対応した応力拡大係数範囲)が適切であることを提唱した ¹³⁾。更に、豊貞は Newman の

Dugdaleモデルの応用によるき裂開閉口モデル¹⁴⁾を任意の応力分布に対して拡張し、この

モデルによる ΔKRPの数値的な解析法を確立した ¹⁵⁾。この手法を基に開発したプログラム

｢FLARP｣は変動振幅荷重に対する疲労き裂の伝播寿命を高精度に推定できることが実験的 に確認されている。しかしながら、この手法は単純な形状に対する定式化であるため、複 雑な形状をした実構造に対しては直接適用できない。そこで本研究ではCP-Systemに結合 可能なき裂開閉口モデルとして、有限要素解析によって得られるき裂先端応力場パラメー

タ等を利用した新たなき裂開閉口モデルを考案した。考案した手法をCP-Systemに結合し、

中央板厚貫通き裂の試験片に対する試解析を行い、その結果をFLARPと比較した。また水 圧荷重が作用する船体構造モデルについて、荷重履歴を嵐モデル ¹⁶⁾と仮定してき裂伝播シ ミュレーションを行い、遭遇海象の履歴影響による疲労き裂伝播寿命の変動を示した^17）。

２．半無限き裂に対する疲労き裂開閉口モデル

本研究ではき裂を有する有限体の問題を扱うが、き裂先端の塑性域寸法が物体の代表寸 法(き裂長さ、リガメント長さ等)に比べ十分小である場合、そのき裂先端近傍の挙動は半無 限き裂に対するものとほぼ同様と考えられる。そこで、以後では解くべき問題のき裂を半 無限き裂で置き換えて考える。

（１）半無限き裂における応力拡大係数の変化率を考慮したき裂結合力モデル

疲労き裂開閉口モデルを定式化する前段階として、無限板中における無限な長さを持っ たき裂(半無限き裂)に対するき裂結合力モデルを応力拡大係数の変化率を考慮して定式化 する。Fig.1(a)に示すように、半無限き裂の先端に長さaの引張塑性域が生じている問題を 考える。このとき、物体の材料は弾完全塑性体を仮定し、塑性域内のき裂面に対して垂直 方向の応力はλσYで打ち切られるものとする。ただしλ は塑性拘束係数、σYは降伏応力で ある。Fig.1(a)の問題は、Fig.1(b)に示す塑性域先端までき裂先端を伸ばすことで、仮想き 裂先端でモード I の応力拡大係数_k⁾が作用している問題と、Fig.1(c)の仮想き裂部で結合力 -λσYが作用する問題の重ね合わせとして考えることができる。

Fig.1 Superposition of the strip yield model in an infinite solid: (a) Original problem, (b) A semi-infinite crack having stress intensity factor k), (c) A semi-infinite crack subjected to internal load.

ここで、解くべき問題のき裂先端近傍における弾性応力場を、Fig.2の実き裂先端を原点 とする座標系に対して

), 2 (

) 2 0 , (

), 2 (

) 2 0 , (

), 2 (

) 2 0 , (

II II

I I

I I

x x O x b

x k

x x O x b x k

x x O b T x x k

xy y x

+ +

=

+ +

=

+ + +

=

π π τ

π π σ

π π

σ (1)

と近似する。ここに、kI, kIIはモードI及びIIの応力拡大係数、Tはき裂に平行な一様応力 成分、bI, bIIはき裂先端からの距離の平方根に比例するモードI及びII成分の係数である。

このき裂が塑性域の長さaだけ直進した後の応力拡大係数がk) となるので、文献²⁾を参考に すると

{

I^/2 I

}

^,

I b k a

k

k)= + + (2) となる。ただし_kI は有限境界の影響を表すパラメータである。(2)式を見ると、その右辺の 括弧内はき裂進展長さに対する応力拡大係数の変化率を表していることがわかる。仮想き 裂先端では応力の特異性が存在しないことから Fig.1(b)及び(c)の応力拡大係数の和を零と することができるので、塑性域長さaは

{

^{8 ( 2 )}

}

^,

/ ² _Y² _{I I I}

2

I b k k

k

a=π λσ −π + (3) と得られる。ただし、aは微小値として２次の項は省略して計算している。

Fig.2 Crack tip coordinate system.

（2） き裂開閉口モデルの定式化 最大荷重時の計算

き裂に最大荷重が作用し、それによって生じる塑性域の先端が過去に生じた塑性域先端 よりも前方にある状態を考える(最大荷重時に形成された塑性域が過去に形成された塑性域 内に留まる場合、後述する最小荷重時の計算を全て最大荷重に対応したパラメータで置き 換えた収束計算が必要となる)。このときのき裂結合力モデルは Fig.1 において仮想き裂先

端のk)を最大荷重に対応したk)_maxで置き換えたものになる。最大荷重に対応する kI, bI, kI

をそれぞれkmax, bmax, kmaxとすると塑性域長さaは(3)式から

{

^{8 ( 2 )}

}

^,

/ ² _Y² _{max max max}

2

max b k k

k

a=π λσ −π + (4) となり、_k⁾maxは(2)式から

{

max^/2 max

}

^,

max

max k b k a

k) = + + (5) となる。最大荷重時のき裂開口変位vmax(x)は

. ' ln

2 '

) ( 2 ) 2

( 0

max

max ^{x =} _E^{−xk −} _E^Y

∫

^{a +}_{− −}⁻_x^xd

v ξ

ξ ξ π

λσ π

) (6)

として与えられる。ただし

⎩⎨

⎧

= −

. )

1

' /( ₂

strain plane for E

stress plane for E E

ν

(7)

ここにEはヤング率、νはポアソン比であり、座標原点x=0は仮想き裂先端に移動した。(6) 式右辺の第１項は仮想き裂先端でk)_maxが作用する問題の変位、第２項はき裂結合力が作用 する問題の変位にそれぞれ対応する。ここでFig.3(a)のようにき裂先端塑性域及び残留引張 変形層域に n 個の棒要素を配置する。残留引張変形層域の棒要素は実き裂先端からの距離 がhまでの有限な領域のみに配置する(本報の解析では想定される塑性域に対して十分大き

くとりh=50mmとしている)。ただし、これらの棒要素長さは過去の計算ステップで累積的

に計算されるものであるので、最初のステップでは零である。Fig.3(b)に示すように、最大 荷重時の塑性域内の棒要素には λσY の引張応力が作用しているので、本来の棒要素長さ li(i=1,...,k)は

), ' / 1 /(

)

max(x E

v

l_i= _i +λσ_y (8) となる。塑性域内の棒要素長さを更新した後、最小荷重時の計算へ移行する。

Fig.3 Schematic illustration of the crack opening/closing model; (a) Configuration of bar elements, (b) Contraction of a bar element.

最小荷重時の計算

Fig.4(a)に示すように、最小荷重時ではき裂先端に圧縮塑性域が形成され、き裂面はそれ までに形成された残留引張変形層のために接触し、圧縮応力が生じる。この状態は、仮想 き裂先端においてk)_minが作用する問題(Fig.4(b))と、最小荷重時に仮想き裂部及び実き裂面 に作用する応力を分布力対として作用させた問題(Fig.4(c))の重ね合わせとして表すことが できる。このときの仮想き裂部の長さ a*は、過去の履歴を考慮して現時点で最も前方にあ る引張塑性域の長さである。

最小荷重時に生じる棒要素の応力及び変位は以下のようにして計算する。最小荷重に対 応するkI, bI, kIをそれぞれkmin, bmin, kminとすると、k)_minは(2)式から

{

min^/2 min

}

^*,

min

min k b k a

k) = + + (9) となる。第j要素に働く一定直応力をσjとすると、第i要素の変位v_min(xi)は

, ) , ' (

) ( 2 ) 2 (

1 min

min

∑

=

− −

= ⁿ

j

j i j i

i vx x

E k x x

v σ

π

) (10)

となる。ここにv(xi, xj)はき裂面の区間[bj,b_j+1]に単位応力対が作用した時の第i要素の中心位 置(座標値xi)における変位であり、

, ' ln

) 2 ,

(

∫

−^{− +1} − −

−

= ^j +

j b

b i

i j

i d

x x x E

x

v ξ

ξ ξ π

(11)

で与えられる。一方、実き裂部と仮想き裂部の弾性変形しかしない棒要素では

, ) ' / 1 ( )

min(xi i E li

v = +σ (12) が成り立つので、(10)式と(12)式を等値して書き直すと

, ) , ' ( / ) , ' (

) ( 2 2

1 min

⎭⎬

⎫

⎩⎨

⎧ +

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

−

− −

=

∑

≠= i i

i i n

i j j

j i j i

i vx x

E l l x x v E

k

x σ

σ π

)

(13) となり、これをGauss-Seidel法で解くことで各棒要素に働く応力を得ることができる。た だし、棒要素は弾完全塑性体であり、実き裂内では引張応力を受け持たないことから

・xi<-a*(実き裂内)では

σi>0のときσi=0, σi<-λσYのときσi = -λσY

・-a*<xi<0(仮想き裂内)では

σi >λσYのときσi =λσY, σi <-λσYのときσi = -λσY

と収束計算中に置き換えを行う必要がある。得られた棒要素応力を再び(10)式に代入するこ とで、最小荷重時の棒要素変位vmin(xi)が得られる。最小荷重時に圧縮降伏した棒要素の長 さは、

), ' / 1 /(

)

min(x E

v

l_i = _i −λσ_y (14) として再設定される。

Fig.4 Superposition of the strip yield model at the minimum load: (a) Original problem, (b) A semi-infinite crack having stress intensity factor k)_min, (c) A semi-infinite crack sub- jected to internal load.

ΔKRPの計算

負荷過程において、き裂先端に再び引張塑性域が形成される瞬間の荷重を RPG 荷重

(Re-tensile Plastic zone Generated Load)¹³⁾¹⁵⁾と呼ぶ。k_RPGは、き裂が完全開口してき裂先 端の棒要素(第k要素)の応力がλσYとなった瞬間の実き裂先端に対する応力拡大係数として 定義される(仮想き裂先端ではk)

RPGとなる)。仮想き裂部の要素の応力 σi(i=1,...,k-1)は、実き 裂部の要素の応力を零、実き裂先端の要素の応力をλσYとして(13)式を書き直し

(15) となる。一方、実き裂先端の要素において応力が λσYとなる瞬間の応力拡大係数k)’_RPG (k) は

' , ) ( 2 / 2 ) , ( ) , ' (

1 )

( ¹

1

' ⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −

⎥⎦

⎢ ⎤

⎣

⎡ ⎟+ +

⎠

⎜ ⎞

⎝⎛ +

=

∑

⁻

=σ π

λσ λσ

E x x x v x

x v E

l k

k ^k

k

j

j k j k k Y Y k RPG

)

(16) また、実き裂部の要素においてき裂面変位と棒要素長さが一致する瞬間の応力拡大係数

k)’RPG(i)(i=k+1,...,n)は

) ,..., 1 ( ' ,

) ( 2 / 2 ) , ( )

, ( )

( ¹

1

' i k n

E x x

x v x

x v l

i

k ⁱ

k

j

j i j k

i Y i

RPG = +

⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧ −

⎭⎬

⎫

⎩⎨

⎧ + +

=

∑

⁻

= σ π

) λσ

(17) である。(16)、(17)式から得られたk)’_RPG(i)(i=k,...,n)の中で最大となるものをk)

RPGとすると、

このときき裂が完全開口かつ実き裂先端要素の応力がλσYとなる条件を満たす。 (15)~(17) 式の計算を応力が収束するまで繰り返し、得られたk)_RPGからkRPGは

), / ( _{RPG max}

max

RPG k k k

k ) )

= (18) となり、ΔKRPは

RPG,

max k

k K_RP= −

Δ (19) と計算される。

き裂進展時の計算

本計算モデルでは、各計算ステップでΔKRPを計算した後にき裂を進展させる。各ステッ プでのき裂進展長さΔcは

, ) ( K_RP ^m C N c=Δ ⋅ Δ

Δ (20) とする。ΔNは１計算ステップあたりの荷重の繰り返し数に対応し、荷重履歴の状態に応じ て変化させる。C及びmは実験的に得られる定数であり、一般的な鋼材では

C=3.514×10^-11, m=2.692, (unit: ΔKRP in MPa m, da/dN in m/cycle) である¹⁵⁾。

き裂進展時にき裂面に残留する棒要素の長さは、塑性収縮 ¹⁵⁾を考慮して以下のように定 める。

), ) ( '( / 1

1

min i i

Y

i v x

l E κδ

λσ ⁻

= − (21)

) 1 ,..., 1 ( , ) , ' ( / ) , ( )

, ' (

) ( 2

2 ¹

1

−

⎭ =

⎬⎫

⎩⎨

⎧ +

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

−

−

− −

=

∑

⁻

≠=

k i x x E v l l x x v x

x E v

k x

i i i i k

i j j

j i j k

i Y RPG i

i λσ σ

σ π

)

ただし

⎪⎩

⎪⎨

⎧

≥

= <

), 1 ) / ( ( 1

) 1 ) / ( ( ) / (

の場合 の場合

n pi e

n pi e n pi e

γ γ α

γ γ α γ γ

κ α (22)

である。δiはき裂面の接触応力が全て解放されたと仮定した場合の塑性収縮量、γeは現在最 も前方にある塑性域長さ、γpiは現荷重サイクルでの引張塑性域長さ、α, n は定数である。

き裂進展部の棒要素長を設定した後、最大荷重時の計算に戻る。以上の計算をき裂がある 設定した長さに達するまで繰り返す。

３．き裂開閉口モデルの CP-System への結合

CP-System では板厚貫通き裂を局所的に２次元平面のき裂進展問題とみなし、逐次有限

要素解析によりその伝播挙動をシミュレートする。またスーパーエレメント ¹⁸⁾を利用する ことで大規模構造物への適用が可能であり、複数き裂の同時進展に対応できる等の特徴が ある。CP-Systemに前節で定式化したき裂開閉口モデルを結合し、自動プログラムを作成 した。シミュレーションの手順を以下に示す(Fig.5参照)。

(1)データ入力：GUIを用いた専用プリプロセッサーによりき裂伝播領域、初期き裂形状、

材料定数、溶接残留応力等を入力する。荷重履歴は外部ファイルから入力する。き裂が伝 播しない周辺領域はスーパーエレメントとして取り込む。このとき負荷荷重の情報はスー パーエレメントに含まれることになる。

(2)メッシュ分割：き裂伝播領域において改良Paving法により四辺形有限要素を自動生成す

る。

(3)有限要素解析：いくつかの解くべき問題について有限要素解析を行い、き裂先端要素の 応力を算出する。

(4)き裂先端応力場パラメータ等の算出：有限要素解と解析解の重ね合わせ法によりき裂先 端応力場パラメータkI, kII, T, bI, bII及び有限境界影響を表すパラメータkI, kIIを算出する。

(5)き裂開閉口シミュレーション：き裂開閉口シミュレーションによって前ステップから現 在までのき裂伝播寿命を求める。き裂開閉口シミュレーションに必要となるき裂先端応力 場パラメータ等は前ステップと現ステップにおいて有限要素解析で得た値を線形内挿して 用いる。

(6)状態のチェック：KIIの KIに対する割合、K_Iの前ステップに対する変化率、複数き裂間

のき裂伝播寿命の誤差をチェックし、許容値を超えた場合はき裂進展長を再設定して(2)に 戻る。

(7)き裂伝播経路の推定：き裂伝播経路を第一次摂動法と局所対称性規準の組み合わせ ²⁾に

よって求める。

(8)き裂進展経路に沿ってき裂を伸ばし、(2)に戻る。

ドキュメント内 <4D F736F F D CA48B8688CF88F589EF95F18D908F C8E86312C322C368FCD E31816A> (ページ 72-87)