Crack modelling using
zero-thickness interface elements
Nguyen Vinh Phu
The University of Adelaide
Introduction
There exists problems in which the
crack path is known in advance...
delaminated
composites
inter-granular cracking in
polycrystals
Introduction (cont.)
...then the use of interface elements is recommended
• easy to implement (2D, 3D)
• available in major FE packages: ABAQUS, LS-DYNA
Alternatives:
XFEM/GFEM
Cohesive crack model
Governing equations
(strong form)
Constitutive equations
deformation
separation
Cohesive crack model
Weak form
new term
where
(skipped for static problems)
Discrete equations
upper face
lower face
N1 = 0.5(1 ⇠)
N2 = 0.5(1 + ⇠)
Solid elements
Discrete equation
Discrete equations (cont.)
assembled
Static problems
ext int coh
f =f +f
Linearization (Newton-Raphson)
transformation matrix
Common interface elements
Solid elements Q4/T3 Q8/T6
2D
H8
3D
Numerical integration
It has been observed numerically that integrating the
internal force and stiffness matrix of interface elements
using the standard Gauss rule led to oscillatory response
[de Borst, IJNME, 1993].
Newton-Cotes integration scheme
for interface elements
Cohesive laws
Mode I Bilinear cohesive law (traction-separation law)
1 1 tensile strength
GIc = tc f = k c f
2 2
2GIc tc
= =
f
k c
c
k fracture energy
elastic stiffness k
t = (1 d)k
8
>
< f
if loading
d= max
>
: f
otherwise
max ( f c) + c
crack initiation
c
= crack propagation
f c
Implementation aspects
• How to generate interface element meshes?
F
• Solution control
cannot pass the snap-through
- load control:
NO
- displacement control
- path following control
(arc-length methods)
u
cannot pass the snap-back
Mesh generation
A C++ code was written to
• read a Gmsh mesh file
• double nodes along a path defined by the user
• modify the solid elements involved and
• generate interface elements
Some application examples
Path following method
Riks 1972 ext
load factor
f = g reference load vector
Newton-Raphson (u, ) arc-length/constraint function
where
1 1
uI = K r, uII = K g
correction
Energy control
Z
✏ = Ba f int
= B T
Gutierrez 2004 ⌦ equilibrium
Z Z
1 T 1 T T 1 T int 1 T
V = ✏ = a B = a f = a g
2 ⌦ 2 ⌦ 2 2
Energy release rate V̇ T
ȧ g
G>0
predefined amount of energy
Arc-length function to be released [Nm]
forward Euler
straight line
Assumption: secant unloading!!!
E
Indirect displacement
control [de Borst 1986]
Indirect displacement control
(||uA uB || , l) local quantity!!!
SEN beam
imagine what if
there are 2 cracks???
Advantages of
energy control SEN beam
Energy control
- fast to evaluate
- global quantity
multiple cracks
FP van der Mer
EFM, 2008
Solution procedure
FP van der Mer
load energy control EFM, 2008
control
f int + f coh
Numerical examples
• Simple tests (to debug code)
• Material interface debonding
• Multi-delamination of a composite DCB
• Delamination of the DCB
Simple test (2D)
plane strain
Simple test (3D)
• As previous 2D example
• Thickness: 50
• Solved with Hex8 and Tet4
Debonding of a
material interface
interface: linear 4-node
continuum: T3 elements
Cohesive law: Xu-Needleman
Multi-delamination
Multi-delamination
Matrix kinking
interface elements everywhere
except for the hard inclusion
“Discontinuous Galerkin/extrinsic cohesive zone modeling: Implementation caveats and applications in
computational fracture mechanics”,VP Nguyen, Engineering Fracture Mechanics, 2014.
Discontinuous Galerkin and
cohesive interface elements
k: sufficiently large not to reduce the compliance of solid
Weibull distribution
(two-parameter version)
(three-param. version)
Stresses smaller than sigma_m, no failure: P=0
Probability Density Function (PDF)
Account for the fact that larger samples more prone
to failure
Weibull distribution: implementation
where rand is a random number between 0 and 1.
Weibull distribution: implementation
where rand is a random number between 0 and 1.
double xrand;
double value, f2t;
0.02
std::random_device rd; m=5
std::mt19937 gen(rd()); m=15
m=40
std::uniform_real_distribution<> dis(0, 1); 0.015
int ipCount = 2;
PDF [−]
int elemCount = count / ipCount; 0.01
for ( int i = 0; i < elemCount; i++ )
{ 0.005
xrand = dis(gen);
value = pow(-log(xrand),1/m_) ;
f2t = f2t_ * value + sigmaMin_; 0
0 500 1000 1500
f2ts_.pushBack ( f2t, ipCount ); σf [MPa]
}
each IE has the same strength
Things to explore
• New cohesive laws
• New (better) interface element formulations
(current element technology does not allow
industrial applications to be realized.)
• Mesh topology such that mesh bias can be avoided.
