1. No category

XFEM解析システムと
複合材料構造解析への応用
Development of an analysis system
based on the extended finite element method
and application to structural analyses of composite
長嶋 利夫（上智大）
Toshio NAGASHIMA (Sophia University)
2014年7月8日（火）第3回航空機CAE研究会
1
背景１
拡張有限要素法(XFEM）
メッシュと独立にき裂を定義可能．
要素再分割なしにき裂進展解析が可能．
積層複合材料構造
配管
岩石
by courtesy of CRIEPI
2
背景２
アプローチ
長所
線形破壊力学
•
•
比較的粗いメッシュが利用可能
線形解析
•
•
初期き裂の設定が必要
エネルギー解放率、K値評価が必要
•
•
応力に基づくき裂の発生
エネルギーに基づくき裂の進展
•
•
詳細なメッシュが必要
（材料）非線形解析
Linear Elastic Fracture
Mechanics: LEFM
結合力モデル
Cohesive Zone model:
CZM
短所
結合力き裂を用いた応力解析
~
t2
ti = ti
bi
Γu
x2
o
x1
ui = ui
t Imax
Γσ
A
Γ
~
x2
~
t1
Γ
+
GIC
Γ
~
x1
O
m2δ
O
m2δ
t IImax
u~
F 2
Γ−
O
GIIC
(a) Mode I
GIIC
m1δ O
m1δ
F
u~1
−t IImax
(b) Mode II
3
背景３
NR法で収束解が得られる限界
TCT(Transverse Crack Test)
500
tc t
層間はく離：インターフェース要素
樹脂割れ：XFEMによる結合力き裂
averaged stress [MPa]
400
300
200
100
0
0
0.2
0.4
0.6
0.8
1
Displacement [mm]
非線形有限要素法解析で標準的に用いられるニュートンラプソン法による
陰的解法では収束解が得られない場合がある．
4
内容
XFEMによるき裂、損傷進展解析システムを開発し、
複合材料積層板へ適用する．
•
•
•
•
•
背景、内容
解析手法の概要
開発システム
数値解析例
まとめ・今後の課題
5
本研究で用いる方法(方針）
• 層間はく離など予め発生位置が分かっているき裂は通常の
インターフェース要素でモデル化
• 発生位置がわかないき裂はXFEMでモデル化
• き裂面に結合力モデル（CZM)を導入可能
(解析対象によって様々なCZMを選択的に利用）
• 原則としてヘビサイド関数だけを拡充し、
必要があればTIP要素を導入
• 陰解法あるいは陽解法を利用
• 二次元問題：三角形要素
準三次元問題：五面体要素
三次元問題：四面体要素
6
レベルセット関数による二次元き裂の表現
x
φ >0
n
+
φ
ψ
crack tip
n′ +
x
φ < 0 crack line
ψ <0
ψ >0
Γ
φ (x) = min x − x sign(n + ⋅ (x − x))
x∈Γ
ψ (x) = min x − x sign(n′+ ⋅ (x − x))
x∈Γ′
7
二次元X-FEMにおけるき裂のモデル化
ψ (x)
φ (x)
Enriched node with the Heaviside function
Crack line
Crack tip
8
レベルセット関数φによる要素の分類
normal
φ
+
cut
φ
normal
φ
φ φ−
+
φ
+
−
zero
0
φ+
φ
zero
−
0
φ
−
zero
00
−
φ
cut
φ φ
−
φ−
0
0
0
cut
φ+
φ φ+
−
zero
0
+
−
φ
+
zero
0 φ+
φ−
φ+
0
9
レベルセット関数φ,ψによるCUT要素の分類
φ
cut
φ
+
cut
+
ψ
ψ
−
φ
+
φ
cut
φ
ψ
−
φ
−
φ
−
φ
−
+
ψ−
+
φ
−
0
+
φ
−
−
+
ψ
−
φ
+
φ
+
−
φ
ψ+
ψ
φ
φ
ψ+
φ−
+
ψ+
+
−
ligament
ψ−
−
φ
+
ligament
−
φ
tip
ψ
φ
ψ
+
φ
ligament
+
ψ
tip
ψ−
−
φ
tip
φ+
φ
−
ψ+
0
ψ+
φ
−
0
10
三次元XFEMにおける(内部）き裂のモデル化
x
P
φ
P
nS
x
ψ
nC
=
φ ( x)
min
x∈S ∪ S EXT
Γ
S
φ
S EXT
ψ
( x)
x − x sign(n S ⋅ (x − x )) ψ=
min
x∈Γ
x − x sign(nC ⋅ (x − x ))
11
レベルセット関数φによる
テトラ要素の分類（一部）
N
P
P
P
P
P
N
Z
N
(2 1 1)
P
N
N
(3 0 1)
P
(1 0 3)
N
N
N
Z
(2 1 1)
N
P
P
P
P
N
P
N
(2 0 2)
(2 0 2)
P
Z
(2 1 1)
P
12
半構造三次元五面体要素
Normal
Cut
Tip
extrusion
r3
x2
x1
Enriched node with the Heaviside function
u h (=
x, r3 )
3
∑ L (x)( N
I
I 1
B
( r3 )u I + NT ( r3 )u I +3 ) + ∑ FI (x)( N B ( r3 )a I + NT ( r3 )a I +3 )
I ∈J
FI (x) ≡ LI (x)(H (φ (x)) − H (φ (x I )))
N B ( r3 ) ≡ (1 − r3 ) / 2, NT ( r3 ) ≡ (1 + r3 ) / 2
13
内挿関数
3
2D
3
u ( x) =
∑ LI u I + ∑ LI ( H (φ (x)) − H (φ (x I )))a I
h
=
I 1=
I 1
3
φ (x) = ∑ L I φ I
I =1
3
u (x, r3 ) ∑ LI (x)( N B (r3 )u I + NT (r3 )u I +3 )
=
Quasi-3D
I =1
h
3
+ ∑ LI (x)( H (φ (x)) − H (φ (x I )))( N B (r3 )a I + NT (r3 )a I +3 )
I =1
3
φ (x) = ∑ L I φ I N B ( r3 )= (1 − r3 ) / 2 NT ( r3 )= (1 + r3 ) / 2
I =1
4
3D
4
u ( x) =
∑ LI u I + ∑ LI ( H (φ (x)) − H (φ (x I )))a I
h
=
I 1=
I 1
4
φ (x) = ∑ LI φI
I =1
14
解析手法
離散化式
 (t ) + CU (t ) + Q(U(t )) =
MU
F(t )
陰的静解析
Q(U(t )) = F(t )
K T ∆U = F − Q(U)
時間積分法：Newmark-β法（β=0.25, γ=0.5)
 (t + ∆=
 (t ) + ∆t γ U
 (t + ∆t ) + (1 − γ )U
 (t )
U
t) U
陰的動解析
{
}
 (t ) + ( ∆t )2  1 − β  U
 (t ) + β U
 (t + ∆t ) 
=
t ) U(t ) + ∆tU
U(t + ∆



 2

2
 − CU

(K T + (1/ β ( ∆t ) )M + (γ / β∆t )C)∆U = F − Q − MU
陽的動解析
時間積分法：中央差分法
 1
1 
1
1


M
C
U
F
Q
M
U
U
CU(t − ∆t )
+
(
t
+
∆
=
t
)
(
t
)
−
(
t
)
+
(2
(
t
)
−
(
t
−
∆
t
))
+
2
 ( ∆t )2
2∆t 
2
∆
t
( ∆t )


対角化
条件つき安定： ∆tcrit <
2
ωmax
(
1+ ξ 2 − ξ
)
15
XFEM解析における集中質量マトリクス
Uu 
Na ]  
U a 
3
u ( x) =
∑ LI u I + ∑ LI (H (φ (x)) − H (φ (x I )))a I ≡ [ Nu
h
I ∈J
=I 1
M uu ≡ ∫ ρ NTu N u dA, M ua ≡ ∫ ρ NTu N a dA,
A
Consistent mass matrix
M cons
M uu
= T
 M ua
A
M ua 
M aa 
Lumped mass matrix
M lump
0
 diag (M uu )

=

M
0
diag
(
)
aa 

Menouillard, T., Rethore, J., Combescure A., Bung H.: Efficient explicit time stepping
for the eXtended Finite Element Method (X-FEM), International Journal for
Numerical Methods in Engineering, 68, (2006), pp. 911-939.
16
2D-XFEM code
Code
NLXT2D
Latest version
3.140609
Development language
ANSI-C
Discretization method
eXtended Finite Element Method
Analysis type
Static, Implicit Dynamic, Explicit Dynamic
Element type
Three-node triangle element (Constant Strain Triangle)
Four-node interface element considering CZM
Material type
Isotropic, Orthotropic
Enrichment type
Heaviside
Method use to solve the
system equation
Direct method:
Skyline method
PARDISO (Intel Math Kernel Library)
Evaluation J-integral and SIF
Domain Integral Method, M-integral
17
Quasi 3D-XFEM code
Code
NLXP3D
Latest version
3.140707
Development language
ANSI-C
Discretization method
eXtended Finite Element Method
Analysis type
Static, Implicit Dynamic, Explicit Dynamic
Element type
Six-node pentahedral element
Six-node interface element considering CZM
Material type
Isotropic, Orthotropic
Enrichment type
Heaviside
Method use to solve the
system equation
Direct method:
Skyline method
PARDISO (Intel Math Kernel Library)
Evaluation J-integral and SIF
N.A.
18
3D-XFEM code
Code
NLXT3D
Latest version
3.140702
Development language
ANSI-C
Discretization method
eXtended Finite Element Method
Analysis type
Static, Implicit Dynamic, Explicit Dynamic
Element type
Four-node tetrahedral element
Six-node interface element considering CZM
Material type
Isotropic, Orthotropic
Enrichment type
Heaviside, Asymptotic basis
Method use to solve the
system equation
Direct method:
Skyline method
PARDISO (Intel Math Kernel Library)
Evaluation J-integral
Domain Integral Method
19
３D-XFEMによる解析例
20
開発システム
Crack definition file
Finite element model file
*.dat
*.crk
NLXT2D/NLXP3D/NLXT3D
*.out
Analysis results file
*.log
Log file
*.vtf
Analysis results file
for visualization
21
Material properties and
cohesive zone parameters
Laminate
Cohesive zone
EL [GPa]
161
GIC [N/mm]
0.2
ET [GPa]
11.38
GIIC[N/mm]
1.0
GLT [GPa]
5.17
σImax [MPa]
60
GTT [GPa]
3.98
σIImax [MPa]
90
νLT
0.32
kI [N/mm3]
1.0 x106
νTT
0.436
kII [N/mm3]
1.0 x106
ρ [Ns2/mm4]
1.6 x 10-9
α
1.0
2
Conditions for crack initiation
Conditions for crack propagation
2
 < σ >   τ 
1

 +
 =
σ
σ
 Im ax   II m ax 
 GI

 G IC
α

G
 +  II

 G IIC
α

 = 1

22
TCT(Transverse Crack Test)
GII =
σ2
tc
t
4 EL t − tc
σ cr = 2 GII EL
Tension load
t − tc
ttc
Tension load
tc
t
Developments in the science and technology of
composite materials,ECCM4 (1990)
σ cr = 2 1.0 ×1.61×105
4−2
= 401.24 MPa
4 ×2
23
Numerical results for TCT specimens
S11 [MPa]
d
Number of nodes 23958
Number of pentahedral elements 30720
Number of interface elements 7680
24
Numerical results for TCT specimens
Displacement vs. average stress
700
2D/Static
2D/Explicit
3D/Static
3D/Explicit
averaged stress [MPa]
600
500
Ref.
σ cr = 401.24 MPa
400
300
200
100
0
0
0.2
0.4
0.6
0.8
1
Displacement [mm]
25
Numerical results for TCT specimens
(Tetra model)
S11 [MPa]
d
Number of nodes 23958
Number of tetrahedral elements 92160
Number of interface elements 7680
26
Numerical results for TCT specimens
Displacement vs. average stress
(2D/Q3D/3D)
700
2D
Q3D(PENTA)
3D(TETRA)
averaged stress [MPa]
600
500
400
300
200
100
0
0
0.2
0.4
0.6
0.8
1
Displacement [mm]
27
Energy Balance (MS 104)
Energy [Nmm]
3 10
4
2.5 10
4
2 10
4
1.5 10
4
Uint
Ukin
1 104
5 10
3
0
0
5 10-2
1 10-1
1.5 10 -1
2 10-1
time [sec]
28
Fracture analyses for composite specimens
Hallett, S. R. et al, Modelling the interaction between matrix cracks and delamination damage in scaled quasiisotropic specimens, Composites Science and Technology 68 (2008) 80–89.
Jiang,W-G, et al, A concise interface constitutive law for analysis of delamination and splitting
in composite materials and its application to scaled notched tensile specimens,
Int. J. Numer. Meth. Engng 69 (2007) 1982–1995.
29
CFRP laminate specimen under tension load
Length：120 mm
Width: 32 mm
Thickness: 4 mm
[45/90/-45/0]s
Number of nodes 31944
Number of pentahedral elements 30720
Number of interface elements 23040
30
Modeling of matrix cracks by enrichment nodes
45 deg Layer
90 deg Layer
-45 deg Layer
31
Stress distributions and deformation
S11 [MPa]
800
Implicit/PENTA
Explicit/PENTA
d
Average stress [MPa]
700
600
500
400
300
200
100
0
0
0.5
1
1.5
2
LPD [mm]
Number of nodes 31944
Number of pentahedral elements 30720
Number of interface elements 23040
32
Comparison of numerical results
800
Implicit/PENTA
Explicit/PENTA
Average stress [MPa]
700
600
500
400
Hallett, S. R. et al, Modelling the interaction between matrix cracks
and delamination damage in scaled quasi-isotropic specimens,
Composites Science and Technology 68 (2008) 80–89.
300
200
100
0
0
0.5
1
1.5
2
LPD [mm]
Iarve, E. V. et al, Mesh-independent matrix cracking and
delamination modeling in laminated composite, IJNME (2011)
33
FE models with various scales
m
Stacking sequence
dimension
Thickness
[mm]
1
[451/901/-451/01]s
30 x 8
1
2
[452/902/-452/02]s
60 x 16
2
4
[454/904/-454/04]s
120 x 32
4
8
[458/908/-458/08]s
240 x 64
8
m =1
m=2
m=4
m=8
Number of nodes 31944
Number of pentahedral elements 30720
Number of interface elements 23040
34
Comparison of numerical results
for various scaling factors
1000
894
Averaged stress [MPa]
800
631
600
480
400
383
m=1
m=2
m=4
m=8
200
0
0
0.5
1
1.5
2
LPD [mm]
MS 104 , Δt=1.0x10-7, t=0.1sec, 40mm/sec
Hallett, S. R. et al, Modelling the interaction between matrix cracks
and delamination damage in scaled quasi-isotropic specimens,
Composites Science and Technology 68 (2008) 80–89.
35
Finite element models of composite specimen
with an open hole under tension load
Number of nodes 6792
Number of pentahedral elements 6280
Number of interface elements 4710
36
Stress distribution and Damage propagation
S11
[MPa]
d
37
FE models with various thicknesses
t=4 mm
t=2 mm
Input data for NLXP3D
t=1 mm
4 0.0
0.5 1
0.5 1
0.5 1
0.5 1
0.0
-45.0
90.0
45.0
1
1
1
1
2
2
2
2
38
Comparison of numerical results
for various dimensions
500
explicitH4
explicitH2
explicitH1
averaged stress [MPa]
400
300
200
100
0
0
0.2
0.4
0.6
0.8
1
Strain [%]
Jiang,W-G, et al, A concise interface constitutive law for analysis of delamination and splitting
in composite materials and its application to scaled notched tensile specimens,
Int. J. Numer. Meth. Engng 69 (2007) 1982–1995.
39
400
explicitH4
350
averaged stress [MPa]
300
250
200
150
100
50
0
0
0.2
0.4
0.6
0.8
1
Strain [%]
40
まとめ・今後の課題
XFEMに基づき開発した，二次元，準三次元，および
三次元き裂進展解析システムの概要と
それらを用いたCFRP積層構造についての数値解析
例を示した．
今後の課題




任意の位置で発生するき裂のモデル化
繊維破断のモデル化（連続体損傷モデルの導入）
リスタート機能
高速化（並列化、低減積分、陰解法と陽解法との組み合わせ、
合理的なマススケール設定）
 非線形性（幾何学的非線形）の考慮
41
 実問題への適用および検証
ご清聴ありがとうございました．
42