COMPUTATIONAL FRACTURE MECHANICS
WHY?
- to compute fracture mechanics parameters (SIF, G) in 2D
and 3D configurations;
- to compute J integral and CTOD in elastic-plastic
analyses;
- to simulate crack growth (under general mixed-mode
conditions);
- to solve special problems: dynamic fracture, ductile
fracture, cohesive fracture, fracture at interfaces, ….
NUMERICAL METHODS:
To determine the distribution of stresses and strains in a
cracked body subject to external loads or displacements
(when close-form analytical solutions are not available)
- The Finite Difference Method;
- The Finite Element Method;
- The Boundary Element Methods.
Computational FM 1
TECHNIQUES FOR COMPUTING FRACTURE
MECHANICS PARAMETERS
POINT MATCHING APPROACHES:
1) STRESS MATCHING
2) DISPLACEMENT MATCHING
(The Displacement Correlation Technique)
3) STRESS FUNCTION MATCHING
ENERGY APPROACHES:
4) THE GLOBAL ENERGY RELEASE METHOD
5) THE STIFFNESS DERIVATIVE TECHNIQUE
(for FEM)
6) CONTOUR INTEGRATION
7) THE ENERGY DOMAIN INTEGRAL METHOD
8) THE MODIFIED CRACK-CLOSURE INTEGRAL
TECHNIQUE
Computational FM 2
POINT MATCHING APPROACHES:
Idea: to correlate numerical solutions for stresses or
displacements at specific locations with analytic solutions
depending on the stress intensity factors.
Hp: homogeneity, isotropy
linearly elastic material
1) STRESS MATCHING
hp: mode I loading.
KI =
lim
(
σ 22 2π r
r →0
)
(θ = 0)
2) DISPLACEMENT MATCHING
hp: mode I loading.
KI =
κ = 3−4ν
κ = (3−ν)/(1+ν)
2G lim 
2π
 u 2
κ +1 r → 0 
r



(θ = π )
(plane strain)
(plane stress)
better than 1): higher precision in the numerical calculation
of displacements (continuity between elements)
Computational FM 3
3) MATCHING STRESS FUNCTIONS
(The Boundary Collocation Method)
- The Airy stress function is expressed in terms of complex
polynomials. The coefficients of the polynomials are inferred
from nodal quantities. SIF's are inferred from the stress
functions. (highly cumbersome)
Drawbacks:
- High degree of mesh refinement is required for
engineering accuracy even for simple geometry, loading
and a single crack.
- Use of special elements at the crack tip that exhibit the
1/√r singularity is necessary (enriched elements, quarter
point elements).
Computational FM 4
2) THE DISPLACEMENT CORRELATION TECHNIQUE
(Ref. Chan, Tuba and Wilson, Eng. Fract. Mech., 2(1),1970)
Hp: homogeneity, isotropy
linearly elastic material
General form of displacement function near crack tip:
KI
θ 
θ
r
sin (κ + 1) - 2cos 2  + ....
2G 2π
2
2
K
r
θ 
θ
cos (κ - 1) - 2 sin2  + ....
- II
2
2G 2π
2 
u2 =
KI
r
θ 
θ
cos (κ - 1) + 2sin2  + ....
2
2
2G 2π
K
r
θ 
θ
sin (κ + 1) + 2 cos 2  + ....
+ II
2G 2π
2
2
u1 =
κ = 3−4ν
κ = (3−ν)/(1+ν)
(plane strain)
(plane stress)
Consider a FE mesh at crack tip:
Computational FM 5
Pure Mode I and plane stress:
u2 =
u1 =
KI r
θ  2
θ
sin 
- cos 2  + ....
G 2π
2 1+ ν
2
KI r
θ  1 -ν
θ
cos 
+ sin2  + ....
G 2π
2 1 + ν
2
taking θ = ± π and r = rAB:
K *I =
2π G(1+ ν )
(u2B - u 2C )
rAB
4
Pure Mode II and plane stress:
u2 =
K II
G
r
θ ν -1
θ
cos 
+ sin2  + ....
2π
2 1 + ν
2
u1 =
K II
G
θ  2
θ
r
sin 
+ cos 2  + ...
2
2π
2 1+ ν
taking θ = ± π and r = rAB:
K *II =
2π G(1+ ν )
(u1B - u1C )
rAB
4
Mixed Mode and plane stress:
Mode I and mode II problems uncouple at θ = ± π.
Consequently previous equations are still applicable.
Computational FM 6
For Quarter Points Singular Elements:
(Shih, de Lorenzi and German, 1976, IJF, 12, 647-651)
- If FE mesh is composed of quarter point elements at crack
tip:
K *I =
2π G
[4(u2B - u2C ) - (u2D - u 2E )]
rAD κ + 1
K *II =
2π G
[4(u1B - u1C ) - (u1D - u1E )]
rAD κ + 1
(from shape functions of element and analytical expressions
for COD and CSD).
Computational FM 7
FINITE ELEMENT MESH DESIGN
- In homogeneous, isotropic, elastic material, the stress field
at crack tip exhibits a 1/√r singularity.
- The FE mesh must be such that the stress field
characteristics are reproduced in the numerical solutions.
Solutions:
- Mesh refinement at crack tip
(slow convergence of local parameters)
- Use of “enriched” elements at crack tip (Wilson element;
Tracey element, …):
stress / displacement variations around crack tip are
embedded in the shape functions of element
(Tracey, 1971, Eng. Fract. Mech., 3, 255-266)
(edge displacement incompatibility with surrounding
elements)
- Use of Quarter Point Singular Elements:
in 8-nodes isoparametric quadrilateral elements, placing the
mid-side nodes on or outside the ¼ of the side causes the
Jacobian of transformation to become non positive definite.
The 1/√r singularity is obtained.
(Barsoum,1976, IJNME, 10, 25; Henshell and Shaw, 1975, IJNME, 9, 495507)
(only nodal coordinate input data are altered. Element
satisfies essential converge criteria)
Computational FM 8
- If the 8-nodes element is first degenerated to an
isoparametric triangle with 3 nodes at crack tip and then
distorted by moving the side nodes, improved results are
obtained: 1/√r singularity also within element not only along
edges.
- In 3D, analogous results obtained by distorting three
dimensional 20-nodes isoparametric elements
(Ingraffea and Manu,1980, IJNME, 15, 1427-1445)
Computational FM 9
QUARTER POINT ELEMENTS: PROPERTIES
Computational FM 10
Computational FM 11
ENERGY APPROACHES:
1) THE GLOBAL ENERGY RELEASE METHOD:
Hp: 2D linear elastic body.
From energy balance:
G =-
dW
da
W = total potential energy per unit width.
- Perform two analyses, one for crack length a, the second
for a +∆a, for constant load. Compute G from:
G =-
∆W ∆L ∆U ∆U
=
=
∆a ∆a ∆a ∆a
L = potential of the applied loads per unit width in a 2D
body;
U = elastic strain energy stored in the body per unit width
- For pure mode I, compute SIF from:
K I = EG
(plane stress)
EG
KI =
(plane strain)
1-ν 2
(hp: isotropy, homogeneity)
Computational FM 13
Drawbacks:
- two analyses required for KI computation. At least three
analyses for KI and KII computations. Four analyses for KI,
KII and KIII.
- SIF's not directly computed. Postprocessing of results is
required.
- For mixed-mode fracture it is very difficult to separate G
into its mode I and mode II components.
Advantages:
- no special crack-tip elements are necessary. Relatively
coarse meshes can be used.
- When used in conjunction with the Stiffness Derivative
Technique, only partial reanalysis is necessary.
Computational FM 14
2) THE STIFFNESS DERIVATIVE TECHNIQUE
(Parks, IJF, 10, 1974, 487-502)
- Formulated in terms of FE stiffness matrix. Not compatible
with boundary element analysis.
- Total potential energy of FE solution:
W = 1/2 uT[K]u - uTP
[K] = structure stiffness matrix;
u = vector of nodal displacements;
P = vector of applied nodal forces.
- Differentiate to obtain G:
∂uT
([K]u − P) - 1 uT ∂[K] u + uT ∂P
G =2
∂a
∂a
∂a
where first term on right hand side must be zero.
- If applied load is independent of crack length and pure
mode I:
K I2
1
∂[K]
u
G=
= - uT
E'
2
∂a
Computational FM 15
-If only elements surrounding the crack tip are modified by
crack advance ∆a (see figure), then [K] remains unchanged
in the region outside:
K I2
1
G=
= - uT
E'
2
∂[K ic ]
1 T
u
≅
u
∑ ∂a
2
i=1
nc
nc
(
)
1
∑ ∆a [Kic ]a+ ∆a − [K ic ]a u
i=1
[Kci] = stiffness matrix of the ith element surrounding the
crack tip (within Γ1 in the figure);
nc= number of such elements;
[Kci]a element stiffness matrix for original crack configuration;
[Kci]a+∆a element stiffness matrix for new crack configuration;
Computational FM 16
6) CONTOUR INTEGRATION
The J integral can be evaluated numerically along a contour
(in 2D) or a surface (in 3D) surrounding the crack tip. Then
in an elastic body:
G =J
Advantages:
- Path independence enables the user to calculate J at a
remote contour (surface) where numerical accuracy is
greater.
- No mesh refinement is required.
- J can be separated into its mode I and mode II
components to deduce KI and KII.
Drawbacks:
- Contour/surface integrals are difficult to implement
numerically and the accuracy is not high.
Solution:
- The Energy Domain Integral: the J integral is formulated in
terms of area (for 2D) or volume (for 3D) integrals.
Better accuracy. Easier to implement.
Computational FM 17
7) THE ENERGY DOMAIN INTEGRAL METHOD
(Shih, Moran and Nakamura, 1986, IJF, 30, 79-102)
The J integral for a 2D body is:

∂u
∂u 
J = ∫  Uddx 2 - T
ds  = ∫ Udn1 - σ ij i n jds
∂x1
∂x1  Γ
Γ
- Consider the close path Γn = Γ0 + ΓS+ + ΓS- - Γ .
- Introduce a weight function q(x1,x2) that is equal to unity on
Γ and zero on Γ0 and Γs.
The J integral along the new path is:

∂u 
J = - ∫  Udδ 1j - σ ij i  q n jds
∂x1 
Γn 
Computational FM 18
- Apply the divergence theorem to transform the integral
along the close contour into a domain integral:

 ∂Ud ∂  ∂ui  
∂ui ∂q 
∂q

σ
  q dA
J = - ∫ Ud
-σ
dA - ∫ 
 ∂x1 ij ∂x1 ∂x j 
 ∂x1 ∂x j  ij ∂x1  


A
A
where the second term vanishes for elastic problems.
- Implement the domain integral into a finite element code
(e.g. Franc).
Mode separation:
(Bui, J.M.P.S., 31(6), 439-448, 1983)
- The displacement field in the crack-tip region is
decomposed into symmetric and antisymmetric
components: u = usymm. + uantisymm.
- The J integral is evaluated twice: once with the symmetric
displacements to find J1 and a second time with the
antisymmetric displacements to find J2.
- For elastic problems:
GI = J1
GII = J2
SIF’s are found from GI and GII.
Computational FM 19
8) THE MODIFIED CRACK-CLOSURE
INTEGRAL TECHNIQUE
(Rybicki and Kanninen, 1977)
- The strain energy release rate G is estimated in terms of
the work done by the stresses ahead of the crack tip over
the displacements produced by the introduction of a virtual
crack extension (Irwin concept of crack-closure integral):
lim 2 ∆ 1
G =
σ 22 (x1,0) u2 ( ∆a - x1,0)dx 1
∫
I
∆a → 0 ∆a 0 2
a
lim 2 ∆ 1
G =
σ 12 (x1,0) u1( ∆a - x1 ,0)dx 1
II
∆a → 0 ∆a 0∫ 2
a
- Two numerical analyses are required to obtain the stress
field ahead of the crack before propagation and to compute
the displacement field after a virtual crack extension is
introduced.
Computational FM 20
- Simplification proposed by Rybicki and Kanninen:
the displacement field ahead of the crack is approximated
by the the displacement field behind the crack tip.
The problem is then solved with one single analysis step.
- In the finite element approach the crack closure integrals
are rewritten in terms of equivalent nodal forces and relative
nodal displacements (different expressions for each crack
tip element).
Advantages:
- no assumption of isotropy or homogeneity around the
crack tip is necessary.
Computational FM 21
MIXED-MODE CRACK GROWTH
- The direction of crack growth is defined using the
Maximum Circumferential Stress Criterion.
(Erdogan and Sih, 1963, J. Basic Engrg., 85, 519-527)
- The crack is propagated of ∆a along the calculated
direction.
- Remeshing: a new mesh is generated around the new
crack tip. Different techniques.
Franc3D code:
Solid Modeling Techniques to represent the desired
geometrical aspects of the structure explicitly. The
discretization is mapped onto this original solid model.
Remeshing: the solid model remains valid and the new
mesh is remapped into the solid model and inherits the
necessary attributes.
(see Franc3D manual for more details and references)
- New fracture parameters at a + ∆a are calculated.
Computational FM 22