Transcript IntroCompFracMech

J.Cugnoni, LMAF-EPFL, 2014
Stress analysis
Fracture mech.




Stress based criteria
(like Von Mises) usually
define the onset of
“damage” initiation in
the material
Once critical stress is
reached, what happens?
In this case, a defect is
now present (ie crack)
The key question is
now: will it propagate?
If yes, will it stop by
itself or grow in an
unstable manner.
Stress concentrator:
Critical stress is reached…
A crack is formed…
Will it extend further?
If yes, will it propagate abruptly
until catastrophic failure?

Stress singularity
Assumed
Crack extension dA
Crack propagation : stress intensity
factor
◦ Stress intensity factors KI, KII, KIII measure
the intensity of stress singularity at crack
tip.

“Stress intensity factor”:
◦ Constants of the 1/sqrt(r) term in stress
field at crack tip:

“Critical Stress intensity factor”:
◦ Maximum K that a material can sustain,
considered as a material property and
indentified in standard fracture tests.
Units: Pa/sqrt(m), symbol: KIc KIIc KIIIc
Fracture test: measure force at
failure and calculate KIc
From analytical solutions or FE

Crack propagation occurs if
K > Kc

with K=Kic +KIIc+KIIIc
In many materialy, propagation is
mode I dominated:
KI>KIc
Potential energy: P0=U0-V0

Crack propagation : an “energetic” process
◦ Extend crack length: energy is used to create a
new surface (break chemical bonds).
◦ Driving “force”: internal strain energy stored in
the system
Assumed
Crack extension dA

“Energy release rate”:
◦ Change in potential energy P (strain energy and
work of forces) for an infinitesimal crack
extension dA. Units: J/m2, Symbol: G
◦ measure the crack “driving force”
Potential energy: P1= P0 –Er
And Er = G*dA

“Critical Energy release rate”:
◦ Energy required to create an additionnal crack
surface. Is a material characteristic (but depends
on the type of loading). Units: J/m2, symbol: Gc
New crack surface dA:
Dissipates Ed=Gc*dA

Crack propagation occurs if G > Gc

Using numerical simulation method we can:
◦ Apply Finite Element, Finite Volume/Difference,
Boundary Element or Meshfree methods
◦ To obtain approximations of displacement, force,
stress and strain fields in arbitrary configuration.

And then? To apply fracture mechanics, we
need:
◦ to compute fracture mechanics parameters (SIF
K, G) in 2D and 3D configurations;
◦ to compute J integral in elastic-plastic analyses ;
◦ to simulate crack growth (under general mixedmode conditions);
FE simulation of a
compact tension
fracture test

Calculation of K (linear elasticity):
◦ Stress or displacement field matching (single / mixed mode)
◦ Indirectly from G (Interaction integrals in mixed mode)

Calculation of G:
◦ Finite difference of potential energy (linear, evtl. non linear)
◦ Compliance method (linear)
◦ Virtual Crack Closure Technique VCCT (linear, mixed mode)
◦ J-integral (non-linear)

Simulate crack growth
◦ to VCCT criterion, node release, remeshing or XFEM
◦ Cohesive elements / interfaces (Damage mechanics)

Idea:
◦ compare FE stress field at crack tip with theory,
◦ fit KI from the numerical stress value
◦ In LEFM (for r  0):  

(  0)   yy  K I /
2 r  O ( r ) (eq 4.36a)
Step 1: From FE: Extract stress field Syy at =0 (along crack direction)

Fitting method 1: Plot (as a function of r, for r  0)
K I  lim  yy
r0
2 r
Fit region

Fit a line over the quasi constant region of the plot, identify KI either as
average value or as the intercept at r=0
 yy

Fitting method 2: Plot

In theory the plot should be linear as:
 yy  K I (1 /

as a function of
x'1/
2 r
2 r )  O ( r )  K I x ' b
Fit a line over the quasi linear region of the plot, identify KI either as the
slope of the line
Note : the same method can be applied to get Kii from
shear component at theta=0

Idea:
◦ compare FE crack displacement opening field with theory,
◦ fit KI from the numerical displacement value
K I (  1)
◦ In LEFM (for r  0): u  (    )  u y 
4

r / 2
(eq 4.40b)
Step 1: From FE: Extract stress field Syy at =0 (along crack direction)
x
For plane stress:  
3 
1 
For plane strain:   3  4
 = shear modulus= E 2 (1   )
Fitting method 1: Plot (as a function of r=-x, for r  0)
K I (  1)
u y 2 / r 
4

Fit region

Fit a line over the quasi constant region of the plot, identify
either as the average value or as the intercept at r=0
K I (  1)
4

2 / r
K I (  1)
x' 1/

Fitting method 1: Plot uy as a function

In theory, the plot should be linear as : u 
y
4
1

2 / r  

Fit a line over the quasi linear region of the plot, the slope is
K I (  1)
4
x'
K I (  1)
4

Evaluation of K from stress field requires very fine mesh. To
better capture the 1/sqrt(r) stress singularity, one can use
singular quadratic elements with shifted mid side nodes at ¼
of edge length (Barsoum,1976, IJNME, 10, 25; Henshell and Shaw, 1975, IJNME, 9,
495-507)



Using 8 node second order quadrangular
elements, if we move the midside nodes at
¼ of edge, we make the Jacobian
transformation of the element singular as
1/sqrt(r) along the element edges.
By collapsing all nodes of one edge
(degenerate a quadrangle to a triangle),
the 1/sqrt(r) singularity covers then the
whole element.
Can be extended into 3D for 20 node
hexahedrons

For linear elastic isotropic materials, the
following relations link G to K:
GI 
 1
8
For plane stress:  

KI
3 
1 
G II 
2
8
K II
2
For plane strain:   3  4
G III 
1
2
K II
2
 = shear modulus= E 2 (1   )
For Mode I in plane stress and plane strain,
this simplifies to:
plane stress : G I 

 1
KI
E
2
(1   ) K I
2
plane strain : G I 
E
So knowing either K or G, the other can be
determined directly. This can be extended to
orthotropic material.
2
◦ compare the strain energy of several FE models with different crack
length, calculate the ERR as the derivative of potential energy P = U - F
◦ For a linear elastic material: the potential energy F = PD  2U thus
P = U – F = - U with the strain energy U 
◦ Compute ERR G  
P
a

U (a  Da )  U (a )
Da
1
 2
V
ij
 ij dV
or as the slope of U(a)
Different FE models for a0, a1, a2…
a
Equations are for a unit depth, if not the case, compute use U/b instead of U

Idea (similar to experimental test data reduction !)
◦ Calculate the compliance C(a)=D(a)/P(a) for different crack lengths a0,
a1 (several FE models required) and fit it with an appropriate function
◦ For a linear elastic material: the ERR can be obtained from the C(a)
2
2
P  dC 
P  dC 
G 
or
G





2 b  da  P  cst
2 b  da  D  cst
◦ Compute ERR
where dC/da is obtained from the fitted curve
P, D
Different FE models for a0, a1, a2…
a
(eq 3.41 & 3.43b)



Using the “J-Integral” approach
(see course), it is possible to
calculate the ERR G as G = J in
linear elasticity.
If we know the displacement and
stress field around the crack tip,
we can compute J as a contour
integral:
J-integral is path independent for
all continuum materials. But G
must be perpendicular to the
interface if dissimilar materials
are used.
tσ n
W=strain energy density
u = displacement field
s = stress field
G = contour: ending and starting
at crack surface


J-integrals are usually computed from a volume/surface integral in
FE, for example in Abaqus with q = virtual crack extension vector:
With
on G and
on C. Using divergence and
equilibrium equations we can obtain:
J-integral can be extend in 3D by computing J on several “slices” of the model.
However, there are notable 3D effects affecting crack propagation:
in 3D, the external surface is free of normal stress, so in plane stress state
For thick specimens, the center of the specimen is closer to plane strain state
S22 on external surface
Plane stress
S22 on symmetry plane
Plane strain => higher stresses
As a consequence, J-integral and G are not constant across the width.
This will lead to a slightly curved crack front during propagation
-Need to be careful about specimen
size effects when characterizing G or K
- In 2D FE simulation: use plane stress
for very thin specimens, and plane
strain for thick ones.
G is max in the center
=> propagate earlier
G is min on the side
=> propagate last
Results obtained on a relatively coarse 3D mesh to highlight
which methods are less mesh sensitive

External surface, plane stress
Stress fit1
K1
3870
(units: Mpa*sqrt(mm))
K1
Stress fit2
4345
Displacement fit1 Displacement fit2
4380
5079
From Jintegral
4337
K from
Abaqus
4540
Inside, plane strain
(* please note that the FE data for u2 and s22 were extracted on free surface, a source of error)
K from
Stress fit1
Stress fit2
Displacement fit1 Displacement fit2 J-integral
Abaqus
4345
3870
4813
5582
5144
5140
Strain energy,
forward finite
difference
G
(units: mJ/mm2)
322
Strain energy,
backward finite
difference
Strain en.
centered finite
difference
Strain energy fit
282
302
302.5
Compliance
J integ. Avg
method
301.5
Conclusion:
K determination is more sensitive to mesh !
G calculation using compliance, J-integral or strain energy fit are
the most reliable
296




Can be calculated in elasticity
/ plasticity in 2D plane
stress, plane strain, shell and
3D continuum elements.
Requires a purely
quadrangular mesh in 2D
and hexahedral mesh in 3D.
J-integral is evaluated on
several “rings” of elements:
need to check convergence
with the # of ring)
Requires the definition of a
“crack”: location of crack tip
and crack extension direction
Crack tip and
extension direction
Quadrangle
mesh
Crack plane
Rings 1 & 2









Create a linear elastic part, define an “independent” instance in Assembly
module
Create a sharp crack: use partition tool to create a single edge cut, then in
“interaction” module, use “Special->Crack->Assign seam” to define the crack
plane (crack will be allowed to open)
In “Interaction”, use “Special->Crack->Create” to define crack tip and
extension direction (can define singular elements here, see later for more
info)
In “Step”: Define a “static” load step and a new history output for J-Integral.
Choose domain = Contour integral, choose number of contours (~5 or more)
and type of integral (J-integral).
Define loads and displacements as usual
Mesh the part using Quadrangle or Hexahedral elements, if possible
quadratic. If possible use a refined mesh at crack tip (see demo). If singular
elements are used, a radial mesh with sweep mesh generation is required.
Extract J-integral for each contour in Visualization, Create XY data -> History
output.
!! UNITS: J = G = Energy / area. If using mm, N, MPa units => mJ / mm2 !!!
By default a 2D plane stress / plane strain model as a thickness of 1.
See demo1.cae example file

To create a 1/sqrt(r) singular
mesh:
◦ In Interaction, edit crack definition and
set “midside node” position to 0.25
(=1/4 of edge) & “collapsed element
side, single node”
◦ In Mesh: partition the domain to create
a radial mesh pattern as show beside.
Use any kind of mesh for the outer
regions but use the “quad dominated,
sweep” method for the inner most
circle. Use quadratic elements to
benefit from the singularity.
◦ Refine the mesh around crack tip
significantly.
39
Regular mesh
Singular Mesh
38.8
38.6
38.4
38.2
38
37.8
37.6
Contour 1
Contour 2
Contour 3
Contour 4
Contour 5
•Perform a stress analysis
Stress analysis
•Locate stress critical regions
•Assume the presence of a defect in those regions (one at a time)
•Consider different crack lengths and orientation
Crack
analysis
•For each condition, check if the crack would propagate and if yes if it is
stable or not
•Define operation safety conditions: maximum stress / crack length,…
before failure occurs
Design
evaluation
•Define damage inspection intervals / maintainance plan

Abaqus tutorials:
◦ http://lmafsrv1.epfl.ch/CoursEF2012

Abaqus Help:
◦ http://lmafsrv1.epfl.ch:2080/v6.8
◦ See Analysis users manual, section 11.4 for fracture
mechanics

Presentation and demo files:
◦ http://lmafsrv1.epfl.ch/jcugnoni/Fracture

Computers with Abaqus 6.8:
◦ 40 PC in CM1.103 and ~15 in CM1.110