Numerical Fracture Mechanics

Second Paper
This is Chapter 5 of the invited lecture Problems in Fracture Mechanics by

S. Sedmak, J. Jari and M. Berkovi, presented at the 15th Yugoslav Congress of
Theoretical and Applied Mechanics, held in Kupari, 1981. The original paper is
published in Congress Proceedings (in Serbian), Vol. 1, pp. 67-118, and represents a continuation of the First Paper.
Ovde je Poglavlje 5 predavanja po pozivu Problemi mehanike loma autora
S. Sedmaka, J. Jaria i M. Berkovia, izloeno na 15. Jugoslovenskom kongresu
teorijske i primenjene mehanike, koji je odran u Kuparima, 1981. Originalna
verzija rada objavljena je u Zborniku radova sa Kongresa, na srpskom jeziku,
Vol. 1, str. 67-118, i predstavlja nastavak prvog rada.
Mladen Berkovi
NUMERICAL METHODS IN FRACTURE MECHANICS

Original scientific paper, UDC: 539.42:519.6
From the text presented in previous chapters, we can
conclude that the fundamental problem in fracture mechanics is the assessment of stress distribution near the sharp
crack tip. Stress distribution is determined as well as in
solid body mechanics in general, by analytical and numerical methods. Some experimental methods are described in
Chapter 4 with emphasis on standard specimen testing.
Analytical methods described in Chapters 2 and 3, apart
from the qualitative insight of fracture phenomenology,
enable obtaining an accurate solution for a certain number
of simple problems. Numerical methods represent the only
acceptable alternative for most fracture mechanics applications due to the fact that it is related to the solid body of
highly complex forms. Despite the abundance of numerical
methods, from which many are suitable for solving certain
fracture mechanics problems, in this paper we will focus on
the finite element method. The reason for this lies in the
generality and flexibility of this method which, combined
with application of up-to-date computers, has lead to an
industrial revolution in the field of numerical calculation
and design in general.
elastic /1,2/ and elastic-perfectly plastic /3/ problems,

special finite elements with correct stress and strain singularity are used. However, in the general non-linear case,
such approach is pointless and not very useful for special
problems, hence we will not consider it any further. Instead,
we will consider the effective possibilities of determining
certain fracture mechanics parameters by using standard
software, which is more or less available to every finite
element method user.
As we know, the application of this method involves,
above all, the forming of an adequate mesh in which the
observed object is divided into a number of topologically
uniform elements. In case of standard finite element method
programme application on fracture mechanics problems, a
significant refinement of the mesh around the crack edge is
necessary (see First Paper, Fig. 3).
Application of this method on a considered problem
results in displacement of mesh nodes and hence in the
determination of stresses. Once these quantities are known,
we can determine the fracture mechanics parameters by
using Eqs. A1 to A10, given in the Appendix (most of the
Eqs. given here are cited from original paper Proceedings,
please use this reference if necessary).
APPLICATION OF FINITE ELEMENT METHOD IN

DETERMINING STRESS AND DISPLACEMENT
AROUND THE CRACK TIP
Determining stress intensity factors based on given stress

values
First of all, we should emphasize that in fracture mechanics, general finite element method programmes are
being used, and that there are no differences between the
finite element method and the other solid-body mechanics
problems when defining and solving equations as well as in
selecting material models. When applying this method in
INTEGRITET I VEK KONSTRUKCIJA
Vol. 4, br. 2 (2004), str. 6366
Considering Eqs. A1, A3 and A7, we can conclude that

for given coordinates r and , and stresses x, y, z, xy, xz
and yz, the parameters KI, KII and KIII are solvable. For
practical reasons we will assume that = 0, which results
from these equations, in respect:
63
STRUCTURAL INTEGRITY AND LIFE

Vol. 4, No 2 (2004), pp. 6366
Numerical Methods in Fracture Mechanics
The problem is reduced to calculating the J integral,

which is given in Eq. A8 in its vector form. Rices J integral is the component of A8, which can also be shown with
the following expression /5/:
u
(5)
J = W dx 2 ij n j 1i ds
x
K I = y 2 r , K II = xy 2 r , K III = yz 2 r (1)
When effectively determining these values, we usually
calculate several of them along the direction = const, and
then evaluate the most convenient among them by extrapolating for r = 0 (First Paper, Fig. 6). Such a procedure, although simple, requires a certain amount of time and effort
in graphical elaboration of results.
The use of Eq. [5] requires the crack to lie in the x3x1
plane. For a linear elastic material with a symmetrical stress
tensor, the strain energy per unit volume is:
1
(6)
W = ij ui, j
2
where ui,j denotes the displacement gradients:
u
ui, j = ij
(7)
x
Taking into account that the contour surface is parallel to
the x3 axis, the contours normal vector coordinates are:
Determining stress intensity factors based on given displacement values

This procedure is analogous to the above and expressions used are A2, A4 and A6, respectively. From practical
reasons, we will assume that = , which means that we
will observe displacements on the crack surface, hence:
u1 2
u2 2
, K II =
, K III = u3
(2)
+1 r
+1 r
2r
Generally, this approach is a bit more accurate then the
previous due to higher displacement accuracy compared to
the stresses when using FEM.
Determining stress intensity factors based on energy release
rate under crack growth
When this approach is being used, stress intensity factors
may be explicitly determined only for pure strain forms (see
First Paper, Fig. 1); e.g. if KII = KIII = 0 it follows in accordance with /4/ that:
KI =
dx 2
dx1
, n2 =
(8)
ds
ds
Hence the index j in Eqs. [5], [6] and [7] takes the values
j = 1,2, while i = 1,2,3. The previous statement does not conflict with the assumption of stress and strain homogeneity
along the x3 axis. It also allows us to treat the shear strain
(see First Paper, Fig. 1) by using the J-integral.
In most general cases of a linear elastic material which
applies to plane strain, the constitutive equations are given
in the following form:
2 ij kl
ij =
v + (1 2v ) ik jl u( k, l )
(9)
1 2v
In case of plane strain, this expression is reduced to:
2 ij il
ij =
v + (1 v ) ik jl u( k, l )
(10)
1 v
1
where u( k, l ) = uk, l + ul, k . In both expressions j = 1,2
2
and k, l, i = 1,2,3 with the condition that uk,3 = 0. By placing
[6], [8] and [9] or [10] into [5] and using some algebra, we
obtain the expression for J integral under plane and antiplane strain or stress as:
n1 =
8
(3)
G
+1
according to A9 and A10. Concerning the strain energy
release rate G itself, we will determine it by using two
consecutive finite element method analyses, for two close
crack lengths which differ by a. Based on /5/ we can write
approximately:
dU
U
G=
=
hda
ha
where U/h represents strain energy per unit thickness h, for
plane stress state, and U = U2 U1 is the difference of
strain energies for two crack of close lengths under the
same external loading. In case of plane strain, h = 1. Further, in case of solving static problems in linear elasticity
theory by using the finite element method, the strain energy
equals one half of external force work, with a reverse sign:
1
U = RT u
2
where the n-dimensional vector R designates external forces
affecting the structure, and vector u represents their adequate displacements. The disadvantage is the need to conduct two analyses or some complicated interventions in the
program /6/.
KI =
J =
Similarly to Eq. [3] we can write:

8
J
+1

Vol. 4, br. 2 (2004), str. 6366
2
2
2
2
2
2
2
1 u2, 2 u1,1 + u1, 2 u2,1 + u3, 2 u3,1 dx + (11)
2 +2 u u + u u + u u + u u
dx1
(
)
,
,
,
,
,
,
,
,
1
1
1
2
3
1
3
2
2
1
1
2
1
1
2
1
2
2
where, under anti-plane strain:

2 (1 v )
1
1 =
, 2 =
(12)
1 2v
1 2v
as for plane and anti-plane stress, we have:
2
1+ v
1 =
, 2 =
(13)
1 v
1 v
The Eq. [11] is suitable for numerical calculation, since
variables in it are only displacement gradients. However,
these may be determined based on given displacement
fields in an element or, if these are not known, by using the
finite difference method /4/. Having in mind the great
accuracy of displacements in nodes determined with FEM
Determining stress intensity factors based on Rices J integral
KI =
(4)
64

Vol. 4, No 2 (2004), pp. 6366
assuming plane strain. However, this is impossible in certain practically relevant cases, such as a compact specimen.
Since the specimen (Fig. 1) is neither negligibly thin (plain
stress) and neither are its inner points far from the outer
surfaces (plane strain), the analysis must take into account
the three-dimensional stress field character. According to
authors who had solved this problem /8/, the stress intensity
factor calculated in symmetry plane differs by 12% from
the theoretical plane strain value.
as well as the generally irregular network, the Hudec /7/

method can be useful for determining displacement gradients. An elementary case will be taken for presentation. Let
nodes J and K lie on a contour , and L is the node closest
to J and not collinear to the above two. The displacement
gradients, e.g. u3 in point J, are approximately:
2
2
u3KJ xLJ
u3LJ xKL
u3LJ x1KJ u3KJ x1LJ

u
=
,
(14)
,
3
2
2
2
2
2
x1KJ xLJ
x1LJ xKJ
x1KJ xLJ
x1LJ xKJ
where u3KJ = u3K u3J represents the difference between displacements of points K and J. In an analogous way we form
coordinate differences x2LJ, etc. Also in the integral [11], we
will replace dx1 with x1KJ and dx2 with x2KJ, and replace the
sub-integral functions:
u3,1 =
)(
F1 = 1 u2, 2 + u11
u2, 2 u11
,
, +
(
)(
)
+ ( u3, 2 + u3,1 )( u3, 2 u3,1 )
F2 = 2 u11
, u1, 2 + u3,1u3, 2 + u2,1 ( 2 u11
, + 1u 2, 2 )
+ u1, 2 + u2,1 u1, 2 u2,1 +
(15)
with their mean values F1KJ and F2KJ on the segments of .

In this case the integral [11] is reduced to a sum along
segments of the considered contour:
J =
J , K = N 1, N
J , K =1, 2
2
( F1KJ xKJ
+ F2 KJ x1KJ )
(16)
The example of fracture mechanics parameter determination for a plate with a central symmetrical crack
Figure 1. Compact specimen for tensile testing.
THE PROBLEM OF NUMERICAL DETERMINATION OF

ELASTIC-PLASTIC BOUNDARY OF CRACKED BODY
We assume that the plate is in a state of plane stress. The

adopted mesh is shown in Fig. 1 of the First Paper. The
same figure also shows very good agreement of stress
results with the known analytical solution by Vestergard.
Further, in Fig. 3 of the First Paper the iso-lines of Mises
stresses are shown:
= 12 + 22 1 2
There is no basic difference between the elastic-plastic

analysis of objects with or without cracks if we use the
finite element method. If we are interested in the elasticplastic boundary, it is most usual to consider the problem as
non-linear since the strains of common structural materials
are very small, except in the negligible near crack tip area.
Anyhow, even such a simplified scheme comprises the
solving procedure which could include non-linear problems
as well. For this purpose the tangent method (NewtonRaphson) and its different modifications are used in most
cases. However, the complex procedure of calculating the
so called stiffness tangent matrix requires application of
special programmes for non-linear problems which are currently available to a limited number of users of the finite
element method.
Anyhow, it should be mentioned that programmes for
linear analysis by the finite element method can be applied
in a simple way on solving non-linear problems, especially
materially non-linear. The secant method is used in this
case (Fig. 2). The procedure is reduced to determining
secant modulus from the curve after successive iterations. With these modules we calculate new stiffness
matrixes. Experience from the author of this chapter and his
co-workers suggest that after 4-5 iterations, technically
sufficient accuracy can be achieved. From the mathematical
point of view, this convergence is linear, while for the
tangent method it is square. However, the advantage lies in
(17)
where 1 and 2 are main stresses near the crack tip.

Figure 4 in the First Paper shows a deformed configuration near the crack tip, magnified for its better presentation.
Figure 5 in the First Paper shows the graphical determination of stress intensity factors based on displacements
and stresses that were calculated using the standard package
module SMS, developed in the Aeronautical Technical Institute (VTI) in arkovo, and graphically presented in Figs. 3
and 4 in the First Paper.
Table 1 in the First Paper shows results of stress intensity factor determination using various procedures. Obviously for the adopted mesh, all of the described procedures
are accurate enough.
Certain advantage should be given to the numerical calculation using the J integral and based on Eq. [16], whose
programming gives us entirely automatic results.
Three-dimensional analysis of fracture mechanics problems
In every real fracture mechanics application, the stress
state is usually three-dimensional. Anyhow, in most cases
we can approximate it either by plane stress state or by
Vol. 4, br. 2 (2004), str. 6366
65

Vol. 4, No 2 (2004), pp. 6366
uncoupled and coupled thermal and dynamic fracture

mechanics problems /10,11/.
the fact that there is no accumulation of consecutive iteration inaccuracies and that the stiffness matrixes are generally better conditioned compared to the tangent method.
REFERENCES
1. Pian, T.H.H., Crack elements, Proceedings World Cong. FEM
Struct. Mech., Robinson ans Ass., Woodlands, Wimborne,
Dorset, England, 1975.
2. Blackburn, W.S., Calculation of stress intensity factors at
crack tips using special finite elements, The Mathematics of
Finite Elements and Applications, Academic Press, 1973.
3. Borsoum, R.S., Discussion of paper by R.D. Henshell and K.G.
Shaw, Int. J. Num. Meth. Engng. 10, 235-237, 1976.
4. Berkovi, M., Determination of stress intenisity factors using
finite elements method, in the monograph Introduction to Fracture Mechanics and Fracture-Safe Design, Institute GOA,
TMF, Belgrade, 1980, pp. 107-124.
5. Jari, J., Rices contour J integral, in the monograph Introduction to Fracture Mechanics and Fracture-Safe Design, Institute
GOA, TMF, Belgrade, 1980, pp. 69-86.
6. Hellen, T.K., On the method of virtual crack extensions, Int. J.
Num. Meth. Engng, Vol. 9, 187-207, 1975.
7. Standard Method of Test for Plane-Strain Fracture Toughness
of Metallic Materials, ASTM E399-74.
8. Hudec, M., Comparison of irregular mesh method and finite
elements method for plane elasticity problem, Software and
Hardware in Structural Analysis and Computer-Aided-Design,
VTI and MF, 1980.
9. Barsoum, R.S., On the use of isoparametric finite elements in
linear fracture mechanics, Int. J. Num. Meth. Engng, 10, 2537, 1976.
10. Hellen, T.K., Methods for computing contour integrals, PostYield Fracture Mechanics, Applied Science Publishers, London, 1979.
11. Berkovi, M., On the nonlinear transient analysis of the coupled-thermomechanical phenomena, Computers and Structures,
10, 195-202.
Figure 2. Methods for solving non-linear problems.
Once, after a completed iterative procedure, the isostress lines are drawn (see the First Paper, Fig. 3); the line
that corresponds to the yield stress for the considered
material is the elastic-plastic boundary.
THE PROSPECTIVES OF NUMERICAL ANALYSIS OF
FRACTURE MECHANICS PROBLEMS
Some of the most difficult fracture mechanics problems,
such as determining fracture mechanics parameter critical
values have not been considered in this chapter. This problem is in close connection with crack tip blunting during
deformation (First Paper, Fig. 4). Because of this blunting,
the crack tip stresses have finite values, and the next crack
growth occurs only after stress values exceed the ultimate
tensile strength for the material. This problem is, apparently, not only materially non-linear, but also geometrically
non-linear. However, also in this case we can apply the
iterative procedure from the previous section with modifications of geometric data after each iteration, but always
starting with referent configurations, and by using nonlinear strain-displacement relations. Also, there are no
principle obstacles for the numerical analysis to include

Vol. 4, br. 2 (2004), str. 6366
66

Vol. 4, No 2 (2004), pp. 6366

Numerical Fracture Mechanics

Uploaded by

Copyright:

Available Formats

Numerical Fracture Mechanics

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numerical Fracture Mechanics

Uploaded by

Copyright:

Available Formats

Second Paper

This is Chapter 5 of the invited lecture Problems in Fracture Mechanics by

NUMERICAL METHODS IN FRACTURE MECHANICS

elastic /1,2/ and elastic-perfectly plastic /3/ problems,

APPLICATION OF FINITE ELEMENT METHOD IN

Determining stress intensity factors based on given stress

Considering Eqs. A1, A3 and A7, we can conclude that

STRUCTURAL INTEGRITY AND LIFE

Numerical Methods in Fracture Mechanics

The problem is reduced to calculating the J integral,

Determining stress intensity factors based on given displacement values

Similarly to Eq. [3] we can write:

INTEGRITET I VEK KONSTRUKCIJA

1 u2, 2 u1,1 + u1, 2 u2,1 + u3, 2 u3,1 dx + (11)

where, under anti-plane strain:

Determining stress intensity factors based on Rices J integral

STRUCTURAL INTEGRITY AND LIFE

Numerical Methods in Fracture Mechanics

as well as the generally irregular network, the Hudec /7/

u3LJ x1KJ u3KJ x1LJ

with their mean values F1KJ and F2KJ on the segments of .

Figure 1. Compact specimen for tensile testing.

THE PROBLEM OF NUMERICAL DETERMINATION OF

We assume that the plate is in a state of plane stress. The

There is no basic difference between the elastic-plastic

where 1 and 2 are main stresses near the crack tip.

STRUCTURAL INTEGRITY AND LIFE

Numerical Methods in Fracture Mechanics

uncoupled and coupled thermal and dynamic fracture

Figure 2. Methods for solving non-linear problems.

INTEGRITET I VEK KONSTRUKCIJA

STRUCTURAL INTEGRITY AND LIFE

You might also like