Vietnam Journal of Mechanics, VAST, Vol. 28, No. 2 (2006), pp.

83 – 93

CALCULATION OF FRACTURE MECHANIC

PARAMETERS VIA FEM FOR SOME CRACKED
PLATES UNDER DIFFERENT LOADS
Ngo Huong Nhu and Nguyen Truong Giang
Institute of Mechanics, VAST, 264 Doi Can, Hanoi, Vietnam

Abstract. This paper describes some results in analyzing cracked plates via FEM based on
the procedures in CASTEM 2000 [1]. The basic methods for computing the crack parameters
by the finite element analysis are presented. Some programs written by GIBIAN languages
to solve problems for cracked plates are given. In possible cases, the numerical results
are composed with analytical solution or testing result that gives a good agreement. The
influence of plate configurations, the crack length, the external load type on the crack
characteristic values are considered. The numerical analysis for inclined crack at angle and
in arbitrary position of plate, the crack at hole in the plate, the crack of gravity dams are
realized. The given results and programs can be applied to practical problems for controling
the brittle failure state of a structure.

1. THE MAIN METHODS FOR CALCULATION OF FRACTURE

MECHANIC PARAMETERS VIA FEM

The determination of fracture mechanics parameters as Stress Intensity Factors (SIF),

energy release rate G, J - integral Rice and Crack Opening Displacement (COD) plays an
important role in fracture analysis. By comparing these parameters with critical values
one can estimate the brittle failure state of structures. Nowaday, the crack analysis is
founded on three different following methods.
1.1 Estimation of stress intensity factors by extrapolation of displacements or
stresses
The stress intensity factors (SIF) KI can be defined as a function of r [2]:
r
2µ 2π
KI = v(θ = π),
(χ + 1) r

where v(θ) is Westergard displacement solution:

r   r  
KI 2r θ 3θ KII 2r θ 3θ
v(θ) = (2χ + 1) sin − sin + (2χ − 3) cos + cos ;
8µ π 2 2 8π π 2 2 (1.1)
3−ν
χ= .
1+ν

In FEM the function KI (r), KII (r) are defined from displacements v(θ = π) or v(θ = 0)
at the nodes nearest to the crack tip. The values KI , KII are given by extrapolation KI (r),
KII (r) for points r near to 0 (Fig. 1).
84 Ngo Huong Nhu and Nguyen Truong Giang

Note that the stress intensity factors can be also defined from Westergard stresses
solutions. The failure occurs when KI = KIC , where KIC is a measure of the critical
material toughness.
In the Castem 2000 procedure SIF is based on the above mentioned method. It allows
to calculate the stress intensity factors KI , KII from displacements (1.1) of the sides of
the crack. The values KI , KII are a mean of the displacements of three points, which are
closest to the crack tip.

Fig. 1. The stress field near the tip of a crack

1.2 Calculation of the energy release rate G
The energy release rate G in the two-dimension body with unit thickness is evaluated
by ratio of variation of energy potential and fissure length [2]:

∂Wpot Wpot (a + ∆a) − Wpot (a)

G=− =− , (1.2)
∂a ∆a
where:
Z Z
Wpot = W (ε) − Wext = w (ε) dΩ − Tiui dΣF − the energy potential,
Ω ΣF

Z Z
w (ε)dΩ − strain potential, Tiui dΣF − external force work.
Ω ΣF

The energy potential may be received from FEM formulas:

n
X
Wpot = {ui }t Ki{ui } − {ui}t {R}.
i=1

The displacement ui , deformation and stress fields are evaluated by FEM, {R} is nodal
vector and the energy Wpot is calculated for two problems with fissures a and a + ∆a.
At fracture G = Gc, Gc is material resistant to fracture, defined by combination of
critical stress and crack size for failure [ 3].
1.3 Calculation of the J integral
The stress- strain field at the tip of a crack can be characterizied by J- integral Rice.
The J - integral is identical to G for linear elastic behavior and it can be received from
 Calculation of Fracture Mechanic Parameters Via Fem for some Cracked Plates ... 85

(1.2) for a fissure with the length of a on the axis Ox [2]:

Z  
∂ui
J= w(ε)dy + Ti dΣF , (1.3)
ΣF ∂x

where ΣF is a any contour surrounding the crack tip, which goes from lower crack side
to its upper, w is the strain energy density, Ti is a traction vector in the outward normal
along ΣF , ui is a displacement and dΣF is an element of the arc length along ΣF .
All displacement, deformation and stress values for calculating J- integral are obtained
from analyzing by FEM. The numerical integration (1.3) is calculated in each side of
element by Gauss method and the integration takes place within one layer of elements
around crack tip. Procedure G - theta in Castem 2000 is realized in this basis. It let
to calculate G and J for elastic and inelastic behaviour. Operator COD let to calculate
Crack Opening Displacement from the displacement values on the crack side.
For mixed mode I and II, G and J can be received from stress intensity factors KI
and KII by [3]:
K2 + K2
G = J = I 0 II ,
E
E
where E 0 = E - for the plane stress state, E 0 = - for the plane strain state.
1 − ν2

2. NUMERICAL SIMULATION AND CALCULATING OF FRACTURE

MECHANICS PARAMETERS FOR SOME CRACKED PLATES

In this part some results of calculating fracture mechanics parameters of plates with
fissures are presented. To obtain them, some programs are written by languages Gibian
by using operators of FEM, procedures SIF and G - theta in the Castem 2000 to modelize
different positions of crack in the plates and to solve the crack plate problems with the
different geometry and type of loading. These programs can estimate the Stress Intensity
Factor KI , KII , the energy release rate G, the displacement of the crack tip U x, and Crack
Opening Displacement COD for many problems. They are followings:
• Program RUPSB for a center crack plates in a tension.
• Program RUPTIN for plates with the inclined crack at angle β from 0x for axial and
biaxial loading.
• Program RUPHOL for crack at hole in the plate under biaxial loading.
• RUPDAM for crack gravity dam under hydrostatic loading and its own weight.
The general steps for all programs are:
- Input all material, geometrical characteristics, boundary and loading conditions for
FEM analysis of plate.
- The crack points, the length of crack tips, the crack position and sides of crack will
be input data.
- Meshing the structures and its fissure. The local FE mesh must be refined in the
surrounding or front of crack tips by special technique for each problem.
- Using operators SOLV for linear problem and procedure PASAPAS for nonlinear
problem to give corresponding displacement and stress fields.
- Application of the procedures SIF or G - theta or COD requires following input data:
86 Ngo Huong Nhu and Nguyen Truong Giang

+ The field of material properties

+ The displacement field
+ The crack points, the line describing the side of the crack (mesh),
The second side of crack for loading in mix mode or inclined crack.
- Indicate the outputs for example: Stress Intensity Factors in interested point, the
energy release rate G, the Rice – integral for a plates and Crack Opening Displacement
at any point in the crack side.
So, all these codes request not only the informations for a Fine Element Method but
also the crack situation. These received codes let automatically to mesh for local zone of
crack tips and whole structures with indicated density.
2.1. The solutions for a center cracked plate in a tension
2.1.1 The polynomial expression for a finite plate with center crack under
tensile stress [3]

√ h  πa i 1   a 2  a 4 
2
K1 = σ πa sec 1 − 0.025 + 0.06 , (2.1)
2w w w
2w is a width of the plate, 2a is a length of the fissure, σ is a tensile stress.
2.1.2 The numerical solutions KI for cracked plates with different ratios a/w

Fig. 1a. Modeling of the Fig. 1b. The model of Fig. 1c. The deformation
plate a=100 mm, w=200 mm the quater crack plate of the crack plate
a=100 mm, b=200 mm a=100 mm, w=400 mm
Consider a cracked plate subjected to a tensile stress σ =10 Mpa (Fig. 1a). The width
and length of this plate are 2w and 2h with w= 200 mm; h=600 mm. The crack is located
at the middle of plate and the length of the fissure is 2a, a= 100 mm. The plate is made
from brittle material with E= 20000 MPa, ν =0.3.
The geometry shape and loading of plate are symmetrically, therefore one considers
only1/4 of the plate, which is divided for example into 248 elements and 675 nodes (for
a/w=1/2, Fig. 1b) and 313 elements, 880 nodes (for a/w=1/4, Fig. 1c). The condition
of symmetry is modelized as following: on the line licot displacements U x=0 and U y=0
on the line libas.
 Calculation of Fracture Mechanic Parameters Via Fem for some Cracked Plates ... 87

Using RUPTS to analyze other plates with different values of ratio of crack length
and width a/w. The calculated results are compared with solution (2.1) and are shown in
the Table 1. The units of Stress intensity Factors KI and KII in all tables in this work is
(MPA.mm0.5), G (MPA .mm), 1 MPa = 1 N/mm2.

Table 1. The calculated results of center cracked plates with

a=100 (mm), w= 200, 300, 400, 500 (mm), h=600 mm

a/w Ux COD KI KI Err G KI Err

(mm) (mm) (SIF) (theory) % (G-Theta) %
1/2 -0.05397 0.107372 208.60 210.25 -0.794 2.0105 210.21 0.022
1/3 -0.04831 0.097504 188.73 190.07 -0.714 1.6453 190.16 0.045
1/4 -0.04707 0.095050 183.81 184.16 -0.190 1.5606 185.20 0.561
1/5 -0.04690 0.094550 182.77 181.58 0.650 1.5430 184.16 1.395

The results in the Table 1 show that:

- The calculated results for the KI by RUPSB with both SIF and G - Theta procedures
have a good agreement with analytical results (2.1) for all symmetrical center cracked
plates in tension.
- When the width w of a plate (or loading boundary side) more increases but the crack
length isn’t change, the stress intensity factor KI , the energy release rate G, displacement
of the crack tip U x, and Crack opening displacement COD more decrease.
2.1.3 The numerical solutions K I for cracked plates with different ratios a/h
The numerical results in the Table 2 show that the height of finite plate has a litlle
influence on the KI as a/h ≤ 1/2. The solution (2.1) with all different h is: KI (theory)
= 210.25 (MPA.mm0.5) is compared with numerical solutions as following:

Table 2. The calculated results of center cracked plates with different h

a=100 (mm), w= 200 (mm), h=200, 400, 600, 800 (mm)

h (mm) h/a Ux(mm) COD(mm) KI (Cal.) Err %

200 1/2 -0.06613 0.12627 234.83 10.46
400 1/4 -0.05415 0.10751 208.77 -0.711
600 1/6 -0.05399 0.10737 208.60 -0.794
800 1/8 -0.05399 0.10737 208.60 -0.795

In next part the unsymmetrical cracked plates will be considered.

2.2. The solutions for plates with the inclined crack
Consider plates with the inclined crack subjected to uniaxial and biaxial loading. The
numerical results for some plates with different angles β and any position of crack are
received by application RUPIN. The calculated results are checked with the theoretical
results of SIF.
88 Ngo Huong Nhu and Nguyen Truong Giang

2.2.1 The analytical expression of the stress intensity factors for infinite plate
when a crack is oriented (90◦ − β ) degrees from the applied normal stress [3]
For uniaxial loading:

KI = KI(0) cos2 (β),

(2.2)
KII = KI(0) sin(β) cos(β),

where KI(0) is the Mode I Stress Intensity when β =0. For finite plate KI(0) are taken as
the same as (2.1). The maximum KII occurs at β = 45◦, where the shear stress is also at
a maximum.
For biaxial loading:

KI = KI(0)(cos2(β) + B sin2(β)),
(2.3)
KII = KI(0)(sin(β) cos(β))(1 − B),

where B= σ1 /σ2, tension stresses on the directions Oy and Ox and σ1 > σ2.
2.2.2 The numerical solutions for rectangular plates with an inclined crack
For uniaxial loading
Consider the problem for a plate subjected to tensile stress σ = 10 Mpa, a=100
mm, w=200 mm, E=20000 MPa, ν=0.3, the crack a=100 mm is inclined at an angle
(90◦ − β) = 45◦ to the direction of extension. Using RUPTI, the plate with two crack tips
is divided into 905 elements with 2409 nodes (Fig. 2). Near to the crack tip zones the
meshes have a small density. The boundary conditions are: U y=0 at the point P0 ; U x=0
at the points P4 and AP4 . With different angles β = 0; 20; 30; 45 degrees, the results are
compared with analytical solution (2.2) in the Table 3. The error is estimated as:
The error % = [( Calculated- Ktheory)/Calculated]*100%.

Table 3. The calculated SIF for inclined crack in a plate with uniaxial loading
and different angles β
β0 Ux COD KI KI Err KII KII Err
(mm) (mm) (calcul.) (theory) % (Calcul.) (Theory) %
0 -0.0532 0.1072 207.15 210.25 -1.49 - - -
20 -0.0517 0.1054 184.41 185.66 -0.67 57.935 67.573 -16.64
30 -0.0494 0.0982 158.16 157.69 0.29 80.097 91.041 -13.6
45 -0.0428 0.0809 107.15 105.13 1.89 95.632 105.13 -9.926

For biaxial Loading

Consider above cracked plate w=200 mm, h=600 mm subjected to principal stresses
σ1 = 20 MPa and σ2= 10 MPa and B= σ1 /σ2 = 2, σ1 > σ2. The calculated results
are shown in the Table 4 are compared with theoretical stress intensity factors (2.3) for
different angles.
 Calculation of Fracture Mechanic Parameters Via Fem for some Cracked Plates ... 89

Fig. 2a. The deformation of the plate Fig. 2b. The deformation of the plate
with w = 200, h=600, β=300 with w=200, h=400, β =450

Table 4. The calculated KI , KII for different angles β for biaxial loading
β0 Ux COD KI KI Err KII KII Err
(mm) (mm) (calcul.) (theory) % (calcul.) (theory) %
0 -0.0610 0.2145 414.30 420.5 -1.49 - - -
30 -0.0629 0.2045 376.03 36794 2.15 68.039 91.041 -33.6
45 -0.059 0.1714 330.71 315.38 4.63 88.50 105.13 -19.39
60 -0.0472 0.1235 279.06 262.81 5.82 83.030 91.041 -9.647
75 -0.0264 0.0648 236.21 224.33 5.02 48.433 52.563 -8.52

Remarks:
- The written code RUPTIN for the plates with inclined crack with two crack tips gives
a high precision for KI in the all loading cases. The accuracy for KII can be accepted for
angle β greater than 45 degree.
- The values of COD, U x, KI decrease when β increase. When β = 45o, the calculated
value of KII gets a maximum as the same as in the theory.
- In the uniaxial case, when β=0 the result KI = 207.15 (Table 3) is closed to KI =
208.60 (Table 1) so the written programs give a good agreement.
2.3. The plate has an interior inclined crack in arbitrary position

Fig. 3a. The mesh of plate w=400 mm, Fig. 3b. The cracked plate after
h=600 mm, a=200 mm, β=450 deformation
The RUPTIN also can solve problems for crack in arbitrary position of the plate.
Consider two cases of the cracked plate (Fig. 3) with w= 400 mm, h=600 mm, when
the length of the crack 2a is: a=200 mm and a=100 mm. The middle point of a crack
(point P 11) has arbitrary coordinates (x, y), for example x=500 mm, y=200 mm. The
angle of the crack is inclined at 450 to Ox. The material and biaxial loading characteristics
90 Ngo Huong Nhu and Nguyen Truong Giang

are the same as the previous problem. The deformation of this plate is shown in the Fig.
3b. The boundary conditions are the same as in the above case 2. The crack characteristics
of this plate are shown in Table 5, where these values are given for a one crack tip PF, for
other crack tip they can be easy received.

Table 5. The crack characteristics of the plate with inclined crack with
the different crack length or angles

a (mm), β KI KII COD Ux

a=100, β=45◦
(EL: 1400, Nodes: 3910 ) 287.32 91.029 0.27767 -0.00309
a=200, β=45◦
(EL: 678, Nodes: 1688 ) 532.09 157.89 0.46981 0.01777
a=200, β=30◦
(EL: 682, Nodes: 1700 ) 627.45 117.25 0.56232 0.00893

Note that when the length of the crack increases, all crack characteristics increase,
when the angle of the crack decreases KI , COD increase too. So the RUPTIN can be
apply for any plate with any interior inclined crack of different values a, w, h, β, x and y.
2.4. The problem for the crack at hole in the plate
The examples in this part show capabilities of the RUPHOL for problems with com-
plicated geometry. Consider a rectangular plate with 2w=250 mm and 2h=200 mm, E=
20000 MPa, ν =0.3. The plate has a hole with radius R= 30 mm and center at C(x, y)
(Fig. 4). The crack is point lies at the hole at point P11(x= 100 mm, y=0 mm). The
length of crack a = 50 mm. The crack inclined at angle β to the axis Ox.The plate is
subjected to tension stresses on opposite sides LG and LD with σ = - 10 N/mm2. The
boundary conditions are: at P0: U x =0, U y=0l; at P 2: U y=0. The influence of the hole
position, the crack point position on the crack characteristic parameters are considered
with different C(x, y), P11(x, y).
2.4.1 The hole is not posited in the middle of plate with different angles β
In the case when the center of the hole C(100, −30), P11(x= 100 mm, y=0 mm) (Fig.
4), the results by application of program RUPHOL give:

Table 6. The calculated results for the crack at hole in the plate with different angles β

β KI KII COD(mm) Ux(mm)

300 51.964 -93.994 -1.84846.10−2 7.38299 .10–2
450 100.93 -95.166 -1.95199.10−2 6.24685. 10–2
600 147.11 -78.690 -2.46177.10−2 5.60088. 10−2

The shape of the crack after deformation is shown in the Fig. 4b.
 Calculation of Fracture Mechanic Parameters Via Fem for some Cracked Plates ... 91

Fig. 4a. The mesh of plate with Fig. 4b. The deformation of the plate
crack at hole when β = 45◦ , C(100, −30) with crack at hole
2.4.2 The hole is posited in the different places
To estimate the influence of positions of the hole on the results of crack analysis, the
problems are solved for different centers of hole C(x, y) with the highest point of crack
P 11(125, 30) on the hole. The crack inclined at angle β= 45◦ to the axis Ox. The results
of crack analysis by RUPHOL are presented in the Table 7.

Table 7. The crack characteristics for the crack at hole with different center C(x, y)
C(x,y) KI KII COD(mm) Ux(mm) J or G
(100,-30) 100.93 -95.166 -1.95199.10−2 6.24685. 10–2 8.80967. 10−1
(110,-20) 98.596 -95.728 -1.73148. 10−2 6.09693. 10–2 8.6510. 10–1
(110,-10) 98.019 -97.952 -1.69857. 10–2 5.38747. 10–2 8.79161.10−1
(125,0) 105.66 -101.33 -1.19074.10−2 5.57361. 10–2 9.79797. 10−1

Thus, when the hole in the middle of plate C(125, 0), all the absolute values KI , KII ,
J become more while the values of COD and U x are less than case, when the hole is sited
in the left lower corner of the plate.
2.4.3 The influence of the crack point position

Fig. 5a. The crack shape on the hole in Fig. 5b. The crack shape on the hole in
the middle plate, P11(125, 25) the middle plate, P11(150, 0)

Consider the plate with the hole of radius R = 25 in the middle. Two positions of
crack point P 11 are investigated: the highest position (in direction y): P 11(125, 25) and
lower position P 11(150, 0). The results are shown in the Table 7, Fig. 5a and Fig. 5b.
92 Ngo Huong Nhu and Nguyen Truong Giang

Table 8. The crack characteristics for the crack at hole in

the plate with different position of P 11
P11(x, y) KI KII COD Ux
(mm) (mm) (mm)
(125, 25) 90.372 -92.990 -1.30421.10−2 6.19066. 10–2
(150, 0) 62.826 -69.336 -1.5623.10−2 8.47052. 10–2

Thus, the position of crack point on the hole has also noticeable influence on the values
KI , KII , COD and the crack shape. So, the program RUPHOL can be used to investigate
the plate with crack at hole for many different sizes and positions of the hole.
2.5. Numerical fracture analysis of a gravity dam under self-weight and hy-
draulic loading
With the aim of checking the program, the fracture calculation by RUPDAM is realized
for the experimental test 2D model in [5]. In that model the geometric scale between
prototype and model was taken to be 40. The model has a horizontal notch on the
upstream side at a quarter of the dam height; with a notch depth is 20 cm (Fig. 6a). The
hydraulic load was modeled by the force 1000 KN, which was distributed into concentrated
load at all points in the upstream side.
In the case the material has linear behavior with Young’s Modulus E= 35.700 MPa
the following cases are calculated:
- When external loads are hydraulic loads 1000 KN without self-weight: The calculated
Crack Opening Displacement is COD(calcul.) = 0.0449 mm, while the test result COD
(test) ∼
= 0.055 mm, so RUPDAM gives a acceptable error and the crack shape is shown in
Fig. 6b. The calculated Stress Intensity Factors KI =0.703596 MPA.mm0.5, KII =0.055930
MPA.mm0.5.
- When external loads are hydraulic with self-weight, the calculated values COD =
0.030638 mm, KI = 0.500374 MPA.mm0.5, KII = 0.062564 MPA.mm0.5. Note that these
values are smaller than the values in the case without self-weight and it is reasonable.

Fig. 6a. The scheme of the testing Fig. 6b. The calculated deformation
of the crack dam
 Calculation of Fracture Mechanic Parameters Via Fem for some Cracked Plates ... 93

3. CONCLUSION

In this work, the numerical simulation for crack plates with complicated geometry
under different loads is realized successfully. The crack characteristics as COD, Stress
Intensity Factors, J –integral or the energy release rate G are received for a plates with a
centered, inclined crack under uniaxial and biaxial loading and for crack at arbitrary hole
in the plate. The 2D model of crack gravity dam under hydrostatic loading and self-weight
is also considered. The influence of the plate size or crack size, of the crack point position
at hole on the crack characteristics is investigated.
The developed programs based on the Castem 2000 are checked, and compared with
analytical or testing results, they give a good agreement, and can be applied to other
practical structures.
This work is completed with the partly financial support from the Vietnam National
Council of Natural Sciences.

REFERENCES
1. C. E. A/ D. M. T/L. A. M. S., Castem 2000 Recueil D’exemples Commentes, 1992.
2. Michel Prrat, Philippe Bisch, Alain Millard vv.., Calcul des Ouvrages Generaux de
Construction, AFPC- Emploi des elements finis engineer civil, Hermes, Paris, 1997.
3. T. L. Anderson, Fracture Mechanics. Fundamentals and Applications, CRC Press.
Boca Raton Ann Arborv London Tokyo, 1992.
4. Stanley T. Rolfe, John M. Barsom, Fracture and Fatigue Control in Structures. Ap-
plication of Fracture Mechanics, Prentice- Hall, Inc., Englewood Cliffs, New Jersey,
1979.
5. Alberto Capinteri and Silvio Valente and others. Expremental and numerical fracture
modelling of a gravity dam. Proceeding of the First International Conference on Frac-
ture Mechanics of Concrete Structures (FraMCoS1), Colorado, USA, 1992 Elsevier
Applied Science.

Received October 8, 2005

Revised February 20, 2006

TÍNH TOÁN CÁC THÔNG SÔ ’ A CO. HOC PHÁ HUY

´ CU ` NG
’ BĂ
.
´ ’ . ’
PPPTHH CHO MÔ . T SÔ BAN NÚ T CHI.U CÁC DA. NG TAI KHÁC NHAU

Bài báo dã phân tı́ch các bài toán ba’n có vê´t nú.t bă `ng phu.o.ng pháp phâ ` n tu’. hũ.u
ha.n. Co so’ tı́nh toán các thông sô dă.c tru ng vêt nú t dã du o. c trı̀nh bày. Mô.t sô´ chu.o.ng
. . ´ . ´ . . .
trı̀nh dã du.o..c thiê´t lâ.p trên co. so’. phát triê’n, su’. du.ng các moduyn và toán tu’. cu’a Castem
2000, kê´t qua’ sô´ du.o..c kiê’m chú.ng vó.i nghiê.m gia’i tı́ch thê’ hiê.n dô. tin câ.y cu’a chu.o.ng
trı̀nh. A’nh hu.o’.ng cu’a kı́ch thu.ó.c ba’n, dô. dài vê´t nú.t, vi. trı́ vê´t nú.t, da.ng ta’i tro.ng dê´n
các thông sô´ SIF, J, G, COD vv dã du.o..c xem xét. Các bài toán ba’n có vê´t nú.t nghiêng,
ba’n có vê´t nú.t ta.i lô
˜ hô’ng, mă.t că ´t dâ.p có vê´t nú.t chi.u ta’i thuy’ tı̃nh dã du.o..c phân tı́ch.
Kê´t qua’ nghiên cú.u và các chu.o.ng trı̀nh nhâ.n du.o..c cho kha’ năng ú.ng du.ng trong dánh
giá phá huy’ công trı̀nh bê tông.