Computers and Geotechnics 82 (2017) 157–172

Contents lists available at ScienceDirect

Computers and Geotechnics

journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Mixed mode fracture analysis of semi-circular bend (SCB) specimen:

A numerical study based on extended finite element method
Yousheng Xie ⇑, Ping Cao, Jin Jin, Min Wang
School of Resources and Safety Engineering, Central South University, Changsha 410083, China

a r t i c l e i n f o a b s t r a c t

Article history: Previous studies are mainly concentrated on the use of the semi-circular bend (SCB) specimen for deter-
Received 4 July 2016 mining the entire mixed-mode I-II fracture toughness of rock, while less attention has been paid to its
Received in revised form 22 September mixed-mode fracture process. In this situation, this study investigated mixed-mode fracture behavior
2016
of the SCB specimen using the extended finite element method (X-FEM). The crack growth trajectory,
Accepted 13 October 2016
crack initiation angle and onset of fracture were discussed in detail. This paper is expected to provide
a better understanding of mixed-mode fracture process of the SCB specimen occurring during fracture
initiation and propagation.
Keywords:
SCB specimen
Ó 2016 Elsevier Ltd. All rights reserved.
Mixed-mode fracture
Crack growth trajectory
Crack initiation angle
Onset of fracture

1. Introduction varied crack angle and external loading [1]. Therefore, the study
of mixed-mode fracture of rocks under complex loading states is
Flaws and cracks which inevitably exist in rock mass decrease an important subject for rock engineering design.
significantly the load bearing capacity of rock because of the con- Several theoretical models and experimental methods have
centration of stress at the vicinity of their tips [1]. A commonly been proposed by researchers for studying mixed-mode rock frac-
used fracture parameter called fracture toughness has been pro- ture. Among available mixed-mode fracture criteria, the minimum
posed to describe the critical states of stresses or energy in the strain energy density criterion (Smin ) [5], the maximum tangential
vicinity of the crack tip and to study crack initiation and propaga- stress criterion (MTS, rmax ) [6] and the maximum energy release
tion. This parameter usually arises in the discussions of rock burst- rate criterion (Gmax ) [7] have provided a commonly accepted theo-
ing, rock fragmentation, hydraulic fracturing, tunneling, rock retical foundation for predicting the path of crack propagation and
excavation and rock cutting processes, etc. Therefore, the determi- direction of propagation in rock under mixed-mode loading and
nation of fracture toughness for rocks plays a major role in rock have been widely used by researchers for the numerical modeling
engineering. and experimental investigations. With the exception of three crite-
Fracture in rock mass is usually in a complicated strain or stress ria mentioned above, the cohesive zone model (CZM) [8] and the
environment [2]. Pure Mode I deformation takes place for in-plane modified maximum stress criterion [9] in recent years have been
loading when the crack surfaces open without any sliding, while suggested. In addition, the results of some investigations carried
pure mode II deformation occurs when surfaces slide relatively out by rock or geological researchers show the results predicted
along the crack line [3]. Mode I fracture is frequently encountered by the latter two criteria are more coincident with the observations
failure mode of rocks against fracture, and thus extensive experi- in rock fracture tests, compared to that predicted with three classic
mental and theoretical studies have been performed in the past criteria [10,11].
for studying the mode I fracture toughness [4]. However, in At the same time, various experimental techniques to deter-
practical rock engineering cracks and flaws are more likely to be mine mixed fracture toughness develop continually together with
subjected to various mixed-mode (i.e. combinations of mode I those mixed-mode fracture criteria. To better serve the engineer-
and mode II) loading rather than pure I or pure II loading, due to ing practice and achieve consistency in the laboratory measure-
ments, the International Society for Rock Mechanics (ISRM) has
suggested four methods to measure the mode I fracture toughness
⇑ Corresponding author.
of rocks: the chevron-notched bend (CB) specimen, the short rod
E-mail address: yousheng_xie@163.com (Y. Xie).

http://dx.doi.org/10.1016/j.compgeo.2016.10.012
0266-352X/Ó 2016 Elsevier Ltd. All rights reserved.
158 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

(SR) [12] specimen, the cracked chevron-notched Brazilian disc method (GFEM) [42,43]. For the sake of simplicity, the contents
(CCNBD) [13] specimen and the semi-circular bend (SCB) [14] presented in this section are only related to the basic concept
specimen. While the punch-through-shear with confining pressure and theory of X-FEM, which can be found almost in all X-FEM
experiment [15] is the only ISRM suggested method for the deter- related documents and references. More knowledge about other
mination of Mode II fracture toughness. The centrally cracked cir- aspects of X-FEM, including its discretization of the governing
cular disc (CCCD) [16,17] specimen and the edge cracked semi- equations, can be found in the literature [47–49].
circular bend (SCB) [18,19] specimen are two frequently used disc
type specimens for the determination of mixed fracture toughness 2.1. The Heaviside function for X-FEM
of brittle material, such as rocks. In addition, in order to overcome
the shortcomings of using the classical SCB specimen, two The basic idea of X-FEM is to decompose the displacement field
improved semi-circular bend [1,20] specimens have been pro- into a continuous part ucont and a discontinuous part udisc :
posed. There is reason to pay particular attention to the SCB spec-
imen. The major advantages in using the specimen to determine uh ðxÞ ¼ ucont ðxÞ þ udisc ð1Þ
the mixed-mode fracture resistance of rock are that it can be easily where x is defined as the position vector. The first term on the right-
obtained from rock cores, has a simple geometry and only requires hand side of Eq. (1) is obtained by the conventional FEM, while the
common loading configuration. Therefore, the SCB specimen is second term is the enrichment approximation which considers the
more cost-effective, reliable and versatile compared to other spec- presence of a strong discontinuity like a crack or a flaw [48]. In
imens [21]. order to calculate the displacement of the point x located in the
In recent years, several numerical methods capable of modeling domain containing a crack, Eq. (1) is often rewritten [44,47,48] as:
crack initiation and propagation have been presented with the
development of advanced numerical techniques and constitutive X
n X
m
uh ðxÞ ¼ Ni ðxÞui þ Nj ðxÞwðxÞaj ð2Þ
modeling. The displacement discontinuity method (DDM) [22– i¼1 j¼1
24], the numerical manifold method (NMM) [25–33], the discon-
tinuous deformation analysis (DDA) [34–36], the mesh free meth- where both N i ðxÞ and N j ðxÞ are the standard FEM shape functions, ui
ods (EFG, SHP) [37], the discrete element method (DEM) [38–40], is the standard nodal displacement, aj is the added set of degrees of
the general particle dynamics (GPD) [41] and the generalized finite freedom to the standard FEM, wðxÞ is the enrichment function, n is
element method (GFEM) [42,43] are among available numerical the set of all nodal points of the domain and m is the set of nodes of
methods. Furthermore, the extended finite element method (X- the elements located on the discontinuous boundary. Choosing a
FEM) [44,45] based on partition of unity (PU) method [46] has also suitable enrichment function is a major step for modeling strong
attracted considerable attention. Modeling strong discontinuities, discontinuities. The Heaviside function (HðxÞ) has been deemed to
like cracks and flaws, with the conventional finite element method be one of the most effective enrichment functions for modeling a
(FEM) requires that the geometric of the newly produced crack strong discontinuity in X-FEM. There are two types of Heaviside
coincides with the finite element mesh, possibly resulting in functions [48], including
unpredictable mesh refinement near the crack tip when a growing 
0 if 8k < 0
crack is modeled. The PU method in the X-FEM has a good capacity HðxÞ ¼ ð3Þ
1 if 8k > 0
to describe jumps or discontinuities in the approximation field
without additional remeshing. Consequently, the X-FEM not only and
retains the advantages of the conventional FEM but also overcomes 
the shortcoming of the conventional FEM in modeling
1 if 8k < 0
HðxÞ ¼ ð4Þ
discontinuities. 1 if 8k > 0
In this paper, mixed-mode fracture behavior of the SCB speci- where 8k is defined as the signed distance. Fig. 1 shows a simple
men is investigated using the X-FEM. The path of crack propaga- one-dimensional representation of the two Heaviside functions.
tion, crack initiation angle and the onset of fracture during the Eq. (3) is the common Heaviside step function (also known as
specimen fracture process are analyzed. Comparisons of the simu- the original Heaviside function). The step function is discontinuous
lated results and curves of several mixed-mode fracture criteria are at the interface (k ¼ 0) and its application on a quadrilateral ele-
given. A better understanding of mixed-mode fracture processes ment may result in a discontinuous field [47]. In order to avoid
occurring during fracture initiation and propagation is presented. numerical problems such as instabilities in the numerical solution,
some alternative smoothed functions [66,67] have been proposed.
While Eq. (4) is the Heaviside sign function. Fig. 2 illustrates the
2. Brief introductions to the X-FEM and criterion for crack effect of Heaviside sign function in the approximation field.
initiation
2.2. Blending elements in the X-FEM
In recent years, the X-FEM originally proposed by Belytschko
and Black [44] has emerged as a very flexible method for dealing
Elements used for modeling an entire domain are categorized
with linear and nonlinear fracture mechanics problems. In the tra-
into three groups of standard, enriched and blending elements
ditional formulation of the finite element method (FEM), modeling
according to their different effects in the X-FEM. The first category
discontinuities, such as cracks and flaws, requires that the mesh
is standard finite element, while the second category is fully
conforms to the geometric discontinuities [47]. Consequently, in
order to capture the asymptotic singular fields adequately, consid-
erable mesh refinement is needed near the crack tip, possibly 0 =0 (crack) 0 0 =0 (crack) 0
resulting in a cumbersome calculation especially for modeling a
growing crack. The X-FEM alleviates this shortcoming by using Node 1 Node 2
the local partition of unity techniques to describe a jump, referred
Node 1 Node 2 Node 3 Node 4 Node 3 Node 4
to as the strong discontinuity, in the displacement field. The tech-
nique has been also widely employed in other numerical modeling
of fracture problems, such as the generalized finite element Fig. 1. Two types of Heaviside function HðkÞ.
 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172 159

" #
0 =0 (crack) 0 0 =0 (crack) 0 X
4 X
4 X
4 X
4
l
uh ðxÞ ¼ N i ðxÞui þ Nj ðxÞHðxÞaj þ Nk ðxÞ F l ðxÞbk ð7Þ
i¼1 j¼1 k¼1 l¼1
Node 1 Node 2 Node 3 Node 4 Node 1 Node 2 Node 3 Node 4
l
where bk are the DOF of the node enriched with the crack-tip func-
tions, F l ðxÞ are the crack-tip functions and have the following form
[44,47] :
½F l ðxÞ; l ¼ 1; 2; 3; 4
 
H( 3)= 1 H( 3)= 1 pffiffiffi h pffiffiffi h pffiffiffi h pffiffiffi h
¼ r sin ; r cos ; r sin sin h; r cos sin h ð8Þ
H( 2)= -1 H( 2)= -1 2 2 2 2
where r, h are local polar co-ordinates defined at the crack tip and h
is defined as the crack propagation angle with respect to the tan-
gent to the crack tip, as shown in Fig. 4. The second term on the
right-hand side of Eq. (7) is used for nodes whose shape function
N 2(x)H( ) N 2(x)[H( )-H( 2)] support is cut by the crack interior, while the third term is valid
N 3(x)H( ) N 3(x)[H( )-H( 3)] only for the crack tip.
Fig. 2. Principle of XFEM for strong discontinuities in 1D [47].
2.3. Criterion for crack initiation

governed by the X-FEM approximations and all its nodes are In this study, the cohesive crack model first introduced by Moës
enriched with extra degrees of freedom (DOF). There are internal and Belytschko [45] is used for modeling crack initiation. Analyti-
discontinuities between the standard and enriched elements due cal solutions for cracks or flaws based on the LEFM reveal that
to different element structures. In the X-FEM, blending elements there are infinite strain and stress fields at the crack tip. From a
(also known as partially enriched element) that contain at least physical point of view, however, no material in nature can with-
one enriched node are employed for treating incompatibilities of stand such an infinite stress [48]. Therefore, LEFM is only usable
the solution and interior discontinuities. Fig. 3 shows three differ- when the size of the fracture process zone (FPZ) at the crack tip
ent types of finite elements in the X-FEM. is less than that of the crack and the specimen. It prevents LEFM
Element A has one enriched node and three standard nodes, from being used widely. Conversely, the theory of fracture
while element B has three enriched nodes and one standard node. mechanics shows that the material has to undergo a progressive
The approximation of the displacement field for elements A and B nonlinear behavior in the vicinity of the crack. Nonlinear behavior
can be written respectively as: eliminating the stress singularity at the crack tip can improve
effectively properties of the material against fracture and has been
X
4 X
1
introduced into the FPZ of the cohesive crack model.
uh ðxÞ ¼ Ni ðxÞui þ Nj ðxÞHðxÞaj ð5Þ
The fracture process in the cohesive crack model can be divided
i¼1 j¼1
into four stages in terms of traction-separation relation: (1) elastic,
and (2) crack initiation, (3) crack damage (also called progressive fail-
ure) and (4) complete failure. Schematics of the FPZ and its four
X
4 X
3
uh ðxÞ ¼ Ni ðxÞui þ Nj ðxÞHðxÞaj ð6Þ stages are shown in Figs. 5 and 6.
i¼1 j¼1 Various types of traction-separation equations for the cohesive
crack model have been proposed for defining the traction-
Element C is a special element whose enrichment functions typ- separation behavior. In this paper, a triangular cohesive crack
ically consist of crack-tip functions that capture the singularity model, as shown in Fig. 7, used frequently for brittle fracture is
around the crack tip and a discontinuous function that represents used to predict crack initiation and propagation in the SCB speci-
the jump in displacement across the crack surfaces. Its approxima- men under pure mode I or mixed-mode loading. It is seen that with
tion of the displacement field can be expressed as: increasing separation displacement (dsep ), the traction load (TðdÞ)
first increases linearly, then decreases after the maximum value
C
is reached, and finally decays to zero. In this model the cohesive
fracture energy (U) defined by the area under the traction-
separation curve can be expressed [51] as follows:

A
r
=0

(a) (b)
0 Crack tip
Standard node Heaviside enriched node Crack tip enriched node
rack )
Standard element Blending element Enriched element =0 (c 0

Fig. 3. Standard, blending and enriched elements in the X-FEM: (a) a weak
discontinuity, (b) a strong discontinuity [48]. Fig. 4. A polar coordinate system at the crack tip.
160 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

Material crack tip Virtual crack tip ð1  m2 ÞK 2I

Un  GI ¼ ð11aÞ
E

ð1  m2 ÞK 2II
Ut  GII ¼ ð11bÞ
E

where E and m are the modulus of elasticity and the Poisson’s ratio,
respectively. K I and K II are the stress intensity factors of mode I and
II. It is found from Eqs. (10) and (11) that a length of the fracture
process zone fails completely and no cohesive resistance exists
δ along fracture surface when the cohesive fracture energy attains
the critical energy release rate (GC ).
When the separation displacement rises to the critical separa-
tion displacement, the traction load reaches the ultimate cohesive
strength and then the onset of fracture occurs. In the present paper,
Crack Facture process zone the maximum principal stress criterion is adopted to describe the
ultimate cohesive strength as:
To  
hrmax i
F¼ ð12Þ
oad T n;t
nl
ctio
T ra where T n;t is the tensile or shear strength, hrmax i is the maximum
principal stress. It is worth noting that tensile stresses are consid-
0 o
ered to be positive throughout this paper. Some details about Eq.
(12) require special attention: hrmax i is equal to 0 when rmax 6 0,
Fig. 5. FPZ in the cohesive crack model.
and hrmax i is equal to rmax when rmax > 0. In other words, no crack
will be resulted from the pure compressive stress.
Z dsep
1 o
Un ¼ T n ðdÞdd ¼ T df ð10aÞ
0 2 n 3. Calculation of crack parameters for the SCB specimen
Z d sep
1 o Fig. 8 shows a schematic diagram of the semi-circular bend
Ut ¼ T t ðdÞdd ¼ T df ð10bÞ
0 2 t (SCB) specimen which has been widely used by several researchers
in the past to investigate mixed-mode fracture in brittle materials
where T on and T ot are the normal and tangential cohesive strengths [52–54]. The SCB specimen containing an angled crack of length a
and are generally considered to be the ultimate tensile and shear is a semi-circular disc of radius R and subjected to three-point
strengths of materials [47], respectively. df is the post damage initi- bending. The specimen is loaded by the vertical load P and is
ation relative displacement at completely failure. The cohesive frac- placed on two bottom supports separated by a distance 2S. Differ-
ture energy (U) and the cohesive strength (T o ) are frequently ent combinations of mode I and mode II in the SCB specimen are
utilized due to the difficulty in determining the separation displace- obtained by producing specimens with varying crack angles. When
ment (dsep ). The cohesive fracture energy for brittle materials is the crack inclination angle a ¼ 0 the specimen is subjected to pure
approximately equal to its LEFM energy release rate (G) [51]. There- mode I. By increasing the angle a, the stress intensity factor (SIF) of
fore, the cohesive fracture energy (U) for plane strain may be mode II is increased, while for mode I it gives an opposite result.
expressed as: For a large value of a the SIF of mode I may be negative because

(a) (b)

(c) (d)

Fig. 6. Fracture process in the cohesive crack model: (a) elastic stage, (b) crack initiation, (c) crack damage, (b) complete failure [50].
 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172 161

Elastic limit
Elastic limit Tt Damage occurs
Tn Damage occurs max 1.0
max 1.0
Failure
0.5
G
0.5 Failure
f c
c f
G G
0.5
0.0
c f
1.0 -T t
max

Fig. 7. Triangular constitutive model for the cohesive crack model: (a) normal component, (b) tangential component [51].

P P

R R
a
a

2S 2S

(a) (b)
Fig. 8. Semi-circular bend (SCB) specimen: (a) under pure mode I loading, (b) under mixed-mode loading.

the crack faces experience a closing rather than an opening mode. results reported by Lim et al. [3], shorter crack lengths in the SCB
Pure mode II occurs for a specific angle a. specimen are preferred for fracture testing because the SIF is not very
Prior to using the SCB specimen, several crack parameters such sensitive to any variation in specimen geometry. Therefore, the val-
as K I , K II and T-stress should be known. The critical SIFs and T- ues of a=R = 0.3, 0.4 and 0.5, and S=R = 0.72 are selected for conduct-
stress for the SCB specimen can be written [3,21] as: ing numerical investigations and the crack inclination angles are set
to a = 0°, 10°, 20°, 30°, 40°, 50° and 60°.
Pcr pffiffiffiffiffiffi
K If ¼ paY I ða=R; S=R; aÞ ð13aÞ Varying finite element models of the SCB specimen are analyzed
2Rt using a commercial general-purpose finite element package ABA-
QUS for calculating Y I , Y II and T  . A typical mesh pattern generated
Pcr pffiffiffiffiffiffi
K IIf ¼ paY II ða=R; S=R; aÞ ð13bÞ for simulating the SCB specimen is shown in Fig. 9. In those models
2Rt the following parameters are applied: R = 25 mm, P cr ¼ 1 kN, t = 1
(under plane stress conditions), 2S = 36 mm, Elastic Modulus
Pcr 
T¼ T ða=R; S=R; aÞ ð13cÞ E = 75.5 GPa and Poisson’s ratio m = 0.3. A total number of 6600
2Rt plane stress quadratic elements are used in these simulations. It
where P cr is the fracture load of the tested specimen, t is the specimen is noted that in each model an embedded seam with duplicate
thickness, K If and K IIf are the SIFs of modes I and II corresponding to overlapping nodes that can separate during an analysis is intro-
the onset of specimen fracture. Normalized SIFs Y I and Y II of modes I duced to define the initial crack. The SIFs and T-stress are obtained
and II are functions of a=R, S=R and a. Normalized T-stress (T  ) for the directly from ABAQUS. The values of Y I , Y II and T  calculated from
SCB specimen is also a function of a=R, S=R and a. Lim et al. [3] different models are shown in Fig. 10.
extended previous numerical work [18] to cover a wide range of pos- It is seen from Fig. 10 that by increasing the crack inclination
sible specimen geometries of experimental interest. According to the angle, the normalized SIF (Y I ) for the mode I decreases monoto-

Fig. 9. Typical finite element mesh used for simulating SCB specimen with a=R = 0.5 and a = 40°: (a) before loading, (b) after loading.
162 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

6 a/R=0.3 1.5
a/R=0.4
5 a/R=0.5 1.2

Normalised SIF YII

Normalised SIF YI
(S/R=0.72)
4 0.9
3
0.6 a/R=0.3
2 a/R=0.4
0.3 a/R=0.5
1 (S/R=0.72)
0.0
0
0 10 20 30 40 50 60 0 10 20 30 40 50 60
Crack inclination angle (degree) Crack inclination angle (degree)
(a) (b)

4 a/R=0.3
a/R=0.4
Normalized T-stress (T*)

a/R=0.5
3 (S/R=0.72)

-1
0 10 20 30 40 50 60
Crack inclination angle (degree)
(c)
Fig. 10. Variations of three parameters for S=R = 0.72: (a) Y I . (b) Y II . (c) T  .

nously and the normalized T-stress (T  ) increases monotonously, influence on the constraint forces like the vertical load P and two
while the normalized SIF (Y II ) for the mode II decreases after an ini- bottom support loads. Further, by defining the section one can
tial increase. It is noted that the sign of several normalized T-stress obtain the same constraint force of the 2D model as that of the
values is negative. In the light of Cotterell and Rice’s [55] reports, 3D model with a thickness of 20 mm. A schematic diagram of the
both the sign and the magnitude of the T-stress exert influence on SCB specimen subjected to three-point bending is illustrated in
the stability of crack propagation. A crack subjected to Mode I load- Fig. 11.
ing propagates steadily along its extension line when the sign of its Simulated specimens are described following the sequence
T-stress is negative, while its propagation may become unstable as ‘‘Ratio of the initial flaw length to the radius of SCB specimen-Flaw
the T-stress has a positive value. inclination angle”. For example, ‘‘0.5-30” represents a simulated
SCB specimen with a=R = 0.5 and a = 30°. Following the naming
method described above, a total of 21 2D models are listed in Table 2.
4. Model setup

A total number of 21 two-dimensional SCB specimen configura- 5. Result and discussion

tions are simulated. The 2D SCB specimen has the same sample
geometry than the one used to calculate fracture toughness param- 5.1. Mixed-mode fracture criterion
eters. The mechanical parameters for simulating crack propagation
are listed in Table 1. All numerical experiments are performed Several criteria have been proposed to investigate mixed-mode
under load displacement control at a displacement-rate of I-II brittle crack initiation angles under varying loading condition
0.01 mm/step. The static friction coefficient is set equal to 0.3
between the SCB specimen and the fixed bottom supports. The
Displacement loading
numerical simulations on crack propagation are performed using
ABAQUS. Note that we enter a value of 20 mm for the section plane
stress/strain thickness when creating and editing the section. For a
two-dimensional model, the setting will have no impact on the
stress and strain changes within the model, but has an important Initial crack

25m
Table 1 m
Mechanical parameters of specimen [56].

Parameter Value Units

a
Young’s modulus (E) 75.5 GPa
Poisson’s ration (m) 0.3 /
Tensile strength (rt ) 16.1 MPa 36mm
Fracture toughness (K I ) 0.933 MPa m1/2
Fracture toughness (K II ) 1.05 MPa m1/2
Fig. 11. SCB specimen subjected to three-point bending.
 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172 163

Table 2
A total of 21 2D models.

a = 00° a = 10° a = 20° a = 30° a = 40° a = 50° a = 60°

a=R = 0.3 0.3-00 0.3-10 0.3-20 0.3-30 0.3-40 0.3-50 0.3-60
a=R = 0.4 0.4-00 0.4-10 0.4-20 0.4-30 0.4-40 0.4-50 0.4-60
a=R = 0.5 0.5-00 0.5-10 0.5-20 0.5-30 0.5-40 0.5-50 0.5-60

in various brittle materials like rock [19], polycrystalline graphite Newly produced crack Crack initiation angle
[57] and ceramics [58]. The maximum energy release rate (Gmax )
criterion [7], the minimum strain energy density (Smin ) criterion
[5] and the maximum tangential stress (MTS) criterion [6] are
among the most commonly used mixed-mode fracture criteria.
Of the three criteria, the MTS criterion is the simplest and thus is Flaw inclination angle o

used more frequently by researchers for predicting the crack initi-

ation angle. All criteria mentioned depend on the stress field exist-
ing just before the onset of crack propagation [59]. Defining r and h
as shown in Fig. 12, the near crack tip stress field for homogeneous
and isotropic linear elastic material in polar co-ordinates is defined
Main crack
as [60]:
     Fig. 13. A definition for the crack initiation angles.
1 h h 3 h
rrr ¼ pffiffiffiffiffiffiffiffiffi cos K I 1 þ sin2 þ K II sin h  2 tan
2pr 2 2 2 2
þ T cos2 h þ Oðr 1=2 Þ ð14aÞ
However, some researchers [52,57,58] found that the crack ini-
  tiation angles predicted with the conventional MTS criterion were
1 h h 3
rhh ¼ pffiffiffiffiffiffiffiffiffi cos K I cos2  K II sin h þ T cos2 h not in good agreement with those of experimental observations.
2pr 2 2 2
The principal reasons for this seem to be that only the singular
þ Oðr 1=2 Þ ð14bÞ stress term in Eq. (14b) is used to determine the crack initiation
angle and that the effect of non-singular stress term (i.e. T-stress)
1 h is neglected. In this situation, Smith et al. [9] first introduced the
srh ¼ pffiffiffiffiffiffiffiffiffi cos ½K I sin h þ K II ð3 cos h  1Þ  T cos2 h
2pr 2 generalized MTS (GMTS) criterion, also known as the modified
þ Oðr 1=2 Þ ð14cÞ MTS criterion, which considered the effects of T-stress on the orig-
inal basis. In terms of the GMTS criterion, the direction of fracture
where K I , K II and T denote the SIFs of Mode I/II and T-stress, respec- initiation (ho ) is determined [9] from:
tively. The conventional MTS criterion proposed for brittle materials
by Erdogan and Sih [6] postulates that: (1) crack propagation initi- @ rhh @ 2 rhh
¼ 0; <0 ð17Þ
ates radially from the crack tip in a plane perpendicular to the direc- @h h¼ho @h2
tion (ho ) of greatest tension (i.e. the tangential stress (rhh ) reach
and
maximum or the shear stress (srh ) is zero), (2) crack propagation
occurs when the maximum tangential stress (rhhmax ) reaches a crit- 16T pffiffiffiffiffiffiffiffiffiffi ho
K I sin ho þ K II ð3 cos ho  1Þ  2pr c cos ho sin ¼ 0 ð18Þ
ical value rhhcrit which is considered to be a material constant. 3 2
According to this criterion, by letting srh = 0, ignoring the higher
where rc is a radius of fracture process zone (also called a critical
order terms in r and only making use of the singular terms of stress
distance) and is also considered to be a material constant. The size
in Eq. (14b), the crack initiation angle (ho ) shown in Fig. 13 is given
of r c for brittle materials can be determined in terms of Taylor
[6] by:
et al.’s [61] suggestion by the following formula:
ho  2
cos ½K I sin ho þ K II ð3 cos ho  1Þ ¼ 0 ð15Þ 1 K Ic
2 rc ¼ ð19Þ
2p rt
and
By substituting the tensile strength (16.1 MPa) and mode I frac-
K II
sin ho þ ð3 cos ho  1Þ ¼ 0 ð16Þ ture toughness (0.933 MPa m1/2) listed in Table 1 into Eq. (19), rc is
KI determined to be 0.534 mm. Eqs. (13a)–(13c) are substituted into
Eq. (18) and then it can be rewritten in terms of Y I , Y II and T  as:
rffiffiffiffiffiffiffi
16T  2r c ho
Y I sin ho þ Y II ð3 cos ho  1Þ  cos ho sin ¼ 0 ð20Þ
3 a 2
rr
The results of mixed-mode fracture toughness are commonly
r
y shown in the normalized form of K II =K IIc and K I =K Ic [9]. According
to the GMTS criterion and Ayatollahi et al.’s [52] study, the ratios
K I =K Ic and K II =K IIc for the SCB specimen subjected to mixed-mode
r Crack tip Crack loading may be written in terms of Y II , Y II and T  as:
"   rffiffiffiffiffiffiffi  #1
x KI ho 2 ho 3 Y II 2rc T 2
¼ cos cos  sin ho þ sin ho ð21Þ
K Ic 2 2 2 YI a YI

Fig. 12. Stress field in polar co-ordinate system. and

164 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

"   rffiffiffiffiffiffiffi  #1

K II ho Y I ho 3 2rc T 2
Edges where the next crack tip can possibly emerge
¼ cos cos2  sin ho þ sin ho ð22Þ
K Ic 2 Y II 2 2 a Y II

The corresponding ratios of K I =K Ic and K II =K Ic are determined by Centroid

substituting the calculated crack initiation angle (ho ) determined Mesh point
from Eq. (20) and the values of Y II , Y II and T  into Eqs. (21) and (22).

5.2. Path of crack propagation

Crack propagation in engineering components has been a topic

of great interest to researchers for many years. Various methods Crack Crack tip
have been proposed to predict the path of crack propagation. Maiti
and Prasad [62] used the criteria of zero shear stress, minimum
strain energy density and maximum tangential stress to predict
the unstable crack trajectory from the stress field existing just Fig. 15. A diagram of prediction of the path of crack growth using X-FEM.
before the onset of extension. Sumi [63] developed a computer
program for the numerical prediction of curved crack growth paths
under proportional loading conditions. Oliver et al. [64] presented
two new numerical techniques (i.e. the crack-path-field and the Specimen 0.3-00
strain injection concept) to study cracks in fracture mechanics or Specimen 0.4-00
slip-lines in soil mechanics. Ayatollahi et al. [10] and Aliha et al. Specimen 0.5-00
[53] applied an incremental crack growth method to predict the At the moment of
fracture path in the SCB and CCCD specimens. Attention has been crack propagation
paid to the incremental crack growth method because its numeri-
cal prediction technology of crack propagation is similar to that
used in X-FEM. For brevity, only the salient features of this method
are described. First four parameters of K I , K II , T and rc are calcu-
lated for the pre-existing flaw using finite element analysis, then
the direction of crack propagation for the next increment is deter-
mined by using the GMTS criterion and four calculated parameters,
and then the crack extends a small distance along the calculated
direction and the same step is repeated for the next increment till
the whole crack propagation is completed. A flow diagram shown
in Fig. 14 describes the extending process of a crack. Fig. 16. Numerically simulated force–displacement curves of the SCB specimen
There are some differences in predicting the fracture path subjected to pure mode I loading.

between the incremental crack growth method and X-FEM. In

the X-FEM method, the stress or strain calculated at the element
centroid ahead of the crack tip is used for judging whether it meets force is nonlinearly increased with load point displacement till
the crack initiation criterion and for determining the crack growth the maximum load Pmax is reached and then decreased steeply.
direction, as shown in Fig. 15. Special attentions, however, should be paid to some differences
The X-FEM is one of the few numerical methods capable of pre- in the fracture process of the specimen between the results
dicting the fracture path, but some limitations should be noted: (1) observed experimentally and those obtained from the numerical
the crack tip cannot stay inside an element and thus stay only on simulation. In laboratory tests, no visual crack occurs on the sur-
its edges, (2) an enriched element is penetrated by only one crack, face of the tested specimen except invisible sub-surface crack
(3) a crack is not allowed to turn more than 90° in one increment before the maximum load is reached. The specimen fractures sud-
during an analysis. denly due to the release of the elastic strain energy when the max-
Three typical curves of force-load point displacement are shown imum load is reached, and then is split in two. However, results of
in Fig. 16. This figure shows at a constant displacement-rate the numerical simulation in this paper show a crack initiates when

Fig. 14. A flow diagram of prediction of the path of crack growth using the incremental crack growth method.
 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172 165

tensile stress in the vicinity of the pre-existing flaw tip attains ten- Numerical results in those figures also show that as the SCB
sile strength of the simulated specimen and has already propa- specimen is subjected to pure mode I loading, a crack emanates
gated before the maximum load is attained. from the tip of the flaw and propagates exactly straight along the
The fractured SCB specimens under varying combinations of extended line of the flaw till the specimen fails completely. With
modes I and II are shown in Figs. 17–19. The fracture process of the increasing flaw inclination angle (a), mode II is becoming a
a simulated specimen appears to be divided into three stages dominant effect. In this circumstance, the crack initiates at an
according to locations of crack propagation within the specimen: angle to the original flaw line around the tip of a flaw, and propa-
crack initiation stage, crack propagation stage and complete failure gates along a curvilinear path that eventually aligns with the cen-
stage. In the study, only last stage is shown. It is seen from those terline of the SCB specimen. Furthermore, it is found that crack
figures that only a single crack occurs in a tested specimen. Unlike propagating trajectories for specimens with a length ratio of
in a rock plate containing an inclined flaw under compression, no a=R = 0.3 observed in this paper are in good agreement with previ-
quasi coplanar or oblique secondary cracks are observed. It may ous experimental observations [52] shown in Fig. 20, although
be concluded that the use of the SCB specimen for determining varying S=R are employed.
mode I or mixed-mode fracture toughness may have more advan- It is noteworthy that not every crack initiates at the flaw tip.
tages in terms of stability and reliability compared to the cracked The majority of cracks initiate at the flaw tip, while some cracks
chevron-notched Brazilian disc (CCNBD) specimen or the cracked initiate behind the flaw tip. A similar observation was found in
straight through Brazilian disc (CSTBD) specimen in which two Lim et al.’s [19] works on a water-saturated synthetic mudstone.
cracks nucleate around the tips of a pre-existing flaw. All cracks with a length ratio of a=R = 0.3 emanate from their flaw

(a) 0.3-00 (b) 0.3-10

(c) 0.3-20 (d) 0.3-30

(e) 0.3-40 (f) 0.3-50

(g) 0.3-60
Fig. 17. Simulation of path of crack propagation for SCB specimens with a=R = 0.3 and S=R = 0.72.
166 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

(a) 0.4-00 (b) 0.4-10

(c) 0.4-20 (d) 0.4-30

(e) 0.4-40 (f) 0.4-50

(g) 0.4-60
Fig. 18. Simulation of path of crack propagation for SCB specimens with a=R = 0.4 and S=R = 0.72.

tips. While fractures of the specimens 0.4-50, 0.4-60, 0.5-40, 0.5-50 It is seen from Fig. 21 that there are some discrepancies
and 0.5-60 tend to initiate behind theirs flaw tips with the increas- between the simulated angles and the conventional MTS curve,
ing length ratio and flaw inclination angle. Furthermore, it is also while a relatively good agreement is shown between the simulated
found that with the increasing flaw inclination angle, a distance angles and the GMTS curve when the flaw inclination angle (a) is
becomes longer between the location of fracture initiation and low. Compared to the MTS criterion which only uses the singular
the flaw tip. It may be concluded that the location of fracture ini- stress term and ignores the effects of the non-singular stress term
tiation is resulted from the joint effects of the length and inclina- (i.e. T-stress) in a calculation, the GMTS criterion overestimates the
tion angle of a flaw. Some possible explanations for abnormal crack initiation angles (ho ) specially when the flaw inclination
initial fracture point will be discussed in detail in the subsequent angle (a) is large, whereas compared to the simulation results, it
sections of this paper. underestimates the crack initiation angles (ho ).
However, when the flaw inclination angles (a) exceed a certain
5.3. Prediction of crack initiation angle value, for example, being 50 degrees for a=R = 0.4 and being 40
degrees for a=R = 0.5, both curves cannot accurately predict the
A comparison of the results predicted by the conventional MTS direction of crack initiation. As mentioned earlier, fractures of the
criterion and GMTS criterion with the simulated results is shown in specimens 0.4-50, 0.4-60, 0.5-40, 0.5-50 and 0.5-60 tend to initiate
Fig. 21 where the x axes can be represented the mixity parameter behind theirs flaw tips. Their directions of crack initiation are very
M e [9] as: close to the direction perpendicular to the pre-existing flaw. It is
  worth noting that there is no contradiction between the third lim-
2 K II
Me ¼ tan1 ð23Þ itation of X-FEM and those simulated angles of approximate 90
p KI
 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172 167

(a) 0.5-00 (b) 0.5-10

(c) 0.5-20 (d) 0.5-30

(e) 0.5-40 (f) 0.5-50

(g) 0.5-60
Fig. 19. Simulation of path of crack propagation for SCB specimens with a=R = 0.5 and S=R = 0.72.

degrees, because a pre-existing flaw does not belong to the new The maximum principal stresses at the time when the crack
crack produced by X-FEM. starts to propagate near the crack tip for specimens 0.4-50, 0.4-60,
0.5-40, 0.5-50, 0.5-60, 0.3-00, 0.4-00 and 0.5-00 are shown in
Fig. 22. It is seen that the value of tensile stress is about 16.1 MPa
5.4. Onset of fracture
that is very close to the tensile strength of the specimen. This illus-
trates that the distribution region of tensile stress changes gradually
According to the second assumption of the MTS criterion [6],
with increasing the flaw inclination angle. For specimens 0.3-00,
brittle fracture takes place when
0.4-00 and 0.5-00 subjected to pure mode I loading, the shown stress
rhhmax ¼ rhhcrit ð24Þ contour is located symmetrically in the front of the flaw. While
under mixed loading condition the region is gradually away from
The parameter (rhhcrit ) is often considered to be the tensile strength the crack tip and moves to the upper surface of the flaw. Conse-
(rt ) for brittle materials like rock, graphite and ceramics. Several quently, fracture naturally occurs when tensile stress in the region
experimental studies have already pointed out that fracture of brit- attains the ultimate tensile strength (rt ). This seems to be a possible
tle material frequently occurs in the plane where the maximum explanation for why some fractures initiate behind the crack tip. It is
tensile stress is reached. Ayatollahi and Torabi [65] recently showed worth mentioning that Lim et al. [19] gave a similar explanation for
that the final fracture for a component will occur as the bonds that this phenomenon. Also, they held that an unusual crack initiation
hold atoms together is broken by tensile or shear loading. In addi- point in the experiment with the SCB specimen would not affect sig-
tion, similar viewpoints appeared in the literature [60] that a mate- nificantly the fracture resistance to be measured.
rial fractures when sufficient stress caused by the macro- tensile or The fracture resistance of the SCB specimen is commonly calcu-
macro-shear stress is applied at the atomic level to break the bonds. lated using the peak load (P max ) recorded during testing and the
168 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

Simulated path
Simulated path
Experimental path Experimental path

(a) α =0 (b) α =10

Simulated path
Experimental path
Simulated path
Experimental path

(c) α =20 (d) α =30

Experimental path

Simulated path
Simulated path
Experimental path

(e) α =40 (f) α =50

Fig. 20. Comparison of the simulated fracture trajectory with the experimentally observed path of crack growth [52].

obtained from pure mode I test and corresponding result are

Numerical result (a/R=0.3)
shown in Table 3.
Numerical result (a/R=0.4)
It is seen that the relative deviation between the mode I fracture
Numerical result (a/R=0.5)
toughness calculated by using the peak load and the initial flaw
length and that of the sample is between 16% and 20.6%. Clearly,
a large scatter in the fracture toughness data will inevitably occur
in this study, if the conventional calculation method for the frac-
ture resistance continues to be used. Some remedial measures for
this should be taken.
Considering the fracture resistance experiment with the SCB
MTS specimen in practice, when the peak load is reached the crack ini-
GMTS (a/R=0.3) tiates simultaneously, and then a fracture propagates rapidly till
GMTS (a/R=0.4) the specimen completely fails. An attempt is made to calculate
GMTS (a/R=0.5) the fracture resistance using the initial flaw length (a) and the
applied load (Pinitial ) at the onset of fracture, instead of the peak
load (P max ). The applied load (P initial ) is extracted from the recorded
curve of load-displacement and then is used to calculate parame-
ters related to fracture toughness. The applied load (P initial ) and cor-
Fig. 21. Comparison of the results predicted by the MTS and GMTS criteria with the responding calculated results are listed in Table 4. It is found that
simulated angles. the average mode I fracture resistance of the specimen is measured
to be about 0.949 MPa m1/2, compared to the value of
0.933 MPa m1/2 listed in Table 1. The results differ only by about
initial flaw length (a). However, the recorded curve of load- 1.7%, which indicates that the calculated result using the applied
displacement shows there is already a length of newly produced load (P initial ) is in good agreement with pure mode I fracture tough-
crack in the specimen before the maximum load is reached, as ness (KIC) of the sample.
was pointed out in part 5.1 of this paper. A question arises of A comparison of the numerical results for varying combinations
whether it is appropriate to continue using the peak load (Pmax ) of modes I and II in terms of K I =K Ic versus K II =K IIc and the GMTS cri-
and the initial flaw length (a) to calculate the fracture resistance terion curves is shown in Fig. 23, where KIC is set equal to
of the SCB specimen. An inspection of variations of the normalized 0.933 MPa m1/2. It is shown that numerical results are not com-
SIFs shown in Fig. 10 shows that the pure mode II fracture tough- pletely consistent with theoretical predictions by the GMTS crite-
ness (K IIC ) cannot be obtained from the existing results. An alterna- rion, particularly when flaw inclination angle (a) increases.
tive method to solve the question above is to compare the mode I It is remarkable that three classical mixed-mode fracture crite-
fracture toughness calculated by using the conventional calcula- ria mentioned above have provided a stable theoretical basis for
tion method with that listed in Table 1. The peak load (Pmax ) predicting the direction of crack propagation and the onset of brit-
 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172 169

(a) 0.3-00 (b) 0.4-00 (c) 0.5-00

(d) 0.4-50 (e) 0.4-50 (f) 0.5-40

(g) 0.5-50 (h) 0.5-60

Fig. 22. Distribution region of tensile stress near the crack tip.

Table 3 was completed. The calibration process was achieved by iteratively

Peak load (Pmax ) obtained from pure mode I test and corresponding result. changing the cohesive zone model parameters including Young’s
Specimen 0.3-00 0.4-00 0.5-00 modulus (E), cohesive fracture energy (U) and cohesive strength
Peak load (N) 1706 1288 992 (T o ), until the FE model results fitted the experimental results. In
KIf (MPa m1/2) 1.102 1.082 1.125 this paper, one of the reasons for the premature crack growth
Deviation from KIC (%) +18.1 +16.0 +20.6 within the SCB specimen may be that the parameter calibration
process is not carried out.
Damage evolution of the cohesive zone model can be character-
tle fracture. Fig. 24 shows a comparison between the simulated ized by either specifying cohesive fracture energy (U) or the post
results and three criterion curves. It is seen that among the three damage initiation relative displacement at completely failure (df ).
criteria, the conventional MTS criterion is closest to the simulated The energy criterion is applied in this paper. The applied crack ini-
results. This should be a reasonable result because the crack initi- tiation criterion and cohesive fracture energy (U) jointly influence
ation criterion employed in this study is the maximum principal the peak load value and its location. The fracture properties listed
stress criterion, as described earlier in the previous chapter. in Table 1 applied in calculating the cohesive fracture energy (U)
The X-FEM has already been used successfully to investigate the were determined using the CCCD specimen manufactured from
cracking behavior in the SCB specimen made of asphalt mixtures marble [56]. Clearly, heterogeneities and natural fractures within
[68–72]. Those works showed that X-FEM model simulations were the rock material were inevitably involved in the determination
in consistent with SCB test results after a model calibration process of the fracture properties. However, the SCB specimen used in this
170 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

Table 4
Applied loads (Pinitial ) and corresponding calculated results.

Specimen code a (mm) a (°) P initial (N) K If (MPa m1/2) K IIf (MPa m1/2) Initiation point
0.3-00 7.5 0 1412 0.912 0 Tip
0.3-10 7.5 10 1437 0.895 0.110 Tip
0.3-20 7.5 20 1474 0.824 0.204 Tip
0.3-30 7.5 30 1607 0.754 0.282 Tip
0.3-40 7.5 40 1836 0.681 0.344 Tip
0.3-50 7.5 50 2303 0.641 0.408 Tip
0.3-60 7.5 60 2849 0.569 0.437 Tip
0.4-00 10 0 1127 0.946 0 Tip
0.4-10 10 10 1151 0.919 0.123 Tip
0.4-20 10 20 1216 0.839 0.225 Tip
0.4-30 10 30 1380 0.755 0.308 Tip
0.4-40 10 40 1819 0.727 0.407 Tip
0.4-50 10 50 2273 0.618 0.447 Behind
0.4-60 10 60 3027 0.535 0.469 Behind
0.5-00 12.5 0 873 0.990 0 Tip
0.5-10 12.5 10 905 0.960 0.135 Tip
0.5-20 12.5 20 986 0.865 0.242 Tip
0.5-30 12.5 30 1237 0.805 0.344 Tip
0.5-40 12.5 40 1621 0.701 0.427 Behind
0.5-50 12.5 50 2348 0.585 0.511 Behind
0.5-60 12.5 60 3470 0.413 0.537 Behind

0.8 6. Summary and conclusions

GMTS (a/R=0.3)
GMTS (a/R=0.4) In this paper, an advanced numerical method, called extended
0.6 GMTS (a/R=0.5) finite element method (X-FEM), is used to study mixed-mode frac-
ture process of the SCB specimen. This paper has carried on the
detailed analysis to the crack propagating trajectory, crack initia-
KII /KIC

0.4 tion angle and onset of fracture. Several mixed-mode fracture cri-
teria are used to analyze the simulated results. However, it should
be noted that there are some model limitations in this paper such
0.2 as linear elastic, isotropic, two-dimensional, etc. Further researches
are needed to investigate the effect of the factors above on mixed-
Numerical result (a/R=0.3) mode fracture process of the SCB specimen.
Numerical result (a/R=0.4) The following conclusions may be obtained as:
0.0 Numerical result (a/R=0.5)

0.4 0.6 0.8 1.0 (1) A single crack emanates from around the tip of a pre-existing
KI /KIC flaw and propagates along a straight or curvilinear path dur-
ing fracture process of the SCB specimen.
Fig. 23. Comparison of the numerical results and the GMTS criterion curves. (2) The majority of initial fracture points are located at the crack
tip, while a crack begins to initiate behind the crack tip with
increasing the length of flaw and crack inclination angle. The
1.0 location of fracture initiation may be resulted from the joint
Numerical result (a/R=0.3) effects of the length and inclination angle of a flaw.
Numerical result (a/R=0.4) (3) A comparison shown in Fig. 21 between the simulated
0.8 Numerical result (a/R=0.5) angles of crack initiation and curves of two mixed-mode
fracture criteria shows a relatively good agreement is
0.6 obtained between simulated angles and the GMTS curve
when the flaw inclination angle (a) is low. However, crack
KII /KIC

0.4 initiation angle is close to 90 degrees when a crack initiates

behind the crack tip.
MTS (4) The curve of load-displacement recorded during the numer-
0.2 G max ical experiment shows there is already a length of newly
Smin (v=0.3, plane strain) produced crack in the SCB specimen before the maximum
0.0 load is reached. Therefore, the usual method for computing
Smin (v=0.3, plane stress)
the fracture resistance of a tested SCB specimen is aban-
0.4 0.6 0.8 1.0 doned and then the applied load at the time of the crack ini-
tiation is extracted for calculating parameters related to
KI /KIC
fracture toughness. The calculated result using the extracted
Fig. 24. Comparison of the numerical results and three classical criterion curves.
load is in good agreement with pure mode I fracture tough-
ness (KIC) of the sample.
(5) Comparisons between X-FEM simulated results with that
paper is assumed to be an isotropic linear elastic material, which is predicted by four mixed-mode fracture criteria show that
a possible reason why the result determined using the peak load the simulated results are the closest to the conventional
(Pmax ) is greater than the experimental result. MTS criterion curve.
 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172 171

(6) Possible reasons for the discrepancy between the simulated [23] Scavia C. The displacement discontinuity method for the analysis of rock
structures: a fracture mechanic. In: Aliabadi MH, editor. Fracture of
result and the fracture toughness of the sample are that nat-
rock. Boston: WIT/Computational Mechanics; 1999. p. 39–82.
ural fractures and heterogeneities within the rock material [24] Vasarhelyi B, Bobet A. Modeling of crack initiation, propagation and
are not considered and that the parameter calibration pro- coalescence in uniaxial compression. Rock Mech Rock Eng 2000;33(2):119–39.
cess is not carried out. [25] Shi GH, Goodman RE. Generalization of two-dimensional discontinuous
deformation analysis for forward modeling. Int J Numer Anal Methods
Geomech 1989;13:359–80.
[26] Shi GH. Manifold method of material analysis. In: Transactions of the 9th army
conference on applied mathematics and computing. Minneapolis, USA. p.
Acknowledgements 57–76.
[27] Tsay RJ, Chiou YJ, Chuang WL. Crack growth prediction by manifold method. J
We would like to extend our sincere gratitude to Dr. Wang Eng Mech 1999;125:884–90.
[28] Li SC, Li SC, Cheng YM. Enriched meshless manifold method for two-
YunYun and Dr. Wu QiuHong, for their instructive advice and use-
dimensional crack modeling. Theor Appl Fract Mech 2005;44(3):234–48.
ful suggestions on our paper. Special thanks to the editors and [29] Ma GW, An XM, Zhang HH, Li LX. Modeling complex crack problems with
reviewers for their hard work and valuable comments on this numerical manifold method. Int J Fract 2009;156(1):21–35.
article. [30] Zhang HH, Li LX, An XM, Ma GW. Numerical analysis of 2-D crack propagation
problems using the numerical manifold method. Eng Anal Bound Elem
2010;34(1):41–50.
References [31] Zhang GX, Zhao Y, Peng XC. Simulation of topping failureof rock slope by
numerical manifold method. Int J Comput Method 2010;7:167–89.
[32] Wu ZJ, Wong LNY. Frictional crack initiation and propagation analysis using
[1] Ayatollahi MR, Aliha MRM, Saghafi H. An improved semi-circular bend
the numerical manifold method. Comput Geotech 2012;39:38–53.
specimen for investigating mixed mode brittle fracture. Eng Fract Mech
[33] Wu ZJ, Wong LNY. Modeling cracking behavior of rock mass containing
2011;78(1):110–23.
inclusions using the enriched numerical manifold method. Eng Geol
[2] Jing L, Hudson JA. Numerical methods in rock mechanics. Int J Rock Mech Min
2013;162:1–13.
Sci 2002;39(4):409–27.
[34] Gu J, Zhao ZY. Considerations of the discontinuous deformation analysis on
[3] Ayatollahi MR, Aliha MRM. Wide range data for crack tip parameters in two
wave propagation. Int J Numer Anal Methods Geomech 2009;33(12):1449–65.
disc-type specimens under mixed mode loading. Comp Mater Sci 2007;38
[35] Ning YJ, Yang J, An XM, Ma GW. Modelling rock fracturing and blast induced
(4):660–70.
rock mass failure via advanced discretisation within the discontinuous
[4] Dai F, Wei MD, Xu NW, Zhao T, Xu Y. Numerical investigation of the
deformation analysis framework. Comput Geotech 2010;38(1):40–9.
progressive fracture mechanisms of four ISRM-suggested specimens for
[36] Ning YJ, An XM, Ma GW. Footwall slope stability analysis with the numerical
determining the mode I fracture toughness of rocks. Comput Geotech
manifold method. Int J Rock Mech Min Sci 2011;48:964–75.
2015;69:424–41.
[37] Li HQ, Wong LNY. Numerical study on coalescence of pre-existing flaw pairs in
[5] Sih GC. Strain-energy-density factor applied to mixed mode crack problems.
rock-like material. Rock Mech Rock Eng 2014;47(6):2087–105.
Int J Fract 1974;10(3):305–21.
[38] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies.
[6] Erdogan F, Sih GC. On the crack extension in plates under plane loading and
Geotechnique 1979;29:47–65.
transverse shear. J Basic Eng Trans ASME 1963;85(4):519–25.
[39] Lee HW, Jeon SW. An experimental and numerical study of fracture
[7] Hussain MA, Pu SL, Underwood J. Strain energy release rate for a crack under
coalescence in pre-cracked specimens under uniaxial compression. Int J
combined mode I and mode II. Fracture analysis ASTM STP
Solids Struct 2011;48:979–99.
560. Philadelphia: American Society for Testing and Materials; 1974 [p. 22–8].
[40] Zhang XP, Wong LNY. Cracking processes in rock-like material containing a
[8] Elices M, Guinea GV, Gomez J, Planas J. The cohesive zone model: advantages,
single flaw under uniaxial compression: a numerical study based on parallel
limitations and challenges. Eng Fract Mech 2002;69(2):137–63.
bonded-particle model approach. Rock Mech Rock Eng 2012;45:711–37.
[9] Smith DJ, Ayatollahi MR, Pavier MJ. The role of T-stress in brittle fracture for
[41] Zhou XP, Bi J, Qian QH. Numerical simulation of crack growth and coalescence
linear elastic materials under mixed-mode loading. Fatigue Fract Eng Mater
in rock-like materials containing multiple pre-existing flaws. Rock Mech Rock
2001;24(2):137–50.
Eng 2015;48(3):1097–114.
[10] Ayatollahi MR, Aliha MRM. Mixed mode fracture in soda lime glass analyzed
[42] Melenk JM, Babuška I. The partition of unity finite element method: basic
by using the generalized MTS criterion. Int J Solids Struct 2009;46(2):311–21.
theory and applications. Comput Method Appl Mech 1996;139(1):289–314.
[11] Ayatollahi MR, Aliha MRM. On the use of Brazilian disc specimen for
[43] Fries TP, Belytschko T. The extended/generalized finite element method: an
calculating mixed mode I-II fracture toughness of rock materials. Eng Fract
overview of the method and its applications. Int J Numer Meth Eng 2010;84
Mech 2008;75(16):4631–41.
(3):253–304.
[12] Ouchterlony F. ISRM commission on testing methods. Suggested methods for
[44] Belytschko T, Black T. Elastic crack growth in finite elements with minimal
determining fracturetoughness of rock. Int J Rock Mech Min Sci Geomech
remeshing. Int J Numer Meth Eng 1999;45(5):601–20.
Abstr 1988;25:71–96.
[45] Moës N, Belytschko T. Extended finite element method for cohesive crack
[13] Fowell RJ. ISRM commission on testing methods. Suggested method for
growth. Eng Fract Mech 2002;69(7):813–33.
determining mode I fracture toughness using cracked chevron notched
[46] Ohtake Y, Belyaev A, Alexa M, Turk G, Seidel HP. Multi-level partition of unity
Brazilian disc (CCNBD) specimens. Int J Rock Mech Min Sci Geomech Abstr
implicits. In: ACM SIGGRAPH 2005 courses ACM. p. 173.
1995;32(1):57–64.
[47] Khoei AR. Extended finite element method: theory and applications. John
[14] Kuruppu MD, Obara Y, Ayatollahi MR, Chong KP, Funatsu T. ISRM-suggested
Wiley & Sons; 2014.
method for determining the mode I static fracture toughness using semi-
[48] Mohammadi S. Extended finite element method: for fracture analysis of
circular bend specimen. In: The ISRM suggested methods for rock
structures. John Wiley & Sons; 2008.
characterization, testing and monitoring: 2007–2014. Springer International
[49] Pommier S, Gravouil A, Moes N, Combescure A. Extended finite element
Publishing; 2013. p. 107–14.
method for crack propagation. John Wiley & Sons; 2013.
[15] Backers T, Stephansson O. ISRM suggested method for the determination of
[50] Ha K, Baek H, Park K. Convergence of fracture process zone size in cohesive
mode II fracture toughness. In: The ISRM suggested methods for rock
zone modeling. Appl Math Model 2015;39(19):5828–36.
characterization, testing and monitoring: 2007–2014. Springer International
[51] Wang JT. Investigating some technical issues on cohesive zone modeling of
Publishing; 2012. p. 45–56.
fracture. J Eng Mater – Trans ASME 2013;135(1):011003.
[16] Chen CS, Pan E, Amadei B. Fracture mechanics analysis of cracked discs of
[52] Ayatollahi MR, Aliha MRM, Hassani MM. Mixed mode brittle fracture in
anisotropic rock using the boundary element method. Int J Rock Mech Min Sci
PMMA—an experimental study using SCB specimens. Mater Sci Eng A
1998;35(2):195–218.
2006;417(1–2):348–56.
[17] Liu HY, Kous Q, Lindqvist PA, Tang CA. Numerical modelling of the
[53] Aliha MRM, Ayatollahi MR, Smith DJ, Pavier MJ. Geometry and size effects on
heterogeneous rock fracture process using various test techniques. Rock
fracture trajectory in a limestone rock under mixed mode loading. Eng Fract
Mech Rock Eng 2007;40(2):107–44.
Mech 2010;77(11):2200–12.
[18] Chong KP, Kuruppu MD. New specimens for mixed mode fracture
[54] Aliha MRM, Ayatollahi MR, Akbardoost J. Typical upper bound-lower bound
investigations of geomaterials. Eng Fract Mech 1988;30(5):701–12.
mixed mode fracture resistance envelopes for rock material. Rock Mech Rock
[19] Lim IL, Johnston IW, Choi SK, Boland JN. Fracture testing of a soft rock with
Eng 2012;45(1):65–74.
semi-circular specimens under three-point bending. Part 2-mixed-mode. Int J
[55] Cotterell B, Rice JR. Slightly curved or kinked cracks. Int J Fract 1980;16
Rock Mech Min Sci Geomech Abstr 1994;31(3):199–212.
(2):155–69.
[20] Pirmohammad S, Ayatollahi MR. Asphalt concrete resistance against fracture
[56] Awaji H, Sato S. Combined mode fracture toughness measurement by the disk
at low temperatures under different modes of loading. Cold Reg Sci Technol
test. J Eng Mater – Trans ASME 1978;100(2):175–82.
2015;110:149–59.
[57] Ayatollahi MR, Torabi AR. Tensile fracture in notched polycrystalline graphite
[21] Lim IL, Johnston IW, Choi SK. Stress intensity factors for semi-circular
specimens. Carbon 2010;48(8):2255–65.
specimens under three-point bending. Eng Fract Mech 1993;44(3):363–82.
[58] Aliha MRM, Ayatollahi MR. Analysis of fracture initiation angle in some
[22] Scavia C, Castelli M. In: Barla G, editor. Analysis of the propagation of natural
cracked ceramics using the generalized maximum tangential stress criterion.
discontinuities in rock bridges, EUROCK’96. Rotterdam: Balkema; 1996. p.
Int J Solids Struct 2012;49(13):1877–83.
445–51.
172 Y. Xie et al. / Computers and Geotechnics 82 (2017) 157–172

[59] Khan SMA, Khraisheh MK. Analysis of mixed mode crack initiation angles [67] Patzák B, Jirásek M. Process zone resolution by extended finite elements. Eng
under various loading conditions. Eng Fract Mech 2000;67(5):397–419. Fract Mech 2003;70(7):957–77.
[60] Anderson TL. Fracture mechanics: fundamentals and applications. CRC Press; [68] Lancaster IM, Khalid HA, Kougioumtzoglou IA. Extended FEM modelling of
2005. crack propagation using the semi-circular bending test. Constr Build Mater
[61] Taylor D, Merlo M, Pegley R, Cavatorta MP. The effect of stress concentrations 2013;48(11):270–7.
on the fracture strength of polymethylmethacrylate. Mater Sci Eng – A [69] Wang H, Zhang C, Yang L, You Z. Study on the rubber-modified asphalt
2004;382:288–94. mixtures’ cracking propagation using the extended finite element method.
[62] Maiti SK, Prasad K. A study on the theories of unstable crack extension for the Constr Build Mater 2013;47(Complete):223–30.
prediction of crack trajectories. Int J Solids Struct 1980;16(6):563–74. [70] Im S, Ban H, Kim YR. Characterization of mode-I and mode-II fracture
[63] Sumi Y. Computational crack path prediction. Theor Appl Fract Mech 1985;4 properties of fine aggregate matrix using a semicircular specimen geometry.
(2):149–56. Constr Build Mater 2014;52(2):413–21.
[64] Oliver J, Dias IF, Huespe AE. Crack-path field and strain-injection techniques in [71] Ban H, Im S, Kim YR. Mixed-mode fracture characterization of fine aggregate
computational modeling of propagating material failure. Comput Method Appl mixtures using semicircular bend fracture test and extended finite element
Mech 2014;274:289–348. modeling. Constr Build Mater 2015;101:721–9.
[65] Ayatollahi MR, Torabi AR. Failure assessment of notched polycrystalline [72] Hakimelahi H. Extended finite-element modelling of asphalt mixtures fracture
graphite under tensile-shear loading. Mater Sci Eng – A 2011;528 properties using the semi-circular bending test. Road Mater Pavem 2014;15
(18):5685–95. (15):153–66.
[66] Benvenuti E, Tralli A, Ventura G. A regularized XFEM model for the transition
from continuous to discontinuous displacements. Int J Numer Meth Eng
2008;74(6):911–44.