機械の研究

(1)

構造解析：V-XFEM による疲労き裂進展シミュレーション

長嶋利夫（上智大，理研），今井登（理研） 1. はじめに計算固体力学分野においては，構造物内部のき裂やボイドなどの不連続な変位場を有限要素メッシュと独立に定義することができる拡張有限要素法(X-FEM)[1][2]が提案され，破壊力学問題における停留き裂解析，き裂進展解析に用いられるようになってきている[3]-[5] [10]-[13]．X-FEM を用いればき裂進展に伴うメッシュ再分割が不要になるので，き裂進展シミュレーションを従来のFEM よりも効率的に実施できる．一方，理化学研究所で開発されているVCAD[6]-[7]は，現状の三次元 CAD が有する三次元形状情報に加え，物性情報をもつボリュームデータを扱うことができる．著者らはVCAD で用いられるボリュームデータとの入出力が可能な構造解析プログラムV-X3D[8]の開発を行っている．V-X3D においては，ボリュームデータを利用して仮想仕事の原理に基づく支配方程式を離散化するために必要な，体積積分，面積積分，さらには変位境界条件を満足させるための拘束条件方程式の導出が行われる . また，プログラム実装においては VCAD-Framework を利用する. 著者らは構造物の疲労き裂進展解析を効率的に実施するために，V-X3D に X-FEM によるき裂進展解析機能を追加することを検討している．本稿では，そのための基本技術となる X-FEM によるき裂進展解析手法および V-X3D による応力解析の手法の概要を紹介する． 2. X-FEM によるき裂進展解析ここでは，き裂形状の表現にレベルセット法[9]を用いた X-FEM(レベルセット X-FEM)による三次元弾性体内部の平面状のき裂のき裂進展解析手法を示す． 2.1 レベルセット関数によるき裂の表現三次元等方性弾性体内部に平面状のき裂を想定する．図1に示すようにき裂面はx1-x2平面上にあるものとし，き裂の前縁形状をレベルセット関数で近似表示する．すなわち，き裂近傍の領域において次式に示すような符号付き距離関数fを定義する．

))

(

)

(

min

)

(

x

n

x

−

=

∈ T

sign

f

Γ (1) ここに，Γはき裂前縁，n(x)はΓ上の点xにおけるΓに直交する単位ベクトルである． fの絶対値は，任意の点からき裂前縁（x1-x2平面）へ降ろした投影点とき裂前縁との距離である．また，その符号はき裂面に対する位置を表している．すなわち，正であればき裂前縁の外側（リガメント），負であればき裂面である．式(1)で定義されるレベルセット関数fの値を節点で定義したうえで，これらを有限要素毎

(2)

に定義される内挿関数で補間することによって解析領域の任意の点でのレベルセット関数値を近似的に評価することができる．また，き裂前縁形状はレベルセット関数のゼロ等値線で表現することができる． Crack surface r θ Evaluation point Crack front：Γ f f 1

x

2

x

1 x 2 x 3 x 3 x Crack surface r θ Evaluation point Crack front：Γ f f 1

x

2

x

1 x 2 x 3 x 3 x 図 1 レベルセット関数による三次元き裂形状の定義 2.2 内挿関数解析領域は三次元 8 節点六面体一次要素を用いて要素分割されているものとする． X-FEM においては，き裂近傍の任意の点 x における変位場 uh を，次式を用いて近似的に表わす． (2.1) I I I k k I h k I I I I I h

_N

_f

_x

_N

_H

_x

b

x

a

x

u

x

u

J C

)

(

)

(

)

),

(

)

(

)

(

)

(

₃ 4 1 3 8 1

∑

∈ = ∈ =

+

=

γ

(2.2)

)

(

)

(

)

(

8 1 I I I h

_N

_f

f

x

∑

x

=

ここにNIは通常の FEM で用いられる 8 個の節点から構成される有限要素の内挿関数，uI， aJ，bIは節点自由度，C はき裂近傍の変位場を近似する 4 個の基底関数γI（I=1,…,4）の 自由度をエンリッチする節点の集合，J はヘビサイド関数 H(x)の自由度をエンリッチする 節点の集合，また，C∩ J=φである． ここでは次式に示すような二次元き裂問題におけるき裂先端近傍の変位場の漸近解を再構成することができる 4 個の基底関数γ_I（I=1,…,4）を用いる．

θ

γ

θ

γ

θ

γ

θ

γ

)

sin

2 cos(

,

sin

)

2 sin(

),

2 sin(

),

2 cos(

2 3 4 1

=

r

=

r

=

r

=

r

(3) ここに，r，θ は評価点のき裂前縁からの距離と角度であり，レベルセット関数を用いて次 式で表わされる． 2 3 2

_x

f

r

=

+

(4.1)

)

/

arctan(

x

₃

f

=

θ

(4.2)

(3)

2.3 数値積分法 X-FEM においても通常の FEM と同様に仮想仕事の原理式に基づき離散化方程式を導出するので，その過程において，要素単位での数値積分が必要になる．き裂前縁近傍のエンリッチされた節点を含む要素における内挿関数は，多項式より高次の関数を含んでいるため，通常の 8 節点六面体一次要素で用いられるようなガウスの 2 点公式を用いた数値積分では精度的に不十分である．X-FEM 解析においてき裂面が要素内に含まれる場合には，領域を複数の単純な立体に分割したうえで部分領域毎にガウス積分が適用される．ここでの解析では，剛性評価の処理を簡単化するためにエンリッチされた節点を含む要素において，ガウス積分点を増加させることによって剛性評価のための数値積分を行う． 2.4 エネルギー解放率の評価方法 X-FEM 解析で用いられる解析モデルにおいては，一般にき裂前縁形状と要素境界とは完全に一致しないので，き裂先端あるいは前縁でのエネルギー解放率を評価するために，従来の FEM と組み合わせて用いられている仮想き裂閉口法を直接適用することは困難である． X-FEM の解析結果からエネルギー解放率を計算する場合，領域積分表示された定義式を用いて計算する[1]-[5][10]-[13]．すなわち，図 2 に示すように，き裂前縁の任意の評価点においてき裂進展方向が

x~

方向となるように局所直交座標系

~

x

1,

~

x

2, とその点を囲む任意の領域V を設定し，その点におけるエネルギー解放率を次式により計算する． 3

x

dV x q x u x u x u x q x u x u x u x q w x u x u x u G V ⎥ ⎦ ⎤ ∂ ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ⎢ ⎣ ⎡ ∂ ∂ − ∂ ∂ + ∂ ∂ + ∂ ∂ =

_∫

3 1 3 3 1 2 ~ 3 2 ~ 1 1 ~ 3 1 ~ 2 1 3 3 2 ~ 1 2 ~ 2 ~ 1 1 ~ 2 ~ 1 ~ 1 1 3 3 1 ~ 1 2 ~ 2 ~ 1 ~ 1 1 ~ 1 ~ ) ~ ~ ~ ( ~ ) ~ ~ ~ ( ~ ) ~ ~ ~ (

σ

τ

σ

τ

σ

(5.1)

∫

= c L dC G G

α

(5.2) ここにα は評価点で１，領域 V の表面でゼロとなるスカラー関数，Lcはき裂前縁，σijは応力，ujは変位，w はひずみエネルギー密度である． ここでは Sukumar ら[5]の研究に倣い，積分領域V として，エネルギー解放率の評価点を 中心とし主軸方向が

~

x

₁

~

x

2

x

3方向となる直方体領域を設定する．その直方体内部に応力，ひずみ，変位の評価点を設ける．重み関数 q としては，評価点で１をとり，直方体の側面で ゼロとなるような次式で表される関数を用いる．

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−

⎟

⎠

⎞

⎜

⎝

⎛

−

⎟

⎠

⎞

⎜

⎝

⎛

−

=

3 3 2 ~ 2 1 ~ 1 3 2 1

2

1 ~

2

1 ~

2

1 )

,

~

,

~

(

L

x

L

x

L

x

q

(6) ここに， ₃は局所座標系方向に沿った直方体の辺の長さである． 2 ~ 1 ~

,

L

,

L

1

~

_x

2

~

_x

3

x

式(5.1)に示す体積積分は，直方体領域においては多数の積分評価点を用いて評価される．

(4)

1 ~ x 1 L L2 3 L 3 ~ x 2 ~ x 1 x 2 x 3 x Evaluation point Numerical Integration Integration point 1 ~ x 1 L L2 3 L 3 ~ x 2 ~ x 1 x 2 x 3 x Evaluation point 1 ~ x 1 L L2 3 L 3 ~ x 2 ~ x 1 x 2 x 3 x 1 ~ x 1 L L2 3 L 3 ~ x 2 ~ x 1 x 2 x 3 x 1 x 2 x 3 x Evaluation point Numerical Integration Integration point 図 2 領域積分法によるエネルギー解放率の評価 2.5 き裂進展解析手法 X-FEM によるき裂進展解析[10]- [12]の流れを図 3 に示す．まず，初期き裂形状を考慮せずに汎用の FEM 用プリプロセッサーを用いて三次元有限要素モデルを作成する．また，別途平面状の初期き裂の前縁形状を点列とそれを結ぶ直線でモデル化しておく．つぎにプログラム処理により，作成した FEM メッシュとき裂の前縁形状に対してエンリッチ節点を自動的に設定し，X-FEM 解析の入力データを作成する．三次元 X-FEM 解析を実施し，き裂前縁でのエネルギー解放率，あるいは応力拡大係数を評価する．さらにこれらの値から，き裂進展則によりき裂進展量を評価し，き裂形状を更新する．このような処理を繰り返すことによりき裂進展解析が実行可能となる． X-FEM Solver

X-FEM Result Displacement, Stress/Strain Energy release rate

Updating front tip

Crack propagation law 3DFEM Data FEM Modeler Geometry of crack front Delamination Modeler X-FEM Data X-FEM Solver

X-FEM Result Displacement, Stress/Strain Energy release rate X-FEM Solver

X-FEM Result Displacement, Stress/Strain Energy release rate

Updating front tip

Crack propagation law

Updating front tip

Crack propagation law 3DFEM Data FEM Modeler 3DFEM Data FEM Modeler 3DFEM Data FEM Modeler Geometry of crack front Geometry of crack front Delamination Modeler X-FEM Data Delamination Modeler X-FEM Data Delamination Modeler X-FEM Data 図 3 X-FEM によるき裂進展解析の流れ 2.6 き裂進展解析例一様分布引張荷重を受ける半楕円き裂（短半径 20 mm，長半径 100mm）を有する厚肉平板 (厚さ 50 mm，長さ 1000 mm，幅 1000mm，ヤング率 207GPa，ポアソン比 0.3)を対象として弾性き裂進展解析を実施する．解析においては対称性を考慮して 1/4 領域をモデル化する[13]．

(5)

図 4 に解析対象の形状と解析モデルを示す．解析モデルとして構造的に 50x20x50 分割したメッシュ（節点数 22,491,要素数 20,000）を用いる．ここでは繰り返し負荷される 100 MPa の一様引張り荷重によるモード I の疲労き裂進展を想定して，き裂進展則を次式のように設定する．

C

G

m

dN

da

₌

_Δ

（7）ここにa はき裂長さ，N は荷重サイクル，ΔG はエネルギー解放率範囲，C, m は実験的に 決定される定数である．与えられた半楕円形状の初期き裂形状（短半径 20 mm，長半径 100mm）に対して，静解析を実施し，き裂前縁においてエネルギー解放率を計算し，式(7)を用いてき裂進展量を計算する．ここでは，式(7)のき裂進展則のパラメータとして m=1.535，C=1.5891×10-5_を与え ΔN=5,000 として静応力解析を 10 ステップだけ繰り返すことによりき裂進展解析を実施した．初期状態と最終段階（10 ステップ後）における対称面上のき裂形状を図 5 に示す．図 5 に示すように，X-FEM を用いることにより，平板内を進展するき裂の前縁形状を有限要素メッシュと独立にモデル化でき，したがってメッシュ再分割をせずにき裂進展解析が実施できることがわかる．また，き裂最深部と表面部のき裂前縁位置と繰り返し回数との関係を図 6 に示す．最深部（90°位置）では 12mm 程度，表面部（0°位置）では 1mm 程度,き裂前縁が進展した結果が得られている．図 4 半楕円表面き裂を有する厚肉平板モデル Initial

Final (After 10step) Initial

Final (After 10step)

(6)

90 degree

0 degree

100 100.5 101 101.5 102 0 1 104 2 104 3 104 4 104 5 104 Po st io n [ m m ]

Number of load cycles 20 25 30 35 0 1 104 2 104 3 104 4 104 5 104 Po st io n [ m m ]

Number of load cycles

90 degree

0 degree

100 100.5 101 101.5 102 0 1 104 2 104 3 104 4 104 5 104 Po st io n [ m m ]

Number of load cycles 20 25 30 35 0 1 104 2 104 3 104 4 104 5 104 Po st io n [ m m ]

Number of load cycles

図 6 き裂前縁位置と繰り返し回数との関係 3. VCAD データ構造をベースとした構造解析プログラム（V-X3D）[8] VCAD で扱われるボリュームデータは，オクツリー構造を有するセル情報を基本とし，セルと形状定義面との切断面・切断点などの幾何情報に加えて物性情報を保持することができる．このような特長を有するボリュームデータは，人工物だけではなく自然物のモデル化にも有効であると考えられる．著者らは，損傷や欠陥を含む機械構造物や分布する特性（密度，ヤング率など）を有する生体骨などを直接モデル化して応力解析を効率的に実施することを目的として，VCAD が提供するオクツリー構造を有するボリュームデータを直接利用した応力解析手法の研究開発を行っている．現在，そのプロトタイプ版の構造解析プログラム V-X3D を開発中である．V-X3D の特徴として： z ボリュームデータを直接利用して応力解析を行うことができること． z X-FEM(拡張有限要素法)に基づいた応力解析手法であること． z V-CAD-Framework を介して V-CAD の機能を利用したアプリケーションであること．などをあげることができる． 3.1 解析手法離散化対象となる支配方程式は，通常のFEM と同様に連続体についての仮想仕事の原理式である．支配方程式を離散化するために，ボリュームデータのセル構造を利用する．図 7 に示すようにボリュームデータはオクツリー構造を有するが，セル分割の基本となる最低次(レベル 0)のセルを 8 節点直方体形状の有限要素に対応させる．形状定義面と交差するセル（境界セル）に対しては，拡張有限要素法(X-FEM)の解析技術を用いて剛性評価を行う．この場合，要素に対応するレベル0 のセルに連結する高次のセルを用いて体積積分を行う．

(7)

表面荷重は境界セルの切断面を利用して面積積分することによって離散化される．ボリュームデータを利用した体積積分，面積積分の概念を図 8 に示す．また，形状定義面において変位境界条件を与えるために，複数個のセル切断点での拘束条件式を導出しグループ化したうえで，節点自由度に関する多点拘束条件式を導出する．図 9 にその概念を示す．結果として得られるシステム方程式は通常のFEM と同様な形式となるので，既存の FEM で用いられる求解処理方法が利用可能である．

Inner cell & Level 0

Leaf of the octree under the boundary cell with Level 0

Level 0 Node generation

conventional FEM Inner cell & Level 0

Leaf of the octree under the boundary cell with Level 0

Level 0 Node generation

conventional FEM

図 7 VX3D におけるモデル化の概念

(a) Volume integral

(b) Area integral

(a) Volume integral

(b) Area integral

図 8 境界セルにおける体積積分と面積積分 Group 1 Group 3 Group 2 Group 1 Group 3 Group 2 図 9 境界セルにおける拘束条件評価点とそのグループ化

(8)

3.2 V-X3D の実装

V-X3D の開発のプログラム実装においては VCAD Framework を利用する．VCAD Framework が提供する API (Application Program Interface) 群を呼び出すことでボリュームデータ構造の詳細を意識することなく，ボリュームデータの生成，変更，取得などの処理が可能となる．提供されている API 関数の機能には • 形状データを入力 • VCARファイルの入出力 • 空間のセル分割 • Kitta Cubeの生成 • 指定したセル位置の幾何情報の取得などがあり，目的の処理に対する API 関数を呼び出すことで，ボリュームデータを利用したアプリケーションを効率的に開発することができる． 3.3 解析例機械部品を想定した形状を対象として開発中の V-X3D を用いて応力解析を行った．対象は工作機械の HSK ツールホルダであり，VOBJ 形式の表面形状データを，VCAD Framework の API 関数を用いてボリュームデータ(VCAR データ)に変換した．その際空間を 41x41x41 分割した．V-X3D は VCAR 形式のデータに入力し，形状を表示するとともに境界条件の付与，解析条件の設定，さらには結果の表示などを行うことができる．図 10 に V-X3D の実行画面の一例を示す．境界条件としては，ツールシャンク部のテーパと端面を完全拘束し，工具把持部内面の下面に下方向に圧力荷重を付与した．解析結果から得られる変形分布およびミ－ゼス応力コンタを表示した例を図 11 に示す．詳細は省くが，通常の FEM 解析結果などと比較してほぼ妥当な結果が得られている．図 10 V-X3D による境界条件の設定

(9)

図 11 V-X3D による解析結果（変形およびミ－ゼス応力）の表示 4．おわりに本稿では X-FEM による疲労進展解析事例および V-X3D による応力解析事例を示した．今後は，V-X3D に X-FEM によるき裂進展解析機能を組み込む予定である，図 12 にボリュームデータでモデル化されたブラケット形状を V-X3D に読み込み，さらに半楕円形状の表面き裂を設定した例を示す．このようなモデルについてのき裂進展解析を実施することが今後の課題である．図 12 V-X3D による半楕円表面き裂のモデル化

(10)

参考文献

[1] Belytschko T, Black T, Elastic crack growth in finite elements with minimal remeshing, Int. j. numer. methods eng.1999;45: 602-620.

[2] Moës, N, Dolbow J, Belytschko T., A finite element method for crack growth without remeshing, Int. j. numer. methods eng. 1999;46: 131-150.

[3] Moës N, Gravouil A, Belytschko T, , Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model, Int. j. numer. methods eng. 2002;53:2549-2568.

[4] Gravouil A, Moës N, Belytschko T, Non-planar 3D crack growth by the extended finite element and level sets-Part II: Mechanical model, Int. j. numer. methods eng. 2002;53: 2569-2586.

[5] Sukumar N., Moës N., Moran B., Belytschko T.: Extended finite element method for three-dimensional crack modeling, International Journal for Numerical Methods in Engineering 2000; 48, 1549-1570.

[6] Kase K., Teshima Y., Usami S., Ohmori H, Teodosiu C. , Makinouchi A., Volume CAD, Volume Graphics 2003 Eurographics / IEEE TCVG Workshop Proceedings, I. Fujishiro, K. Mueller, A. Kaufman (eds.), 2003; 145-150. [7] Kase, K. Teshima, T., Usami, S. Kato, M., Yamazaki, S. Ito, M., Makinouchi, A.: Volume CAD-CW-complexes

based approach, Computer Aided Design, 2005;37:1509-1520.

[8] 長嶋，今井：VCAD-API を用いた X-FEM に基づいた構造解析プログラム（VX3D）の開発，理研シンポジウムものつくり情報統合化研究（第 5 回）， V-CAD システム研究-「ものつくり」ツールの構築と「研究」の基盤ツールへの展開-講演予稿集(2005-6),90-99.

[9] Sethian HA, Level Set Methods and Fast Marching Methods: evolving interface in computational geometry fluid mechanics computer vision and material science, Cambridge, UK: Cambridge University Press, 1999.

[10] 長嶋，寺原，三浦：X-FEM による三次元弾塑性応力解析システムの開発，第 19 回計算力学講演会講演論文集,06-9 (2006-11), 581-582. [11] 作増，長嶋，三浦：X-FEM による残留応力場中のき裂進展解析，計算工学講演会論文集，12(2007-5), 415-416. [12] 作増，長嶋： X-FEM による三次元き裂進展解析, 日本機械学会 2007 年度年次大会講演論文集 Vol.1, 07-01 (2007-9), 21-22. [13] 長嶋，三浦：対称条件を考慮した X-FEM による三次元き裂解析, 機論 72-719, A(2006), 974-981.