Rock Mechanics and Rock Engineering

https://doi.org/10.1007/s00603-020-02139-7

ORIGINAL PAPER

Fracture Assessment of Quasi‑brittle Rock Simulated by Modified

Discrete Element Method
Sohrab Gheibi1   · Rune M. Holt1

Received: 12 June 2018 / Accepted: 18 April 2020

© The Author(s) 2020

Abstract
New developments of an in-house hybrid code, named Modified Discrete Element Method (MDEM) are presented in the
paper. The new developments are on the treatment of pre-existing and propagating fractures in quasi-brittle materials. These
developments are the embedment of Linear Elastic Fracture Mechanics (LEFM) and elastic-softening crack band model
-based methodologies in the MDEM and their application in lab and reservoir scale. Using the first methodology, MDEM
can calculate stress intensity factors, K I and K II using the internal contact forces of particles. K I and K II are calculated inde-
pendent of boundary conditions and geometrical configuration with acceptable accuracy level. The methodology has been
also used in reservoir scale to study the rupture likelihood of faults and fractures due to fluid injection. This methodology
enables the code to model mode I and mode II failures and propagation direction based on the fracturing model proposed by
Rao et al. (Int J Rock Mech Min Sci 40(3): 355–375, 2003). Using the second methodology, the MDEM can model nonlinear
behavior of quasi-brittle materials including or excluding preexisting cracks based on fracture energy. A model was verified
against an experiment of a three point bend test with a notch. The numerically obtained force-crack mouth opening curve
was reasonably comparable to the experimental test. The analysis was repeated for three other mesh sizes and the results are
less mesh size dependent. Finally, it was shown that MDEM has the potential in studying fracture mechanics of quasi-brittle
materials both in lab and large-scale investigations.

Keywords  Hybrid FEM/DEM · Stress intensity factor · Strain softening · Linear elastic fracture mechanics · Quasi-brittle
material
List of Symbols H Softening modulus
a Half a crack size H e Equivalent softening modulus
aij Contribution of the deformation of jth contact K 𝛼 Stress intensity factor being 𝛼 = I, II in mode i
on the force of ith contact and ii, respectively
A Area of a cluster/element K 𝛼e Equivalent stress intensity factor in mode i and
𝐂 Elasticity matrix ii being 𝛼 = I, II , respectively
dt Time-step K Material toughness in mode i and ii being
𝛼C

E Young’s modulus 𝛼 = I, II , respectively

Et Tangential modulus K ∗ Stress intensity factor estimation function
𝐟n Contact normal force 𝐊̄ Cluster’s contacts stiffness matrix
ft Tensile strength 𝐌 Displacement–strain transformation matrix
Gf Fracture energy P Pressure
hc Crack band width Pmin Local minimum load in Brazilian test
he Element size Pmax Maximum load in Brazilian test
R Brazilian disc radius
t Brazilian disc thickness
* Sohrab Gheibi 𝐮n Relative normal displacement of contact
[email protected] xc Ligament size
1 𝛽 Crack angle with respect to the horizontal line
Department of Geoscience and Petroleum, Norwegian
University of Science and Technology (NTNU), S. P. 𝛾 f Fracture energy density
Andersens veg 15a, 7031 Trondheim, Norway 𝛆 Strain

13
Vol.:(0123456789)
 S. Gheibi, R. M. Holt

𝜀𝛼 Elastic and plastic strain being 𝛼 = e, p , calculate the toughness in rock cutting. Moon et al. (2007)
respectively developed a general approach to measure fracture toughness
𝜃 Direction of an arbitrary plane at the crack tip under random packing of non-uniform size particles. They
𝜃 𝛼C Direction of minimum tensile and maximum used the energy balance approach which is based on the
shear stress being 𝛼 = I, II , respectively equilibrium state between strain, friction, and kinetic energy
𝜆 Ratio of maximum to minimum stress as the internal and the total accumulated work done by the
𝛬 Plastic multiplier loaded boundaries as the external energies. Their second
𝜇 Friction coefficient method was collocation method which is determining the
𝜈 Poisson’s ratio coefficients in a series of stress solution based on the com-
𝛔 Stress plex stress function presented by Westergaard (1939) and
𝛷max Maximum dimensionless stress intensity factor expanded by Sanford and Berger (1990). In this method, the
result is sensitive to the number of data points. Although,
Abbreviations
PFC is widely used in modeling of rock failure and fracture
SIF Stress intensity factor
propagation, like any other method or code, it has particu-
LEFM Linear elastic fracture mechanics
lar drawbacks such as low unconfined compressive to indi-
DEM Discrete element method
rect tensile strength ratio, effect of the inherent roughness
MDEM Modified discrete element method
of interfaces forming discontinuities on frictional behav-
CMOD Crack mouth opening displacement
ior. Lisjak and Grasselli (2014) have discussed the major
drawbacks and reviewed the recent studies improving them.
Among them are the paper by Potyondy and Cundall (2004)
1 Introduction where they used cluster of particle together to obtain a more
realistic macroscopic friction values. Cho et  al. (2007)
Failure of quasi-brittle materials such as rocks is associ- showed that clumped-particle geometry lead to a more cor-
ated with the localization of strain into a finite band forming rect compressive to tensile strength ratio. Potyondy (2012)
macroscopic fracture. Accumulation of these microcracks has alternatively, formulated the flat-joint model to capture
due to their initiation, growth and coalescence imposes a the effect of clumped-particles.
nonlinear/softening behavior in the stress–strain curve of Hybrid continuum-discontinuum method has recently
rocks. In the recent decades, several researchers have devel- gained attention among the researchers. Hybrid methods,
oped numerical method (codes) to capture the complicated generally consider the domain as a continuum and as soon as
process of fracturing. Jing and Hudson (2002) have classi- the cracks are about to form, the domain becomes discontin-
fied the most commonly applied numerical methods in Rock uum in the corresponding region. ELFEN (Rockfield Soft-
Mechanics in three categories of (1) Continuum method ware Ltd. 2004) is an example of a hybrid code that is based
including the Finite Difference Method (FDM), the Finite on finite element formulation. The transition from elastic
Element Method (FEM) and the Boundary Element Method continuum to discontinuum is controlled by dissipation of
(BEM), (2) Discrete Element Methods (DEM), (3) Hybrid a definite amount of strain energy during strain localization
continuum/discrete methods. into a crack band. For details about ELFEN the readers are
The DEM was introduced by Cundall (1971) for the referred to Klerck (2000) and Profit et al. (2015). Another
analysis of rock-mechanics problems and then applied to known hybrid method is the combined finite-discrete ele-
soils by Cundall and Strack (1979). Particle Flow Code, ment method developed by Munjiza (2004) and extended by
PFC (Itasca Consulting Group Inc. 2012) as a DEM based Mahabadi et al. (2012) named Y-Geo by improving the limi-
code discretizes the domain by a number of discs in 2D and tation of Munjiza’s FDEM such as including a quasi-static
spheres in 3D which are in contact two by two transferring friction law, Mohr–Coulomb and rock joints shear failure
normal and shear forces in contact bond and also resisting criterion etc. Guo et al. (2016) have also developed a three
rolling induced by moment in parallel bond version (Poty- dimensional version of the Munjiza’s FDEM that has the
ondy and Cundall 2004). Fractures are formed by breakage ability to capture transition from continuum to discontinuum
of the bonds and their coalescence with other broken bonds. and the explicit interaction of discrete fractures. Lisjak and
Classical PFC does not use fracture mechanics theories Grasselli (2014) have given a comprehensive review on
directly, but Potyondy and Cundall (2004) showed that the some of DEM and hybrid codes.
fracture toughness can be related to the properties of paral- Hori et  al. (2005) and Alassi (2008) have separately
lel bond model if a definite size of particles (a characteristic developed methodologies by which the microscopic param-
length) are chosen with two non-dimensional constants α eters Kn (contact normal stiffness) and Ks (contact shear stiff-
and β. Marina
√ et al. (2015) obtained a cubic relationship ness) are calculated using macroscopic Young’s modulus
between 𝛽 𝛼 and K I . Huang et al. (2013) also used PFC to and Poisson’s ratio. The main objective in their works was

13
Fracture Assessment of Quasi‑brittle Rock Simulated by Modified Discrete Element Method

removing the difficult process to calibrate the PFC model to 2 Modified Discrete Element Method
find the desired condition to capture the effect of Poisson’s
ratio. The improper arrangement of particles does not pro- The Modified Discrete Element Method (MDEM) was pro-
duce a lateral deformation due to the axial loading in DEM posed by Alassi (2008) to model fracture developments and
(Hori et al. 2005). However, the lateral deformation caused fault reactivation during fluid withdrawal and injection at
by an axial deformation both globally and locally determines reservoir scale. MDEM is a hybrid code that behaves like a
the stress changes in large scale applications especially in continuum model (e.g., finite element method) before failure
oil–gas depletion or fluid injection such as C ­ O2 storage and and like a discontinuum model (e.g., discrete element method)
Enhanced Oil Recovery (EOR) (Fjaer et al. 2008). There- after failure.
fore, there is no need for the calibration procedure required Figure 1 shows a triangular element formed by connecting
in PFC. Hori et al. (2005) called their method/code FEM-β the centers of three discs which are in contact two by two. The
which was developed in the framework of FEM. In contrast, triangle element is also called a cluster.
Alassi (2008) developed his methodology in the framework The
{ constitutive
}T relationship of the normal forces,
of DEM and called it Modified Discrete Element Method 𝐟𝐧 = {fn1 fn2 fn3 }, and the normal relative displacements,
(MDEM). MDEM has been successfully applied in reser- T
𝐮𝐧 = un1 un2 un3 , of the three contacts of a cluster
voir geomechanics and hydraulic fracturing recently (Lavrov (Fig. 1) is given as Alassi (2008)
et al. 2016; Gheibi et al. 2016, 2017). Like the FEM, MDEM
suffers from mesh/particle size dependency in modeling
̄ 𝐧,
𝐟𝐧 = 𝐊𝐮 (1)
crack propagation. In the local FEM approach, mesh size
⎛ kn1 a12 a13 ⎞
dependency is treated by correcting the softening modulus where, 𝐊 ̄ = ⎜ a21 kn2 a23 ⎟.
for the element size by enforcing a mesh size independ- ⎜ ⎟
⎝ a31 a32 kn3 ⎠
ent fracture energy (Bazant and Planas 1997). However, in
It is assumed that a contact can only develop normal forces,
others, a material length scale is incorporated in the for-
therefore, the shear stiffness of contact are set to zero. Instead,
mulation of finite elements. In the nonlocal damage model
aij are introduced to the matrix of stiffness in Eq. (1). aij rep-
formulation, the stress at a point not only depends on the
resents the contribution of the deformation of jth contact on
strain at that point, but also on the strain at a neighborhood
the force of ith contact.
of that point (e.g. Bažant and Pijaudier-Cabot 1988; Guy
Stress state is related to the internal forces of the contacts.
et al. 2012). The gradient-enhanced approach incorporates
Using Gauss theorem, stresses are retrieved from the internal
the dependency on the second gradient of an invariant of
forces as following (Alassi 2008)
the plastic strain (e.g. De Borst and Mühlhaus 1992). In
the micropolar elastoplaticity, the Cosserat’s continuum 1 𝐓
𝛔= 𝐌 𝐟𝐧 , (2)
is adopted (e.g. Steinmann and Willam 1991). Consider- A
ing rate-dependency of the yield surface is also another { }T
where, 𝛔 = 𝜎xx 𝜎yy 𝜎xy  , A is the area of the cluster (tri-
approach to overcome the instability in strain localization
angle), the internal forces 𝐟𝐧 and 𝐌 is the unit normal vector
(Wang et al. 1997). Alternatively, a finite width of a localiza-
matrix defined as Alassi (2008)
tion band is reduced to zero width (or an interface) known
as a strong discontinuity. Assumed enhanced strain (Simo ⎛ I11
2 2
d1 I12 d1 I11 I12 d1 ⎞
and Rifai 1990) and the extended finite element (Moës et al. ⎜ 2 2
𝐌 = I21 d2 I22 d2 I21 I22 d2 ⎟, (3)
1999) methods are classified among the approaches used ⎜ 2 2 ⎟
⎝ I31 d3 I32 d3 I31 I32 d3 ⎠
assuming a strong discontinuity which uses classical plastic-
( ) ( )
ity theory without introducing a length scale. where Ie1 = cos 𝜃e Ie2 = sin 𝜃e and the angle 𝜃e represents
In this paper, we will discuss the new developments in the normal vector orientation of contact e inside the cluster,
MDEM by presenting two different methodologies that de is the contact length (the distance between the centers of
are based on Linear-Elastic Fracture Mechanics (LEFM) the two particles that are in contact), and e = 1, 2, 3 is the
and Elastic-Softening Fracture Mechanics to deal with the ID of contacts.
fracture problem. The application of MDEM in the lab and Strain of a cluster is calculated using (Alassi 2008)
large scale will be shown by modeling several examples.
The general formulation of MDEM will be given in Sect. 2. Fig. 1  Representation of
In Sect. 3, the application of MDEM in LEFM will be dis- MDEM element (Alassi 2008)
cussed by solving several examples. Section 4 covers the
application of MDEM in elastic softening fracture mechan-
ics of quasi-brittle materials. Finally, the conclusions will be
given in the last section.

13
 S. Gheibi, R. M. Holt

are selected based on the direction of minimum principal stress

𝐮𝐧 = 𝐌𝛆. (4)
i.e. the contacts under lower normal force. Post-peak soften-
The geometrical
{ matrix
}T 𝐌 transforms the strain of a clus- ing/hardening behavior is captured by dividing the strain (𝜀)
ter 𝛆 = 𝜀xx 𝜀yy 𝜀xy to the relative displacements 𝐮𝐧 , of into two elastic (𝜀e ) and plastic parts (𝜀p ) i.e. 𝜀 = 𝜀e + 𝜀p . The
contacts. plastic strain is calculated as Chen and Han (2007)
The stress can be related to the strain by using the mate-
𝜕F
rial conventional constitutive elastic matrix 𝐂 as 𝜀p = 𝛬 , (7)
𝜕𝜎
𝛔 = 𝐂𝛆 (5)
where Λ is the plastic multiplier defined by the plastic power
So, by combining Eqs.  (1)–(5) a relationship can be equivalence (plastic work) and is a function of the harden-
derived between the conventional elasticity matrix and inter- ing/softening modulus H (Chen and Han 2007). Since the
nal contact stiffness matrix 𝐊
̄ as following (Alassi 2008) internal calculation of the code is based on displacements
rather than strains, 𝐮p as the plastic relative displacement
1 𝐓̄
𝐂= 𝐌 𝐊𝐌. (6) is calculated using 𝜀p and 𝐮e = 𝐮 − 𝐮p is the elastic relative
A
displacement. 𝐮e is used to update the forces in Eq. (1). Frac-
The solution scheme is an explicit one like the regular tures appear if the effective plastic strain in a given cluster
discrete element method after 𝐊 ̄ is retrieved from Eq. (6). reaches a critical value. The critical plastic strain is dictated
The contact forces are updated ( 𝐟𝐧new = 𝐟𝐧 + d𝐟𝐧 ) based on by H or equivalently the fracture energy. Figure 2 repre-
the incremental change in the forces ( d𝐟𝐧 = 𝐊d𝐮 ̄ 𝐧 ) due to sents the geometrical relation between fracture energy ( Gf  ),
contact’s relative displacement. d𝐮𝐧 is calculated based on softening modulus (H), tangential modulus (Et ) , Young’s
the difference of velocities of the two particles in contact in modulus ( E  ) fracture energy density ( 𝛾 f  ), tensile strength
time step dt. ( ft ) for a bar breaks in uniaxial tension.
Newton’s second law is used to update the particle’s According to Bazant and Planas (1997) the crack band has
motion [see Cundall and Strack (1979) for more details]. a finite width ( hc ) that is a material constant depending on its
Clusters fail in tension or shear after they reach their peak maximum particle size. From numerical point of view, the size
strength value. For a failed cluster, two of the contacts fail. of the elements should not necessarily be exactly equal to hc
Discontinuities such as joints and faults can be introduced but larger element size can be used in the model. However, the
to the model. A discontinuity is a collection of failed clusters softening modulus (H) or equivalently 𝜀fc should be updated to
which has two contacts failed. This provides flexibility in guarantee the correct energy release rate independent of the
introducing preexisting discontinuities even with complex mesh size. The fracture energy Gf defines the amount of energy
geometries. The direction of the failed contact is updated dissipated per unit cross sectional area of the crack band of
to the direction of the discontinuity. This is equivalent to width hc . For an arbitrary mesh, strain-softening should neces-
the smooth-joint model in DEM introduced by Mas Ivars sarily localize into an element of dimension,he . The equivalent
et al. (2008). softening modulus ( H e ) is determined by enforcing a constant
The failure criterion is Mohr–Coulomb with a tension cut- fracture energy Gf as following
off. After peak strength of an element, the strain localizes into
hc he
a band with a finite width (Bazant and Planas 1997) and a Gf = hc 𝛾f = he 𝛾fe ⇒ 𝛾fe = 𝛾 f ⇒ H e
= H. (8)
finite amount of energy is dissipated until the continuum clus- he hc
ter turns to discontinuum. Two of the three contacts of the Therefore, the results should be less mesh size sensitive.
cluster are allowed to separate or slide. The two failed contacts To avoid a negative value of H e for very large mesh sizes,

Fig. 2  The crack band model, a

system response, b crack band
softening (Bazant and Planas
1997; Klerck 2000)

13
Fracture Assessment of Quasi‑brittle Rock Simulated by Modified Discrete Element Method

the tensile strength can also be modified to ensure constant √ �

l

fracture energy (Bazant and Planas 1997). K ∗I = 𝜋∕ 2xc Fy,i (13)

It should be noted that the examples modeled with elastic- i=1

softening in this paper have been restricted to tensile mecha-

nism and the softening after shear failure was beyond the √ �
l

goal of this paper. However, Alassi and Holt (2011) modeled K ∗II = 𝜋∕ 2xc Fx,i . (14)
i=1
softening behavior under compression with MDEM.
The main difference between MDEM and the regular K I and K II are calculated by extrapolating to xc = 0 of
discrete element method is that the material behaves like the linear approximations of K ∗I and K ∗II as a function of
a continuum before failure, where the conventional elastic xc plots, respectively. The deviation from the analytical
properties are given as input values (Eq. (6)). The mate- K I and K II values is about 1%. The method provides an
rial behaves like a regular discrete element method after the acceptable accuracy for reasonably coarser mesh.
deformation of an element reaches 𝜀fc .

2.1 Calculation of KI and KII in MDEM 2.2 Mixed Mode I–II Cracks and Brittle Fracture
Criterion
For an arbitrarily oriented crack in an isotropic body under
a mixed-mode I–II loadings in plane stress or plane strain Erdogan and Sih (1963) proposed the composite criterion
conditions. The stresses near the crack tip are of minimum circumferential tensile stress (minimum ten-
�√ sile–stress criterion). This criterion holds that a mixed
𝜎yy = K I 2𝜋x (9) mode I–II crack propagates along the corresponding direc-
tion of minimum tensile stress satisfying the following
�√ (modified after Wu et al. 2016)
𝜏xy = K II 2𝜋x. (10)
𝜕𝜎𝜃𝜃 𝜕 2 𝜎𝜃𝜃
= 0, >0 (15)
The total forces over the xc − sized ligament (Fig. 3) can 𝜕𝜃 𝜕𝜃 2
be expressed as de Morais (2007)
√
xc 1− 1 + 8w2 K II
√ 𝜃 IC
= 2 tan −1
, w= I, (16)
∫
Fy = 𝜎yy dx = K I 2xc ∕ 𝜋 (11) 4w K
0 where 𝜃 IC is the corresponding direction of minimum tensile
stress.
xc
√ The corresponding SIF of minimum tensile stress or
∫
Fx = 𝜏xy dx = K II 2xc ∕ 𝜋. (12) equivalent mode I intensity factor is given by (modified
0 after Wu et al. 2016)

And their values can be evaluated from the internal forces 1 𝜃 IC [ I ( ) ]

K Ie = cos K 1 + cos 𝜃 IC − 3K II sin 𝜃 IC . (17)
of the contacts for a cluster at its center of gravity as shown 2 2
in Fig. 3. SIFs are estimated using (de Morais 2007)
The fracture criterion for mode I crack grow/rupture is
Rao et al. (2003)

|K | ≥ K IC ,
| Ie |
| | (18)

where K IC is the fracture toughness in mode I.

To obtain the direction of maximum shear stress, the
following condition should be satisfied
𝜕𝜏𝜃r | ( )|
= 0, |𝜏𝜃r 𝜃 = 𝜃𝛼IIC | . (19)
𝜕𝜃 | |max

This leads to a cubic equation with three roots. We derive

Fig. 3  A pre-existing crack and clusters in front of the crack tip and the following solution as the corresponding direction of
calculation of forces required to approximate K I and K II (after Gheibi maximum (absolute value) shear stress (Gheibi et al. 2018)
et al. 2018)

13
 S. Gheibi, R. M. Holt

( )
(g )√ ( p ) bottom was fixed in y direction. The difference between the
−1 𝜋 w
𝜃𝛼IIC = 2 tan 2 cos + (𝛼 − 2) − + , 𝛼 = 1, 2, 3, exact solution and the numerical one is about 1%.
3 3 3 3
It is also possible to calculate K I and K II for biaxial stress
(20)
boundary condition in any ratio of the applied stress, no
where, p= −  ,
− 13 (w)2 q=7
2
−  ,
2
27
(w)3 2
3
w matter if they are compression–compression, tensile–tensile
⎛ ⎞
and tensile-compression. When the crack surface is closed
g = cos−1 ⎜− 21 � �q �3 ⎟ and 𝛼 is the index of the three roots
⎜ ⎟ (under compression) the corresponding K II is corrected by
⎝ − p3 ⎠
imposing the effect of the friction coefficient. Based on the
and only one of them maximizes the shear stress. fracture criterion proposed by Rao et al. (2003), K Ie and K IIe
Therefore, the corresponding SIF of maximum shear are calculated using Eqs. (17 and 21). Figure 5 shows K Ie ,
stress or equivalent mode II intensity factor is given by | IIe / Ie |
K and the corresponding |K
IIe
K | for a central crack with
(modified after Wu et al. 2016) | |
𝛽 inclination under two uniaxial and compressive–compres-
1 𝜃 IIC [ I ( )] sive stresses with 𝜆 = 0, 0.3 and friction coefficient
K IIe = cos K sin 𝜃 IIC + K II 3 cos 𝜃 IIC − 1 . (21)
2 2 𝜇 = 0, 0.6 . Figure 5 shows that mode I is more likely to
occur for a central crack than mode II, because either
The failure criterion for shear is given by Rao et al. (2003) | IIe / Ie | /
|K K | is lower than K IIC K IC for rocks (Al-Shayea and
| |
K | ≥ K IIC K IC
{ | IIe / Ie | /
|K Khan 2000; Rao et al. 2003; Backers and Stephansson 2012)
| |
K IIe ≥ K IIC
, (22) or K IIe is very low.

where K IIC is the fracture toughness in mode II. 3.2 Non‑Collinear Cracks

It must be noted that K I in Eqs. (17) and (21) is only
meaningful if it is negative. Otherwise it should be set to The methodology can be applied to calculate K I and K II for
zero. preexisting cracks with complex geometries such as curved
and multiple cracks in different geometries. Figure 6a repre-
sents a model of two non-collinear symmetrical cracks with
3 Linear Elastic Fracture Mechanics 2a and 𝜃 inclination with centers that are 2.2a away from
one another under uniaxial tensile stress. Figure 6 displays
3.1 Central Crack Under Uniaxial and Biaxial different solutions [by Sih (1973), Chen and Chang (1989),
Loading Guo and Ma (2011), Guo et al. (2013)] for K I for two sym-
metrical non-collinear cracks with several inclination 𝜃 . In
Figure 4 represents K I and K II numerically obtained for a the figure, the dots are the values obtained by this methodol-
central crack with a length of 2a = 1 inclined by 𝛽 varying ogy for the two tips of the crack inclined by 20° and 0° from
from 0° to 90° from the horizontal line. A vertical, constant the horizontal line.
uniaxial tensile stress was applied on the top edge and the

Fig. 4  Comparison of analytical and numerical a K I and b K II as a function of the central crack inclination 𝛽

13
Fracture Assessment of Quasi‑brittle Rock Simulated by Modified Discrete Element Method

| / |
Fig. 5  a K Ie and K IIe and b |K IIe K Ie | for a central crack under uniaxial (λ = 0) and biaxial compression (λ = 0.3) with friction coefficient μ = 0,0.6
| |

3.3 Reservoir Scale Application ratio of the models is 0.25. The elastic constants are the same
for the reservoir and the surrounding rock with E = 15 GPa
Figure 7 shows a schematic description of the reservoir and 𝜈 = 0.2 . The studied fault with length of 100 m and
model used to study of the stress intensity factor change after inclination of 60° were placed at (0, 0). The mechanical
injection of a fluid such as C­ O2. The reservoir is assumed boundary condition at the boundaries of the models is a
to be closed and the pore pressure is increased inside it to constant stress equal to the far-field stress. The maximum
mimic the fluid injection. It was assumed that the pore pres- principal stress is considered to be in the vertical direction
sure remains constant in the surrounding rock as well as in called normal faulting or extensional stress regime. The ratio
the caprock. The reservoir thickness is 200 m and the aspect of horizontal stress to vertical stress was assumed to be 0.4.

Fig. 6  a Schematic of two collinear cracks, and b several solutions for K I(after Guo et al. 2013)

13
 S. Gheibi, R. M. Holt

independent of friction coefficient. K II increases for increasing

pore pressure so that the mode II failure occurrence becomes
more likely. It is observed that for higher friction coefficient
a greater pore pressure increase is needed to have K II > 0 .
Increase of K II loses its dependency on the friction coeffi-
cient after K I becomes negative. The reason is that for K I < 0
the fault plane is open and normal force is zero so no shear
resistance.
Equations (17 and 21) were used to calculate the equivalent
stress intensity factors to be able to assess the rupture likeli-
hood and plotted vs. pressure increase in Fig. 9. In contrast
to Fig. 8, K I was restricted to be ≤ 0 in Fig. 9, this also means
that only negative K I were used in the calculations in Eqs. (17
and 21). Figure 9a indicates that the mode I rupture can occur
after 5 and 15 MPa pressure increase for the 0.6 and 1 fric-
Fig. 7  Schematic representation of the general model, which includes tion coefficient cases, respectively depending on the tough-
a fault at center ness value ( K IC ). The direction of mode I failure initiation
was calculated to be 70.5° in the fault tip. However, the mode
I rupture will probably change its direction to grow vertically
depending on the level of pressurization. Mode II failure is
more complicated because not only it depends on the K IIC
but also the ratio K IC ∕K IIC . Figure 9b indicates that mode II
is more likely for pressure lower than 30 MPa, because the
K Ie ∕K IIe is lower. The reason is that for pressure > 30 MPa
K I becomes negative (Fig. 8), therefore, K Ie has contributions
from K I and K II(Eq. 17). However, for pressure < 30 MPa, K I is
positive so it does not contribute to K Ie in Eq. (17). The mode
II rupture (if occurs) will propagate along the fault plane for
pressure < 30 MPa, and − 2.2° and − 5.8° (counterclockwise
being positive). For more comprehensive discussion on this
subject the readers are referred to Gheibi et al. (2018).

4 Elastic‑Softening Crack Band model

Fig. 8  Variation of K I and K II vs pore pressure increase 4.1 A Model Without a Preexisting Crack: Brazilian
test
The pressure inside the reservoir section was increased by
5, 10, 15, 20, 25, 35 and 40 MPa in separate models. K I and First example in this section is modeling rock Brazilian test.
K II were measured using the method described in 2.1 in both This test is performed to obtain the indirect tensile strength
initial and after injection states. The friction coefficient of of rocks. Brazilian test is also used to calculate the toughness
the fault plane assumed to be 0.6 and 1. Figure 8 shows the of rock either including (Tang 2017) or excluding (Guo et al.
variation of K for different pressure change for the fault with 1993) a preexisting crack. Guo et al. (1993) proposed a method
different friction coefficients. ­ IC from the force curve of a Brazilian test by
to calculate the K
√ �The stress �intensity values (after Wang and Xing 1999)
were normalized by 𝜎V 50𝜋 MPa m−0.5  . Therefore, the
results can be used for different initial in situ stress values. P
It is important to note that positive K I is not meaningful. K IC = √ min 𝛷max , (23)
However, the positive values were used here to show to how R×t
much degree the fault tip was under compression and how where Pmin is the local minimum load which can be obtained
it unloaded due to increase in the pressure. Figure 8 shows from the force curve, 𝛷max is the maximum dimensionless
the K I decreases for increasing pressure and it becomes nega- SIF and R and t are radius and thickness of the disc, respec-
tive after ΔP > 32 MPa. This means that after this amount of tively. 𝛷max is dependent on the angle of the arc made by
pressure increase mode I failure can occur. As expected, K I is

13
Fracture Assessment of Quasi‑brittle Rock Simulated by Modified Discrete Element Method

Fig. 9  Variation of a K Ie and b K IIe vs pore pressure increase and rupture criteria (Eqs. 17 and 22)

the loading plate and the sample. This angle is considered

to be ~ 10° for the standard Brazilian test. However, experi-
mental and numerical studies have shown that the stand-
ard Brazilian test can lead to a shear failure in the loading
region rather than tensile failure at center of the specimen
depending on the ratio of the shear to tensile strength (Fair-
hurst 1964; Gheibi et al. 2015) and deformability of the disc
(Gheibi et al. 2015). Therefore, the modeled specimen was
loaded by a plate with an arc angle of 30°. The radius of
the disc is 52 mm and thickness was assumed to be 26 mm.
𝛷max was calculated to be 0.443 (Erarslan et al. 2012). The
Young’s modulus, Poisson’s ratio and tensile strength are 10
GPa, 0.25 and 7 MPa, respectively with H varying from 10
to 100 MPa. Figure 10 shows the force-strain curve of the
model with several H. Table 1 shows Pmax ,Pmin , calculated
indirect tensile strength (σt), and the toughness value cal-
culated by Eq. (23). As it is expected, the toughness value
is greater for lower softening modulus. Lower softening Fig. 10  Brazilian test model with several softening moduli
modulus means that the cracks require more energy to pro-
ceed. Also, the estimated indirect tensile strength is higher Table 1  Data of the Brazilian test models
for lower H models. Figure 11 shows the tensile fracture H (MPa) Pmin (kN) Pmax (kN) σt (MPa) KIC
formed for different values of H. The fracture has propagated (MPa × m0.5)
a longer distance for higher H values.
100.00 13.75 20.31 8.06 1.45
4.2 A Model with a Preexisting Crack: Three Points 50.00 22.27 25.12 9.97 2.35
Bend Test 30.00 29.77 32.50 12.90 3.15
20.00 39.21 41.57 16.50 4.14
Three points bend test with a notch is a commonly used test 10.00 60.48 61.52 24.41 6.39
to determine the toughness and fracture energy of quasi-
brittle materials. To verify the capability of the code in
elastic-softening fracture mechanics a three points bend was Berea sandstone. The sandstone specimen was prepared
experiment was modeled. The experiment data was picked with length = 276 mm (span = 254 mm), height = 99.9 mm,
up from Fakhimi et al. (2017). The sample used in this study and thickness = 24 mm and a straight notch was cut with

13
 S. Gheibi, R. M. Holt

Fig. 11  Load-CMOD curve of
experimental and numerical
models with several mesh sizes

length = 15 mm and width = 1.7 mm at the center of the

beam (Fakhimi et al. 2017).
A model of the experiment was built in the code based
on the geometry given above. A 15 mm crack was applied
to the model instead of cutting the notch section in mesh
generation. The elements/clusters were equilateral triangular
with sides equal to 1 mm in front of the notch. The model
matched to the experimental data has a Young’s modulus of
10 GPa (Fakhimi et al., 2017), Poisson’s ratio of 0.25, ten-
sile strength 5.1 MPa and softening modulus, H = 40 MPa.
The bending was simulated by applying a velocity boundary
and the reaction force and Crack Mouth Opening Displace-
ment (CMOD) were recorded numerically. Similar to the
experiment, the model was unloaded at a point close to the Fig. 12  Load-CMOD curve of experimental and numerical models
with several mesh sizes

13
Fracture Assessment of Quasi‑brittle Rock Simulated by Modified Discrete Element Method

50% of the maximum force after peak (Fakhimi et al. 2017). factors, K I and K II using the contact forces of particles. It is
Figure 12 shows the Force-CMOD curve of the experiment able to be used in complex boundary condition and geomet-
and the simulation with a relatively reasonable match with rical configuration such as curved and interacting multiple
the two. Fakhimi et al. (2017) inspected the unloaded beam cracks with acceptable accuracy. The methodology has been
and a crack of about 25 mm in length above the notch was also used in reservoir scale to study the rupture likelihood of
observed with an unaided eye. Interestingly, the propagated faults and fractures. This methodology enables the code to
crack (cohesionless term was used by Fakhimi et al. 2017) model mode I and mode II failures based on the fracturing
was 22 mm in the model. The Force–deflection curve could model proposed by Rao et al. (2003).
not be verified because Fakhimi et al. (2017) did not meas- Another development was embedding elastic-softening
ure the deflection (personal communication). Three other crack band model into MDEM. The code is able to model
similar models with the element sizes of 1.5 mm, 3 mm, the nonlinear behavior of quasi-brittle materials including
and 6 mm were built to investigate the effect of mesh size and excluding preexisting cracks. A Brazilian (without a
sensitivity. Figure 12 shows that the Force-CMOD curve notch) test was modeled by the second methodology and
was reproduced with an acceptable accuracy for different effect of softening modulus was studied on obtaining the
values of mesh size. Figure 13 shows the discretization of models mode I toughness. Also, an experiment of a three
the model with 1 mm and 6 mm element sizes overlapping points bend test with a notch was modeled in the paper. The
one another. It should be noted the mesh dependency of numerically obtained force-crack mouth opening displace-
strain localization is not restricted to only the size of the ele- ment was reasonably comparable the experimental test. The
ments (Guy et al. 2012), but its geometry and its relation to model was repeated for three different meshes 1.5, 3 and
the fracture path. In this example, the geometry and loading 6 time larger than the initial mesh and the results are less
condition is such that the fracture formed is in the vertical mesh size sensitive. Finally, it was shown that MDEM has
direction, therefore, the size of the spatial discretization play the capability in studying fracture mechanics of quasi-brittle
a more important role compared to its orientation. Jirásek materials both in the lab size and large-scale investigations.
and Bauer (2012) showed that the effective width of a crack
band depends on the shape of the element and its orientation Acknowledgements  Open Access funding provided by NTNU Nor-
with respect to the overall direction of the band. wegian University of Science and Technology (incl St. Olavs Hospital
- Trondheim University Hospital). This publication has been produced
with partial support from the BIGCCS Centre, performed under the
Norwegian research program Centres for Environment-friendly Energy
5 Conclusion Research (FME). The authors acknowledge the following partners for
their contributions: Gassco, Shell, Statoil, TOTAL, ENGIE, and the
Research Council of Norway (193816/S60). Authors are also grateful
Modified Discrete Element Method (MDEM) and its new to SINTEF Petroleum Research providing the MDEM code. S.G is
developments in fracture problem were presented in the thankful of Hadi Zadeh Haghighi and Hamed Homaei Rad.
paper. A Linear Elastic Fracture Mechanics (LEFM) based
methodology was adopted to calculate stress intensity Open Access  This article is licensed under a Creative Commons Attri-
bution 4.0 International License, which permits use, sharing, adapta-
tion, distribution and reproduction in any medium or format, as long
as you give appropriate credit to the original author(s) and the source,
provide a link to the Creative Commons licence, and indicate if changes
were made. The images or other third party material in this article are
included in the article’s Creative Commons licence, unless indicated
otherwise in a credit line to the material. If material is not included in
the article’s Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will
need to obtain permission directly from the copyright holder. To view a
copy of this licence, visit http://creat​iveco​mmons​.org/licen​ses/by/4.0/.

References
Alassi HT (2008) Modeling reservoir geomechanics using discrete ele-
ment method: application to reservoir monitoring. PhD thesis,
NTNU
Alassi H, Holt R (2011) Modeling fracturing in rock using a modified
discrete element method with plasticity, in proceedings key engi-
Fig. 13  Comparison of the discretization of the model with 1- and neering materials. Trans Tech Publ 452:861–864
6-mm mesh sizes

13
 S. Gheibi, R. M. Holt

Al-Shayea N, Khan K (2000) Effects of confining pressure and tem- Guy N, Seyedi DM, Hild F (2012) A probabilistic nonlocal model for
perature on mixed-mode (I-II) fracture toughness of a limestone crack initiation and propagation in heterogeneous brittle materi-
rock formation. Int J of Rock Mech and Rock Sci 37(4):629–643 als. Int J Numer Methods Eng. https​://doi.org/10.1002/nme.3362
Backers T, Stephansson O (2012) ISRM suggested method for the Hori M, Oguni K, Sakaguchi H (2005) Proposal of FEM implemented
determination of mode II fracture toughness. Rock Mech and with particle discretization for analysis of failure phenomena. J
Rock Eng 45(6):1011–1022 Mech Phys Solids 53(3):681–703
Bažant ZP, Pijaudier-Cabot G (1988) Nonlocal continuum damage, Huang H, Lecampion B, Detournay E (2013) Discrete element mod-
localization instability and convergence. J Appl Mech Trans eling of tool-rock interaction I: rock cutting. Int J Numer Anal
ASME 55:287–293. https​://doi.org/10.1115/1.31736​74 Meth Geomech 37(13):1913–1929
Bazant ZP, Planas J (1997) Fracture and size effect in concrete and Itasca Consulting Group Inc (2012) PFC2D (particle flow code in 2
other quasibrittle materials. CRC Press, Boca Roton dimensions). Minneapolis, USA
Chen W-H, Chang C-S (1989) Analysis of two-dimensional mixed- Jing L, Hudson JA (2002) Numerical methods in rock mechanics. Int
mode crack problems by finite element alternating method. Com- J Rock Mech Min Sci 39(4):409–427
put Struct 33(6):1451–1458 Jirásek M, Bauer M (2012) Numerical aspects of the crack
Chen W-F, Han D-J (2007) Plasticity for structural engineers. J. Ross band approach. Comput Struct 110–111:60–78. https​: //doi.
Publishing, Richmond org/10.1016/J.COMPS​TRUC.2012.06.006
Cho N, Martin CD, Sego DC (2007) A clumped particle model for Klerck PA (2000) The finite element modelling of discrete fracture in
rock. Int J Rock Mech and Min Sc 44(7):997–1010 quasi-brittle materials. University of Wales Swansea, Swansea
Cundall PA (1971) A computer model for simulating progressive, large Lavrov A, Larsen I, Bauer A (2016) Numerical modelling of extended
scale movement in blocky rock systems. In: Proceedings Symp. leak-off test with a pre-existing fracture. Rock Mech Rock Eng
ISRM, Nancy, France, Proc.1 971, Volume 2, pp 129–136 49(4):1359–1368
Cundall PA, Strack ODL (1979) A discrete numerical model for granu- Lisjak A, Grasselli G (2014) A review of discrete modeling techniques
lar assemblies. Géotechnique 29(1):47–65 for fracturing processes in discontinuous rock masses. J Rock
De Borst R, Mühlhaus H-B (1992) Gradient-dependent plasticity: Mech Geotech Eng 6(4):301–314
formulation and algorithmic aspects. Int J Numer Methods Eng Mahabadi OK, Lisjak A, Munjiza A, Grasselli G (2012) Y-Geo: new
35:521–539. https​://doi.org/10.1002/nme.16203​50307​ combined finite-discrete element numerical code for geomechani-
de Morais AB (2007) Calculation of stress intensity factors by the force cal applications. Int J Geomech 12(6):676–688
method. Eng Fract Mech 74(5):739–750 Marina S, Derek I, Mohamed P, Yong S, Imo-Imo EK (2015) Simula-
Erarslan N, Liang ZZ, Williams DJ (2012) Experimental and numerical tion of the hydraulic fracturing process of fractured rocks by the
studies on determination of indirect tensile strength of rocks. Rock discrete element method. Environ Earth Sci 73(12):8451–8469
Mech Rock Eng 45(5):739–751 Mas Ivars D, Potyondy DO, Pierce M, Cundall PA (2008) The smooth-
Erdogan F, Sih GC (1963) On the crack extension in plates under plane joint contact model. In: Proceedings of WCCM8-ECCOMAS, v.
loading and transverse shear. J Basic Eng 85(4):519–525 2008, p. 8
Fairhurst C (1964) On the validity of the ‘Brazilian’ test for brittle Moës N, Dolbow J, Belytschko T (1999) A finite element method
materials. Int J Rock Mech Min Sci Geomech Abstr 1(4):535–546 for crack growth without remeshing. Int J Numer Methods Eng
Fakhimi A, Tarokh A, Labuz JF (2017) Cohesionless crack at peak load 46:131–150. https​://doi.org/10.1002/(SICI)1097-0207(19990​
in a quasi-brittle material. Eng Fract Mech 179:272–277 910)46:1%3c131​:AID-NME72​6%3e3.0.CO;2-J
Fjaer E, Holt RM, Raaen A, Risnes R, Horsrud P (2008) Petroleum Moon T, Nakagawa M, Berger J (2007) Measurement of fracture tough-
related rock mechanics. Elsevier, Amsterdam ness using the distinct element method. Int J Rock Mech Min Sci
Gheibi S, Holt RM, Lavrov A, Mas Ivars D (2015) Numerical modeling 44(3):449–456
of rock brazilian test: effects of test configuration and rock hetero- Munjiza AA (2004) The combined finite-discrete element method.
geneity. In: 49th US Rock Mechanics/Geomechanics Symposium- Wiley, Hoboken
ARMA, American Rock Mechanics Association (ARMA), p 9 Potyondy DO (2012) A flat-jointed bonded-particle material for hard
Gheibi S, Holt RM, Vilarrasa V (2016) Stress path evolution during rock. American Rock Mechanics Association, Alexandria
fluid injection into geological formations. In: Proceedings 50th Potyondy DO, Cundall PA (2004) A bonded-particle model for rock.
US Rock Mechanics/Geomechanics Symposium 2016, American Int J Rock Mech Min Sci 41(8):1329–1364
Rock Mechanics Association Profit M, Dutko M, Yu J, Cole S, Angus D, Baird A (2015) Comple-
Gheibi S, Holt RM, Vilarrasa V (2017) Effect of faults on stress path mentary hydro-mechanical coupled finite/discrete element and
evolution during reservoir pressurization. Int J Greenh Gas Con- microseismic modelling to predict hydraulic fracture propagation
trol 63:412–430 in tight shale reservoirs. Comput Part Mech 3:1–20
Gheibi S, Vilarrasa V, Holt RM (2018) Numerical analysis of mixed- Rao Q, Sun Z, Stephansson O, Li C, Stillborg B (2003) Shear fracture
mode rupture propagation of faults in reservoir-caprock system in (Mode II) of brittle rock. Int J Rock Mech Min Sci 40(3):355–375
­CO2 Storage. Int J Greenh Gas Control 71:46–61 Rockfield Software Ltd. ELFEN 2D/3D numerical modelling package
Guo Z, Ma H (2011) Solution of stress intensity factors of multiple Rockfield Software Ltd., Swansea, UK (2004)
cracks in plane elasticity with eigen COD formulation of boundary Sanford R, Berger J (1990) The numerical solution of opening mode
integral equation. J Shanghai Univ (English Ed) 15(3):173–179 finite body fracture problems using generalized Westergaard func-
Guo H, Aziz NI, Schmidt LC (1993) Rock fracture-toughness determi- tions. Eng Fract Mech 37(3):461–471
nation by the Brazilian test. Eng Geol 33(3):177–188 Sih GC (1973) Handbook of stress-intensity factors. Lehigh University,
Guo Z, Liu Y, Ma H (2013) Solution of stress intensity factors for Institute of Fracture and Solid Mechanics, Bethlehem
2-D multiple crack problems by the fast multipole boundary ele- Simo JC, Rifai MS (1990) A class of mixed assumed strain methods
ment method. Bound Elem Other Mesh Reduct Methods XXXVI and the method of incompatible modes. Int J Numer Methods Eng
56:417 29:1595–1638. https​://doi.org/10.1002/nme.16202​90802​
Guo L, Xiang J, Latham J-P, Izzuddin B (2016) A numerical investi- Steinmann P, Willam K (1991) Localization within the Framework
gation of mesh sensitivity for a new three-dimensional fracture of Micropolar Elasto-Plasticity. In: Brüller OS, Mannl V, Najar J
model within the combined finite-discrete element method. Eng (eds) Advances in Continuum Mechanics: Springer Berlin Hei-
Fract Mech 151:70–91 delberg, Berlin, Heidelberg, pp 296–313

13
Fracture Assessment of Quasi‑brittle Rock Simulated by Modified Discrete Element Method

Tang SB (2017) Stress intensity factors for a Brazilian disc with a Wu LZ, Li B, Huang RQ, Wang QZ (2016) Study on Mode I-II hybrid
central crack subjected to compression. Int J Rock Mech Min Sci fracture criteria for the stability analysis of sliding overhanging
93:38–45 rock. Eng Geol 209:187–195
Wang Q-Z, Xing L (1999) Determination of fracture toughness KIC by
using the flattened Brazilian disk specimen for rocks. Eng Fract Publisher’s Note Springer Nature remains neutral with regard to
Mech 64(2):193–201 jurisdictional claims in published maps and institutional affiliations.
Wang WM, Sluys LJ, Borst R de (1997) Viscoplasticity for instabili-
ties due to strain softening and strain-rate softening. Int J Numer
Methods Eng 40:3839–3864. https:​ //doi.org/10.1002/(SICI)1097-
0207(19971​030)40:20<3839::AID-NME24​5>3.0.CO;2-6
Westergaard H (1939) Bearing pressures and cracks. J Appl Mech
18:49–53

13