See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/30509122

Discrete fracture in quasi-brittle materials under compressive and

tensile stress states

Article  in  Computer Methods in Applied Mechanics and Engineering · July 2004

DOI: 10.1016/S0045-7825(04)00119-7 · Source: OAI

CITATIONS READS
79 168

3 authors, including:

Ewan Sellers
The University of Queensland
41 PUBLICATIONS   375 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Blast Design View project

HSBM-Hybrid Stress Blasting Model View project

All content following this page was uploaded by Ewan Sellers on 13 October 2018.

The user has requested enhancement of the downloaded file.

 Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056
www.elsevier.com/locate/cma

Discrete fracture in quasi-brittle materials under

compressive and tensile stress states
P.A. Klerck a, E.J. Sellers b, D.R.J. Owen a,*

a
Department of Civil Engineering, University of Wales, Singleton Park, Swansea, Wales SA2 8PP, United Kingdom
b
Division of Mining Technology, CSIR, Johannesburg, South Africa
Received 20 March 2003; received in revised form 3 October 2003; accepted 3 October 2003

Abstract

A method for modelling discrete fracture in geomaterials under tensile and compressive stress fields has been
developed based on a Mohr–Coulomb failure surface in compression and three independent anisotropic rotating crack
models in tension. Extension fracturing is modelled by coupling the softening of the anisotropic rotating crack failure
criterion to the compressive plastic strain evolution. Modifications were introduced into an explicit discrete element/
finite element code with an explicit Lagrangian contact algorithm to enforce non-penetration of the surfaces created
when the tensile strength is depleted. The model is applied to triaxial and plane strain tests, as well as punch tests and
borehole breakouts to show that the model is able to quantitatively predict the appropriate load–displacement response
of the system in addition to the observed evolution of discrete fracturing in situations comprising a variety of com-
pressive and tensile stress states.

2004 Elsevier B.V. All rights reserved.

1. Introduction

In the analysis of the stability of mine excavations, it is vital to be able to determine the degree of
fracturing in the rock and the response of the fractured rock mass to subsequent changes of the stress state,
for a variety of excavation shapes and sequences. The constitutive relationship for quasi-brittle materials in
compression is generally determined by performing conventional triaxial tests ðr1 ¼ r2 P r3 Þ, extension
tests ðr1 P r2 ¼ r3 Þ and true triaxial tests ðr1 P r2 P r3 Þ on small material specimens. Typical experimental
results for rock in compression are presented and discussed by Hallbauer et al. [1], Tapponnier and Brace
[2], Franklin [3], Janach [4], Mogi [5–7], Yumlu and Ozbay [8], Hoek and Bieniawski [9], Brace et al. [10],
Wang and Shrive [11] amongst others. The final failure planes observed in unconfined compressive strength
(UCS) tests or extension tests are generally parallel to the direction of maximum compressive stress,
whereas in standard triaxial tests with confinement they are generally at some specific angle. Similar
observations are made in the mining environment, where the absence of confinement in material adjacent to

*
Corresponding author. Fax: +44-1792-295-252.
E-mail address: d.r.j.owen@swansea.ac.uk (D.R.J. Owen).

0045-7825/$ - see front matter
2004 Elsevier B.V. All rights reserved.
doi:10.1016/j.cma.2003.10.015
3036 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

excavations or in mining pillars causes fracturing parallel to the maximum compression direction resulting
in so-called slabbing, as shown in Fig. 1(a) and (b). The ultimate aim of the model is to predict the
development of fracturing around an advancing stop face in a deep level mine. As shown in Fig. 1(c), the
rock mass is extremely fractured, with fractures forming in uniaxial compression at the face and in
extension as the rock moves into the excavated region. In addition, in confined regions ahead of the mining
face the rock fails in a mechanism similar to that of a conventional triaxial test with oblique failure planes
resulting in shear deformation evolve making specific angles with the maximum compressive stress direc-
tion. The model, therefore, must be able to predict the simultaneous formation of fractures due to all these
different mechanisms. It is noted that the extension test, conventional triaxial test and the uniaxial com-
pression (UCS) test only provide insight into the material response to stress states on the boundary of the
full stress domain. Unfortunately, these experiments do not truly isolate the material response, but reflect
the response of the complete experimental system. It is therefore necessary that due consideration be given
to such factors as specimen boundary conditions (end effects), stiffness of the testing apparatus and loading
rate [5,12].

Fig. 1. Equivalence of experiment and in situ mining stress states. (a) The extension test and the unconfined walls of excavations, (b)
the uniaxial compression test (UCS) and unconfined rock pillars and (c) the triaxial test and the highly confined compression zone
ahead of the deep level mining face, after Ortlepp [21].
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3037

The theory of the nucleation and growth of microcracks from pre-existing flaws is attributed to Griffith
[13] and is fundamental to the understanding of the quasi-brittle fracture process. Cracks initiating from
flaws grow and become stable, a process that is repeated for less critical flaws as the stress is raised, resulting
in stable pre-peak inelasticity [10,14–18]. Dilatancy is the most significant manifestation of quasi-brittle
fracture in compression and is defined as the increase in volumetric strain relative to the expected elastic
response [2]. Using acoustic emission Scholz [19] confirmed the observation of Brace et al. [10] that
dilatancy is directly proportional to microfracturing. Crack growth orthogonal to the direction of dilation
(parallel to the maximum compression stress) does not immediately produce a mechanical instability, as
observed in tensile fields. The direction of fracture propagation is necessarily orthogonal to the maximum
extensional strain and generally coincides with the direction of maximum applied compressive stress [15]. It
is this stable fracture process in compression that results in the large differences between tensile and
compressive strengths of quasi-brittle materials. This process ultimately leads to mechanical instability in
the post-peak region resulting in the formation of a macroscopic failure plane from the coalescence and
complex interaction of microcracks and the heterogeneous microstructure [1,2,15,20]. The post-peak drop
in strength that occurs is associated with the rupture of material linkages and the mobilisation of the
macroscopic failure plane.
This creates a global deformation mechanism that is generally associated with significant energy release and
may be the cause of so-called rockbursts, occurring when failure planes extend into mining excavations [21].
Many complex numerical models invoking the concepts of elastoplasticity, damage mechanics and sta-
tistical methods have been presented in the computational literature [22,24–26], but many are not widely
used due to the difficulty in determining material parameters and the high level of technical expertise re-
quired by the user. The aim of this paper is to present a pragmatic and relatively simple numerical model for
the degradation and subsequent discrete fracturing of quasi-brittle materials in compression. To this end the
widely used Mohr–Coulomb failure criterion is adapted and coupled with a fully anisotropic tensile smeared
crack model. The Mohr–Coulomb criterion is able to recover the salient features of the quasi-brittle response
within engineering accuracy, including dilation and has the advantage of being based on parameters that are
easily determined experimentally. Fracturing due to dilation is accommodated by introducing an explicit
coupling between the inelastic strain accrued by the Mohr–Coulomb yield surface and the anisotropic
degradation of the mutually orthogonal tensile yield surfaces. The introduction of discrete fractures into the
material is achieved by modification of the finite–discrete element code ELFEN [23,27], which allows the
finite element continuum to fracture into discrete blocks based on the selected failure criterion. The existing
fracture and contact logic was designed for tensile fracture with little post-fracture interaction of the frac-
tured rock, and had to be modified for the compressive stress regime where fractured subregions interact and
slide relative to each other after the initiation of fracturing. The proposed model represents a phenome-
nological approach in which micromechanical processes are only considered in terms of the average global
response. Isotropy of strength in compression is justified by assuming uniform material heterogeneity, while
accrual of inelastic strain and associated degradation of the tensile strength is necessarily anisotropic and
dependent on the loading direction. The effectiveness of the proposed numerical model will be assessed by
application to problems such as triaxial and plane strain rock tests, borehole breakouts and strip punch tests.

2. Modelling of discrete fracture

2.1. The discrete/finite element method

The spatially discretised form of the dynamic equilibrium equations at time N t is given by [23,24,28],
M€ _ N tÞ ¼ f ext ðN tÞ
f int ðuðN tÞÞ:
uðN tÞ þ Cuð ð1Þ
3038 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

The diagonal mass matrix, M is defined by,

X Z
nelem nnode
M¼ A e
qNTi ½N1 ; N2 ; . . . ; Nnnode  dV ð2Þ
e¼1 eV
i¼1

and C is the diagonal damping matrix ðC ¼ aMÞ defining the viscous damping term with constant factor a.
It is apparent that at discrete time N t equations (1) represent a coupled system of second order linear
differential equations with constant coefficients. The relevant class of initial boundary value problems is
recovered by invoking (2) and the necessary boundary conditions and initial values. Using a central dif-
ference time integration scheme and adopting a diagonal damping matrix C recovers the algebraic equation
for the velocity of the degree of freedom ui at time N þ1 t, giving,
N þ12 2Mii
Cii N Dt N
1 N þ12 fiext ðN tÞ
fiint ðuðN tÞÞ
u_ i ¼ 2u_ i þ 2 Dt : ð3Þ
2Mii þ Cii N þ1 Dt 2Mii þ Cii N þ1 Dt
N þ1
The displacement at time t is given by,
N þ1 1
ui ¼ N ui þ N þ2 u_ i N þ1 Dt: ð4Þ
The amount of damping selected depends on whether the problem is dynamic or quasi-static. For quasi-
static analyses, a dynamic relaxation approach is taken with a higher damping value of 0.3 leading to
convergence to the static solution after sufficient timesteps. It is noted that the use of difference equations
requires special consideration for initial conditions and prescribed quantities. The objective Green–Naghdi
stress update procedure [29], rotates the spatial state variables to the reference configuration in which
classical (material-non-linear-only) elastoplasticity is invoked to update the stress rate. It has been observed
that the Green–Naghdi stress rate is suitable for materials under small elastic strain conditions [24,30,31].
The kinematics of a general quasi-brittle discrete system is that of large displacements, large rotations and
small strain. The clear separation between kinematic non-linearity and material non-linearity, permits the
use of a small strain inelastic stress update in the reference configuration. The updates for the compression,
tension and extension stress regimes are described in the following sections.

2.2. Mohr–Coulomb model for compression

The Mohr–Coulomb criterion is used for the failure envelope in compression as a first order approxi-
mation equating the shear failure stress s to the sum of the internal friction and the inherent material
cohesion c0 , giving,
jsj ¼ c0
rn tan / ¼ c0
rn l; ð5Þ
where / is the internal friction angle and l is the coefficient of friction. It is noted that the normal stress rn
acting on an inclined plane is defined here to be negative in compression.
The Mohr–Coulomb failure criterion in principal stress space is given by,
1
ðr
2 max

rmin Þ þ 12ðrmax þ rmin Þ sin / ¼ c0 cos /; ð6Þ

where rmax and rmin are the maximum (most tensile) and minimum (most compressive) principal stresses
respectively, and are shown in Fig. 2. The linear approximation of the Mohr–Coulomb model is poorest at
low confinement, where the experimental strength-confinement curves exhibit greatest non-linearity (con-
vexity), but for many materials is a reasonable approximation for confinements below the brittle–ductile
transition. The Mohr–Coulomb failure criterion is independent of the intermediate principal stress r2
resulting in the failure plane manifesting itself necessarily parallel to this direction. The poor approximation
of the tensile region of the conventional Mohr–Coulomb yield surface to the mutually orthogonal tensile
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3039

Fig. 2. The compressive fracture model. (a) The isotropic Mohr–Coulomb yield surface with softening anisotropic tensile planes and
(b) Pi-plane representation indicating possible return-mappings.

planes is a major shortcoming in terms of application to quasi-brittle fracture and is addressed later in this
paper.
It is necessary to introduce so-called mobilised material parameters that realise hardening/softening with
respect to some hardening parameter such that the permissible elastic domain depends on the current state
of inelastic strain as well as the history of evolution. The effective plastic strain
ep is adopted here and at
time t is defined by,
Z t rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
p 2 p p

e ðtÞ ¼ e_ : e_ dt; ð7Þ
0 3
where e_ p is the inelastic strain rate tensor. The form of the hardening/softening curves for the mobilised
cohesion c0 ð
ep Þ and the mobilised friction angle /ð
ep Þ are established by considering laboratory test data. A
non-associative implementation of the Mohr–Coulomb elastoplasticity model is required for the recovery
of the correct physical dilation response in compression. Replacing the friction angle / in the Mohr–
Coulomb yield function (6) by the so-called dilation angle w, such that w < /, recovers the widely adopted
plastic flow potential W given by,
Wðr;
ep Þ ¼ 12ðrmax
rmin Þ þ 12ðrmax þ rmin Þ sin w
c0 cos w; ð8Þ
where c0 ¼ c0 ð
ep Þ, w ¼ wð
ep Þ. The relationship between the dilation angle w and observed dilatancy is
defined by the constant dilatancy observed near the peak stress. After the formation and subsequent
mobilisation of a macroscopic failure plane the dilation tends to zero.
Stress states falling outside the elastic domain are non-permissible and are returned to the yield surface
through the so-called return-mapping procedure that is associated with the accrual of inelastic strain. It is
the dependency of the yield surface on the inelastic strain that renders the return-mapping non-linear and
necessitates the adoption of an iterative Newton–Raphson solution scheme. The flow rule defines the
evolution of inelastic strain through the concept of the subgradient of a function and is referred to as
Koiter’s rule for the case of multisurface plasticity [32], giving,
X
m
e_ p ¼ k_ a or Wa ðr; jÞ; ð9Þ
a¼1

where Wa ðr; jÞ is the plastic flow potential function and ka is the corresponding plastic consistency
parameter (or plastic multiplier) for the yield function fa ðr; jÞ, with a 2 ½1; 2; . . . ; m. For non-associative
plasticity Wa ðr; jÞ 6¼ fa ðr; jÞ for some a 2 ½1; 2; . . . ; m. The Kuhn–Tucker complementary (loading/
3040 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

unloading) conditions are the optimality conditions for the closest point projection minimisation problem
and for all a 2 ½1; 2; . . . ; m are given by,

k_ a P 0; fa ðr; jÞ 6 0; k_ a fa ðr; jÞ ¼ 0: ð10Þ

The consistency (persistency) condition is given by,

k_ a f_a ðr; jÞ ¼ 0: ð11Þ

If madm 6 m is the number of constraints at a certain point ðr; kÞ 2 oEr , with associated indices belonging to
the set Jadm ¼ fb 2 ½1; 2; . . . ; mjfb ðr; jÞ ¼ 0g, then the actual active constraint set Jact is defined by,

Jact ¼ fa 2 Jadm jf_a ðr; kÞ ¼ 0g: ð12Þ

It is the determination of the actual active constraint set Jact that is one of the principal algorithmic concerns
in multisurface plasticity.

2.3. The tensile fracture model

The multiple fixed crack model [33] and the microplane model [34] are attempts to overcome the limi-
tations of the fixed crack [34] and rotating crack [37] concepts. The differences between these models are
discussed in detail elsewhere [24,36]. In fixed crack models, the elastic properties degrade or plastic strain
accumulates across a pre-defined plane inducing a coupling between the shear and normal stresses. The
rotating crack concept has no memory of the crack direction and damage accrues in the direction of the
current principal stresses [36].
The fixed crack model is observed to be overly constrained with induced shear stresses unable to invoke
effective reorientation of crack directions, as is the physical manifestation. Conversely, the rotating crack
model is under constrained, exhibiting much-reduced shear stresses (for coincident rotation of principal
stress and strain). The crack direction reflects the current state of damage and not the damage history. In
oscillating stress fields that realise large rotations of the principal directions, it is possible that the intro-
duced discrete fractures will not reflect the previous history of directions in which damage was accrued. The
multiple fixed crack model and the microplane model allow simultaneous cracking in multiple orientations
which relieves excessive shear constraints while introducing strong path dependence. However, these
models are overly complex for the simplifying assumptions on which they are based and generally exhibit
arbitrary coupling between damage directions and spurious energy dissipation.
The stabilised rotating crack model is now proposed that realises crack band orientations defined by the
weighted average of the extensional strain directions invoking damage. The crack band remains fixed
during unloading and reloading, with the average orientation reflecting the cracking history and thus
constraining spurious rotation. Reorientation of the crack band in response to a build up of shear stress
relaxes the fixed crack constraints during loading and more closely resembles the physical manifestation.
To determine the crack direction, consider a two-dimensional configuration at time N t consisting of a
single crack band with normal n at an average angle Nt hav to an appropriate reference frame (global or local
element) given by,
PNt t t
Nt t¼0 D xnn he1
hav ¼ P Nt ; ð13Þ
t
t¼0 D xnn

where D t xnn is the increment of scalar damage in direction n at time t and t he1 is the angle of the first
principal strain direction with respect to the reference frame at time t. Generally, the orientation of
maximum applied strain he1 does not coincide with the crack band normal orientation hav , but with damage
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3041

accrual in a given direction they are convergent. In 3-D the average orthotropic orientation is more
complicated, but can be derived in terms of direction cosines.
It is apparent that the stabilised rotating crack model violates the form invariance condition due to the
general non-coincident rotation of the principal axes of stress and strain. However, for systems which
experience large arbitrary rotations of crack bands due to dynamic effects, unloading or post-failure
interaction, the error introduced due to the omission of shear coupling terms is negligible compared to that
introduced through spurious crack band orientations. It is noted that fracture is an extensional process and
the modelling of the normal behaviour is of primary importance.
The manifestation of deformation during fracturing is accomplished with a crack band model [34,35,38].
This model exhibits a constant strain distribution and a linear variation of displacement across the width of
the band, whereas a cohesive crack model [39–41] exhibits a Dirac delta (spike) strain function and a
corresponding displacement jump across the crack. The crack band description is therefore more realistic at
the initiation of softening, where material damage is dispersed, while a cohesive crack is more realistic at the
end of softening, where dispersed damage coalesces to form a localised discrete crack. The deformation in a
crack band is equivalent to that of a cohesive crack if the fracturing elongation DLf ¼ ef hc is identified with
the cohesive crack opening displacement w, such that w ¼ ef hc . Identical stress-elongation curves for the
crack band model and the cohesive crack model require a correspondence between the constitutive func-
tions. It would appear that although the cohesive crack model and the crack band model are equivalent, the
crack band model has an extra parameter in the crack band width hc . Recognising that the crack band
width hc is used in defining the specific fracture energy, the crack band model is equivalently defined in
terms of a tensile strength ft and a specific fracture energy Gf that defines the amount of energy dissipated
per unit cross-sectional area of the crack band.
The crack band model can be formulated with many different strain-softening curves to describe the
energy dissipation. In general, exponential curves appear to yield better results when applied to concrete
systems [4,13,18,41], although this is not necessarily true for brittle fine-grained rocks and ceramics.
Nonetheless, the linear softening curve is able to recover the salient features of the quasi-brittle structural
response and is applied here as a useful first order approximation. For an arbitrary finite element mesh
strain-softening will necessarily localise into a single element of local length scale (average dimension) hðeÞ
c
and mesh objective dissipation of the specific fracture energy Gf is achieved here using a linear strain-
softening curve defined by,
 
r ¼ /ðef Þ ¼ ft 1
ef =efðeÞ
c ; ð14Þ
where efðeÞ
c is the element critical fracturing strain given by,
2Gf
efðeÞ
c ¼ ðeÞ
: ð15Þ
hc ft
For the crack band model with linear strain softening, the constant softening modulus H of the element has
a bi-linear behaviour given by,

dr hðeÞ f 2
H¼ p
¼
c t for efðeÞ 6 efðeÞ
c ;
de 2Gf ð16Þ
H ¼ 0 for efðeÞ > efðeÞ
c ;

where hðeÞ
c is the local element lengthscale (average dimension) and Gf is the specific fracture energy.
The selected implementation of the anisotropic rotating crack band model consists of three orthogonal
tensile planes
firi ðrÞ ¼ ri
fti
H hDepi i ¼ 0 ð17Þ
3042 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

that act independently in each of the principal stress ri directions and can thus realise anisotropic softening
in the material. Fig. 2(a) shows the coupling of the orthogonal, anisotropic tensile planes and the isotropic
Mohr–Coulomb yield surface, while Fig. 2(b) shows the Pi-plane representation of the composite yield
surface. The incremental update of the tensile strength in each of principal stress direction is given by,
nþ1
fti ¼ n fti þ H hDepi i for i ¼ 1; 2; 3; ð18Þ
where ÆÆæ are the McAuley brackets returning zero for negative (compressive) inelastic strain increments
Depi , and 0 fti ¼ ft . Tensile strengths associated with the principal stress directions at the end of the previous
timestep are associated with the current principal stress directions according to the relative proximity of
orientation. The assumptions on which the rotating crack band model is based restrict approximate validity
to only small rotations of the stress field. Any violation of tensile constraints, whether by the elastic trial
state or induced during the standard Mohr–Coulomb return-mapping, necessarily introduces the possibility
of a consistent return-mapping to any feasible combination of the tensile planes and the Mohr–Coulomb
yield surface. This is a consequence of the anisotropy of the tensile response, that no longer ensures the
validity of the constraint r1 P r2 P r3 for the updated principal stresses.

2.4. Explicit coupling between degradation in compression and tension

Starting from the assumption that quasi-brittle fracture is extensional in nature, any phenomenological
yield surface is effectively divided into regions in which extensional failure can be modelled directly, as in
the case of tensile stress fields and indirectly, as in the case of compressive stress fields. The violation of
either region of the yield surface should result in the realisation of inelastic extensional strain (cracks
opening). In both cases the definition of this inelastic strain is identical, it is the micro-mechanism that
produces the strain that differs. Fig. 3 depicts the dilation response of compressive failure in quasi-brittle
materials and clearly indicates the extensional inelastic strain directions associated with fracture. The vector
of principal increments of inelastic strain obtained from the single vector return-mapping to the Mohr–
Coulomb main-plane is given by [24],
8 p9 8 9
< De1 = < 1 þ sin w =
D
ep ¼ Dep2 ¼ Dka 0 ; ð19Þ
: p; : ;
De3 sin w
1
where the increment of inelastic strain Dep1 in the first principal stress direction is extensional. An explicit
coupling is therefore proposed between the extensional inelastic strain associated with the dilation response
of compressive failure and the tensile strength in the dilation direction, even in the absence of tensile stress.

Fig. 3. (a) Compressive loading with confining stress, (b) relationship between applied and volumetric strain and (c) compressive
failure with associated lateral extensional (positive) inelastic strain causing fracture and dilation.
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3043

Various methods have been applied for modelling of biaxial stress states e.g. [25,26,36]. The proposed
implementation of an explicit coupling between compressive stress induced extensional strain and tensile
strength degradation permits the realisation of discrete fracturing in purely compressive stress fields. Such
increments of extensional strain must be associated with tensile strength degradation in the parallel
direction, giving,
fti ¼ fti ðepi Þ where epinþ1 ¼ epin þ Depi ; ð20Þ
where fti is the tensile strength in the ith principal stress direction. Thus, for the example shown in Fig. 3c,
i ¼ 1. In pure, triaxial, tension, the strength degradation will occur in all three principal stress directions.
For a biaxial state of stress with one direction in tension, and the orthogonal direction in compression, the
degradation will occur perpendicular to the tensile stress, with the possibility of additional damage
occurring in the same direction if the compressive stress exceeds the strength.
The update of the tensile strength is only performed once consistency has been achieved for a return-
mapping including a compressive yield surface. It is noted that the tensile strength in a particular direction
is not updated if the associated tensile yield plane is included in the consistent return-mapping. The
principal stress return-mapping is a specialisation of the standard elastic predictor/plastic corrector schemes
widely employed in computational plasticity [29], with the return-mapping procedure (associated with the
plastic corrector phase) carried out in the principal stress space. This differs markedly from the Mohr–
Coulomb implementations generally encountered in the literature that adopt the invariant representation of
the model. The principal stress approach adopted here results in a far simpler and more efficient compu-
tational implementation. It appears that most non-linearities of the return-mapping are confined to the
spectral decomposition, whereby the principal stresses are evaluated as non-linear functions of the stress
tensor. The feasible return-mappings are numbered in Fig. 2(b) and detailed by Klerck [24]. The two
algorithms proposed by Simo et al. [42] for calculating the return to the active yield surface assume that the
final successful return-mapping includes the initial set of trial constraints. This is not always the case for the
current model and so an algorithm proposed by Peric and de Souza Neto [44] is adopted that assumes
nothing about the nature of the successful return-mapping and exhibits the ability to attempt all permissible
constraint combinations. Generally it is impossible to predetermine which return-mapping will be suc-
cessful and it is therefore necessary to design an algorithm that is able to attempt all return-mappings in
turn until consistency is achieved. In an attempt to maximise computational efficiency it is possible to order
the return-mappings with respect to feasibility and eliminate return-mappings that are known not to be
feasible. A discrete crack is introduced when the tensile strength in a principal stress direction reaches zero
and is orientated orthogonal to this direction.

2.5. The topological update

A discrete crack is introduced when the tensile strength in a principal stress direction reaches zero and is
orientated orthogonal to this direction, as already described with respect to the crack band model. The
insertion of discrete cracks into the quasi-brittle continuum follows three steps.

Step 1. Create a non-local failure map (weighted nodal averages).

Step 2. Determine fracture feasibility and the order of discrete crack insertion.
Step 3. Perform the topological update (remeshing).

The most meaningful quasi-brittle damage indicator or so-called failure factor is the ratio of the inelastic
fracturing strain ef to the critical fracturing strain efc . The local fail factor Fk at Gauss point k is given by,
Fk ¼ ðef =efc Þk ; ð21Þ
3044 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

where Fk is associated with a local fracture direction hk , which is normal to the local failure (softening)
direction. Discrete fracture at the Gauss point level is effectively realised on dissipation of the complete
fracture energy (zero strength) and corresponds to a failure factor greater than unity. However, the nodal
basis of the finite element discretisation renders discrete crack insertion simpler if associated with nodal
failure factors. It is with this motivation that a non-local failure map is sought, based on the weighted
nodal averages of immediately adjacent element (Gauss point) failure factors and fracture directions. The
weighted-average failure factor F p and fracture direction
hp at node p are given by,
Nadj
,N Nadj
,N
X X adj X Xadj

Fp ¼ Fk wk
k ¼
wk ; h
k wk
h wk ; ð22Þ
k¼1 k¼1 k¼1 k¼1

where Nadj is the number of immediately adjacent Gauss points and wk is a weighting factor usually taken as
the element volume.
A discrete fracture will necessarily be realised through a nodal point if the nodal failure factor is greater
than unity. The possibility of adjacent nodes reaching failure factors greater than unity within the same
time increment renders it necessary to determine a discrete fracture sequence prior to updating the
topology. It is apparent that the sequence can affect the final fracture pattern. The physical manifestation of
multiple crack propagation is firmly based on energetic feasibility, with the most energetically feasible
fracture propagating first. The magnitudes of the failure factors represent the numerical analogy of ener-
getic feasibility and their relative magnitudes determine the sequence of topological update.
The discrete crack insertion and the associated topological update commences with the definition of a
fracturing plane through the failed nodal point and in the direction of the weighted average fracture
direction, as shown in Fig. 4(a). The plane represents the actual discrete crack orientation and does not
generally coincide with the existing element sides. Insertion of the actual crack direction therefore requires
local remeshing and the possible creation of new elements, as shown in Fig. 4(b). This intra-element
fracturing may result in the creation of unfavourably dimensioned sliver elements if the underlying dis-
cretisation is simply split and may decrease the time increment stability limit. In this case local adaptive
mesh refinement must be undertaken to achieve an acceptable element topology. Alternatively, an element
threshold dimension can be defined below which the discrete crack is snapped to the most favourably
orientated existing element side, as shown in Fig. 4(c). Inter-element fracturing in which no remeshing is
considered is not usually desirable, requiring very fine meshes to capture the fracture directions accurately.
Following discrete crack insertion the failed elements are reinitialised under the assumption that the
element damage is coalesced into the discrete crack formation. Ideally only the damage normal to the dis-
crete crack should be zeroed, but with remeshing and the creation of additional elements the complexity of
state variable mapping is avoided by assuming that all state variables in the failed elements are reinitialised.
The contact response of the newly created surfaces is of vital importance to the stability of the solutions
after fracture initiation. The conventional approach for discrete contact in finite element methods is to

Fp , θ p

(a) (b) (c)

Fig. 4. (a) Weighted-average nodal failure direction, (b) intra-element fracturing and (c) inter-element fracturing (the dashed lines
denote the potential crack planes and the arrows indicate the fracturing strain direction).
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3045

apply the penalty method to approximately satisfy the contact displacement constraints for finite penalty
coefficients [27,43]. This results in the finite penetration of contact surfaces and an effective reduction in the
volume of contacting entities. This causes an instability when the region is in compression resulting in
spurious energy release. Thus, the correct contact description requires the exact satisfaction of the dis-
placement constraint. The explicit Lagrangian method was proposed for the determination of the contact
forces that satisfy the displacement constraints exactly [45,46]. The classical Lagrange multiplier method
evaluates the contact forces exactly, but is not compatible with the explicit integration operator. Explicit
compatibility is established by referencing Lagrange multipliers (contact forces) one time increment ahead
of associated surface contact displacement constraints. The explicit Lagrange multiplier method is not
purely explicit because a coupled system of equations must be solved to obtain the Lagrange multipliers.
The equation solving strategy adopted here is the symmetric Jacobi method, in which non-linear contact
force conditions (normal tension limit and tangential friction) are enforced during each iteration [24].
Tangential displacement constraints (friction) are enforced using the penalty method. For a stable solution,
the explicit Lagrange multiplier method requires a reduction in the critical time step to one-tenth of the
critical timestep for the continuum analysis. The classical non-associative Coulomb friction law defining a
critical tangential force proportional to the normal contact stress, fn and frictional coefficient l is applied
once the fracture has been placed in the finite element mesh.
Discrete fracture is realised after localisation of damage into crack bands and is therefore prone to the
discretisation dependence of the smeared crack models. This affects the positioning and the spacing of
fractures, but not the predicted orientation. The orientation is determined by the non-local averaging of
local element failure directions and is effectively independent of the discretisation. As already mentioned,
discretisation dependence can be minimised by utilising high-density, non-biased meshes of low order
elements. The context in which discrete fracturing has been introduced so far has been that of localised
failure in individual crack bands. This type of failure generally originates from structural features that
concentrate stress. Generally, objective fracture energy dissipation for distributed failure requires the use of
regularisation methods that ensure objective dissipation over finite volumes through the introduction of a
spatial lengthscale, which have been considered elsewhere [24].

3. Numerical examples

3.1. The strip punch test

Consider the strip punch test investigated experimentally by Dede [47] for quartzite rock to simulate the
failure of stabilising pillars in deep level gold mines. Prismatic quartzite specimens of dimensions 80 · 80 · 80
mm were machined on one side to produce a bisecting 1 mm high and 10 mm wide strip, as shown in Fig. 5(a).
The normal loading of the raised strip was performed within a triaxial loading cell, with 4.5 MPa of confining
pressure applied to the specimen sides. To reduce friction all loading surfaces were coated with stearic acid
and so the platen friction coefficient is set to 0.005. The material friction coefficient is selected to be 1.2.
Although the stress state is three-dimensional, the geometry can be well described in plane strain. The
adopted plane strain finite element mesh is shown in Fig. 5(b) and consists of a graded distribution of three-
noded, linear elements with single point Gaussian quadrature. The mesh is observed to be refined in the
region of the punch strip, with elements of dimension 1 mm. The 4.5 MPa confining pressure is applied to
the specimen via discrete steel platens (shown in dark) and the normal strip load or axial load is applied via
a prescribed velocity of 0.1 m s
1 . Material properties for quartzite are defined by Dede [47]. The input
material parameters and the evolution of the cohesion, friction and dilation were found by matching the
results of a single element triaxial test simulation with the data [47] and are given in Table 1. The Young’s
modulus is 68 GPa, the Poisson’s ratio is 0.17, the density is 2800 kg m
3 , the tensile strength is 27 MPa.
3046 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

Fig. 5. (a) Strip punch geometry and loading configuration, after Dede [47] and (b) plane strain finite element model.

As mentioned in Section 2.5, the introduction of discrete fracturing requires a reduction in the critical
time step to a tenth of the continuum value. The timestep for the discrete fracture model is 0.49E)8 s and
the model takes 8.5 h to run on a 1.8 GHz computer, compared to 57 min for the continuum models. Thus,
the effect of varying the fracture energy was studied using continuum analyses to save time and to provide a
comparison to evaluate the effect of introducing the discrete fractures. Fig. 6(a) shows a typical plot of the
experimental normal punch stress versus the specimen axial strain and the plots recovered from
the numerical system. The initiation of purely compression induced microfracturing occurs at the edges of
the punch strip, where high compressive stress concentrations exists. The first significant feature of the axial
stress–strain response defined in Fig. 6(a) is the drop in stress experienced between points A and B. This
stress drop is associated with the localisation and coalescence of material degradation about the periphery
of the conical microfracture zone to form the macroscopic failure zone shown in Fig. 7(a). The hardening of
the axial stress–strain response between points B and C is attributed to the continued formation and
mobilisation of the conical failure zone and the initiation of the downward wedging action. The purely
tensile nature of the stress field propagating the cleavage crack is verified by observing the absence of
significant compressive inelastic strain vectors beneath the conical microfracturing zone, as shown in Fig.

Table 1
Evolution of material parameters with plastic strain for three rock types
Quartzite
Plastic strain 0 0.0005 0.001 0.0015 0.002 0.004 0.008 0.03
Cohesion (MPa) 53 50 45 42 37 28 10 10
Friction angle () 5 28 34.8 38.8 43.2 46 46 46
Dilation angle () 0 28.4 28.4 28.4 28.4 28.8 20 20

Sandstone
Plastic strain 0 0.0005 0.001 0.0015 0.002 0.008 0.01 0.03
Cohesion (MPa) 20 19 18 17 16 4 0 0
Friction angle () 17 31.5 40.5 52.5 52.5 52.5 52.5 52.5
Dilation angle () 0 5.75 11.5 17.25 23 23 18 4

Tyndall limestone
Plastic strain 0 0.0005 0.002 0.0055 0.011 0.02 0.022 0.03
Cohesion (MPa) 14 14 13 11 8 4 4 4
Friction angle () 4.4 22.4 27.2 29.2 29.6 30.4 30.4 30.4
Dilation angle () 0 19.6 19.6 19.6 19.6 19.6 5 5
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3047

600 1400
stress (compressive models) (MPa)

stress (compressive models) (MPa)

500 C 1200
A tensile only model
1000
400
B Gf=3500
800
300
Gf=1500 600
D Gf=150
experimental data A C
200
400
fracture model Gf=150 B experimental data
100 200 D

0 0
0.000 0.002 0.004 0.006 0.008 0.000 0.002 0.004 0.006 0.008
(a) strain (b) strain

Fig. 6. (a) Numerical axial stress vs. axial strain plot for compressive fracture models with different parameters. Experimental data
after Dede [47]. (b) Numerical axial stress vs. axial strain plot for tension only fracture model.

Fig. 7. (a) Initiation of material degradation at the edges of the punch strip, (b) plastic strain vectors describing the conical micro-
fracture zone beneath the punch strip, (c) full mobilisation of the conical microfracture zone beneath the punch strip and completed
formation of the axial cleavage crack, (d) the observed fracture pattern and (e) discrete fracture pattern for a tension only model.

7(b). The formation of the axial cleavage crack and the full mobilisation of the conical zone beneath the
punch strip produces a global deformation mechanism resulting in the complete failure of the specimen
between points C and D with a fracture pattern as shown in Fig. 7(c). The numerical discrete fracture
pattern compares very well with the typical experimental fracture pattern presented by Dede [47] and shown
in Fig. 7(d). The en-echelon fractures about the periphery of the conical zone are clearly observed, as is the
3048 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

axial cleavage crack. Some coalescence of the en-echelon fractures is observed in the numerical fracture
pattern, although it is noted that the spacing of cracks is necessarily set by the mesh and the recovery of a
complete failure surface on which sliding can be realised would require considerable mesh refinement. The
conical zone beneath the punch strip in the numerical model is observed to be longer than that observed in
experimental results. This difference is attributed to the effect of mesh alignment on the local numerical
model. The data presented by Dede [47] shows the overall deformation of the test machine and the loading
platens. The deformations were multiplied by a calibration factor of 0.24 to produce the response across the
sample only.
The graph in Fig. 6(a) also shows three stress–strain curves from continuum models, with different values
of fracture energy Gf , where the discrete fracturing has not been activated. As expected, the lowest value of
Gf , corresponding to the expected value for the quartzite, produces the most brittle response, with no
regaining of the load after failure. Higher values of Gf results in deformation at almost constant stress levels
prior to the final stress drop, which happens gradually. The model with fracturing activated also shows
deformation at constant stress, but a much more sudden stress drop at the final failure. The stress levels and
deformations agree well with the experimental data.
It is observed that the compressive fracture model, exhibiting explicit compressive–tensile coupling and
the softening of the orthogonal tensile planes, produces the most brittle response, as expected.
In order to investigate the effect of including the compressive fracture model on the load deformation
response and the fracture pattern, an analysis was performed with only tensile fracturing active i.e. by
setting the compressive strength to a very high value. Fig. 7(e) shows the final fracture pattern when a
purely tensile rotating crack model is adopted as the material model. It is apparent that only the axial
cleavage crack is correctly recovered, with a complete absence of material degradation in the vicinity of the
punch strip. This explains the excessive peak stress recovered (1300 MPa) and it is concluded that accurate
numerical modelling of all punch or indenter processes requires the accurate modelling of material deg-
radation obtained with the compressive fracture model in both tensile and compressive stress fields.

3.2. The plane strain and triaxial compression tests

Yumlu and Ozbay [8] tested 30 · 30 · 10 mm prismatic specimens of coal, sandstone, norite and quartzite
in plane strain and triaxial conditions, of which only the experiments on sandstone will be considered here.
Fig. 8(a) shows the plane strain finite element model. For the plane strain experiments, the confining
pressure was applied via a stiff system of discretely defined steel platens. The use of a threaded bolt across
the sample (shown in the schematic diagram of Fig. 8(b) to provide confinement in the plane strain test)
provides a constant stiffness and not a constant stress boundary condition. In the model (Fig. 8(a)), the
confinement is applied with an initial applied displacement to the outside of the loading platen and the
stiffness of the platen determines the subsequent system stiffness. Fig. 8(c) depicts the quarter symmetry
triaxial finite element model. The hydraulic confinement of the triaxial test is modelled using a soft, applied
stress condition as shown in Fig. 8(c). In both models, the ends of the specimen are loaded axially via
prescribed velocities of 0.1 m s
1 . A uniform mesh of linear plane strain elements of average dimension 0.5
mm and single point Gaussian quadrature is adopted. Material properties for the sandstone are defined by
Yumlu and Ozbay [8] and the evolution of the cohesion, friction and dilation is given in Table 1. The
Young’s modulus is 28 GPa, the Poisson’s ratio is 0.25, the density is 2800 kg m
3 , the tensile strength is 4
MPa and Gf ¼ 100 J m
2 . The explicit Lagrangian contact algorithm is adopted for discrete contact, with
the platen friction coefficient being 0.01 and the material friction coefficient is 1.3.
The best fit of the Mohr–Coulomb response to the triaxial test is used to determine the evolution of the
cohesion, friction and dilation during the test. Fig. 9 shows the axial stress–strain curves obtained from the
experiments on sandstone and the models. The experiments exhibit an initial non-linearity due to com-
paction of irregularities on the sample platen interface that has not been included in the models, and re-
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3049

Fig. 8. (a) Plane strain finite element model, (b) schematic of plane strain experiment and (c) triaxial finite element model.

200

180

160
plane strain model
140
stress (MPa)

120

100

80
plane strain data
60

40 triaxial model

20 triaxial data

0
0.00 0.01 0.02 0.03
strain

Fig. 9. Observed data (after Yumlu and Ozbay [8] with initial settling-in strains removed), and modelled compressive axial stress–
strain curves for sandstone specimens tested under plane strain and triaxial conditions subject at a confining stresses of 3 MPa.

moved in the comparisons shown in Fig. 10. Considering the rest of the stress–strain response, the triaxial
model shows the pre-failure plastic strain and low residual stress also observed in the experiments. The
plane strain model shows a linear response with post-failure hardening. As shown in Fig. 10, the models are
able to reproduce the correct strength increase with confinement and the strength difference between the
triaxial and the plane strain models. Yumlu and Ozbay [8] observed that the plane strain constraint in-
creases the peak strength and reduces pre-peak-non-linearity, and attribute this to different material re-
sponses for the different stress paths. However, as both models used the same material parameters, the
3050 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

250

y= 10.265x + 121.19
200 R2 = 0.9754

axial stress (MPa) 150 y= 9x + 98.5

R2 = 0.9996

100
plane strain model
triaxial model
triaxial data
50 plane strain data
Linear (planestrain model)
Linear (triaxial model)
0
0 2 4 6 8 10
confinement (MPa)

Fig. 10. The triaxial and plane strain principal stresses at failure, for the data [8] and models.

modelling indicates that the increased strength and residual hardening are not due to material properties,
but are a response to the stiffer loading system in the plane strain test.
The evolution of the discrete fracture pattern post-localisation of material failure in plane strain is shown
in Fig. 11a. The final pattern compares well with the observed failure shown in Fig. 12b. The constraint of
inter-element fracturing does not permit the resolution of exactly vertical en-echelon fractures. The lo-
calisation occurs within one element width as expected for the crack band model, but would need to be
refined to exactly match the observed fracture zone width. The localisation of failure in the triaxial models
is similar to those of the plane strain example, as also noted for the experiments [8].

3.3. The borehole breakout

The fracture distribution around circular cavities has been widely observed in the literature and well
documented by Hoek and Brown [48], Carter et al. [49,50], Sellers and Klerck [51], Lac du Bonnet granite

Fig. 11. (a) The evolution of discrete fracture in the localised failure zone, obtained using the compressive fracture model in plane
strain and (b) the failure observed in the experiment.
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3051

Strain gauges
Remote fractures

Primary fractures

Sidewall fractures

(a) (b)

Fig. 12. (a) Primary, remote and sidewall fracture distribution around a circular cavity in uniaxial compression. (b) The location of
strain gauges to indicate the initiation of different fracture types, after Carter [56].

tested by Lee and Haimson [52] or Berea sandstone tested by Ewy and Cook [53] and Zheng et al. [54]
amongst many others. Fig. 12(a) shows the location of the so-called primary, remote and sidewall (borehole
breakout) fractures around a circular cavity loaded in uniaxial compression [50]. The primary fractures
initiate first from the tensile stress concentrations at the crown and invert of the circular cavity and
propagate away from the cavity parallel to the direction of applied loading. The remote fractures form next,
as shown in Fig. 12(a), initiating approximately adjacent to the termination of the primary fractures and at
some distance from the circular cavity. These remote fractures are associated with the redistribution of
tensile stress by the primary fractures and propagate away from and towards the cavity sidewalls. The four
possible remote fractures do not generally initiate simultaneously and occasionally not all four fractures
appear [50]. The third type of fracture is the sidewall fracture or borehole breakout that initiates with the
compressive failure of the cavity sidewalls. The global failure mechanism observed in experimental speci-
mens is generally due to the link up of the remote and sidewall fractures.
The mechanisms of borehole breakout in weak sedimentary rock, such as Cardova cream limestone
considered by Haimson and Song [55] as well as the stronger and more brittle Lac du Bonnet granite tested
by Lee and Haimson [51] were successfully modelled by Klerck [24]. The change in the fracture pattern
surrounding borehole breakouts in Elsburg quartzite due to the presence of interfaces was also successfully
simulated [51]. In this paper, the uniaxial compression tests performed by Carter and coworkers [50,56] on
specimens of Tyndall limestone (Tyndallstone) are modelled. The aim of the numerical investigation is to
demonstrate that the proposed compressive fracture model is able to recover the correct evolution of the
primary, remote and sidewall fractures and exhibit the observed size effect. To the authors’ knowledge this
has not yet been achieved in the literature and would constitute a significant contribution to the study.
Circular hole diameters ranging between 3.2 and 62 mm, were used to investigate the effect of hole size on
fracture initiation around circular underground openings. All the blocks were 89 mm thick. Specimens were
cut and polished into prismatic rectangular shapes with central circular holes and loaded via steel platens
with teflon inserts. Strain gauges were positioned, as shown in Fig. 12(b), to indicate the initiation of
primary, remote and sidewall fractures.
Numerical plane strain models of specimens tested by Carter et al. [50], with radii (R) of 6.4, 16.3 and 25
mm are illustrated in Fig. 13. Material properties for the limestone are deduced by matching the results of a
single element triaxial test simulation with the triaxial data presented by Carter [56] and Carter et al. [50],
and the evolution of the cohesion, friction and dilation is given in Table 1. The Young’s modulus is 21 GPa,
3052 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

Fig. 13. The finite element models of the uniaxial loading of Tyndall limestone rectangular prisms with centrally located circular holes
of radii (a) 25 mm (black dots mark positions of numerical ‘‘strain gauges’’), (b) 6.4 mm and (c) 16.3 mm.

the Poisson’s ratio is 0.3, the density is 2000 kg m
3 , the tensile strength is 1.5 MPa and Gf ¼ 9 J m
2 . The
platen friction coefficient being 0.05 and the material friction coefficient is 0.58. Symmetry of the experi-
mental system is invoked for convenience in the finite element models proposed in Fig. 13. Uniaxial stress
loading is applied via discrete steel platens.
The results obtained using the compressive fracture model are shown in Fig. 14 for the specimen with
radius R ¼ 25 mm. Fig. 14 shows the evolution of the discrete numerical fracture pattern for the distinct
stages of primary, remote and sidewall fracture. Clearly demonstrated is the lengthening of the remote
fractures towards and away from the borehole walls, leading ultimately to the complete failure of the

Fig. 14. Evolution of discrete fractures in borehole breakout model immediately after the stages of (a) primary initiation, (b) remote
initiation and (c) sidewall initiation as indicated in Fig. 15.
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3053

specimen. The correspondence of the numerical response to the physical response observed by Carter and
coworkers [49,50,56] is excellent.
The size effect of the initiation stresses of primary, remote and sidewall fractures around circular cavities
has been widely observed in the literature, including Carter et al. [50], Lee and Haimson [52], Lajtai [57],
Martin et al. [58], Martin [59]. The general trend is a decrease in fracture initiation stress with increasing
cavity size and is of great importance in that physical experimental modelling is usually aimed at appli-
cations of a significantly larger scale. No single theory has yet been able to account for the observed size
effect, although limited success has been achieved by Ewy and Cook [60] using fracture mechanics theory,
Gonano [61] using a critical strain energy model, Lajtai [57] using stress averaging techniques and Santarelli
and Brown [62] using pressure dependence of the elastic modulus. In order to assess whether the com-
pressive fracture model is able to recover the physically observed size effect it is necessary to determine the
stresses at which the different fractures initiate. This is achieved by monitoring the principal stresses in
elements in the region of the primary, remote and sidewall fractures. The relevant element positions are
marked in Fig. 13a. Initiation of the primary and remote fractures is indicated by a change in curvature of
the axial stress versus the most extensional principal strain curve, while sidewall fracture is indicated by a
change in the curvature of the axial stress versus the most compressive principal strain curve, as shown in
Fig. 15. This method of indicating the initiation of fracture approximates the physical strain gauges used by
Carter and coworkers [50,56]. The stresses agree well with the relevant experimental values. The observed
strains are larger than the modelled strains. However, the monitoring of the stress and strain at individual
element Gauss points is a local indicator, whereas the physical strain gauges are effectively non-local
indicators. Also, the Tyndall Limestone is reported to have a significantly lower modulus in tension [56]
which could account for the higher observed strains, as a single modulus is used for compression and
tension in this analysis. In this experiment, the non-uniform stress state results in different parts of the
model experiencing post-peak softening although the overall load–deformation response is hardening. The
softened parts occur when cracks form and cause additional deformation. The good correspondence of
the modelled and observed crack initiation stress levels and the resulting change in shape of the overall
stress–strain curves at three separate points in the model as shown in Fig. 16, indicates that the model is
correctly predicting the post-peak stress redistribution in the sample.
The experimental and numerical initiation stresses for primary, remote and sidewall fractures are plotted
with respect to the hole radii in Fig. 16. Fracture initiation stresses are observed to decrease with increasing

24
sidewall initiation
20 remote
Axial stress (MPa)

primary
16
remote
sidewall remote
sidewall
12 initiation
stress
primary stress
8
stress
4
primary initiation

0
-6000 -4000 -2000 0 2000 4000 6000
Microstrain

Fig. 15. Experimental (light lines) and numerical (dark lines) stress–strain response indicating initiation of primary (point p in Fig. 13),
remote (point r in Fig. 13) and sidewall (point s in Fig. 13) fractures for R ¼ 25 mm model.
3054 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

20 40

15
Axial stress (MPa)

30

Axial stress (MPa)

experimental experimental
10 20
numerical numerical

5 10

0 0
0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35
(a) Hole radius (mm) (b) Hole radius (mm)

40

30
Axial stress (MPa)

experimental
20
numerical

10

0
0 5 10 15 20 25 30 35
(c) Hole radius (mm)

Fig. 16. The variation of fracture initiation stress with hole radius. (a) Primary, (b) remote and (c) sidewall for modelled results and
data after Carter [56].

hole size and appear to approach different asymptotes for large cavity sizes. The numerical results obtained
using the compressive fracture model appear to closely approximate the observed size effect. This is cer-
tainly understandable with respect to the primary and remote fractures that initiate in purely tensile stress
fields, as the crack band model is effectively invoked with a response normalised by the material fracture
energy. However, the compressive failure model is purely local and yet still appears to recover the physically
observed size effect. This may be due to the fact that failure is not being localised into a single band of
elements, but in finite regions determined by the geometry. Another more likely explanation is that the
correct size effect response of the primary and remote fracture initiation offsets the compressive mode
sidewall fracture initiation.

4. Conclusions

A method for modelling discrete fracture in geomaterials under tensile and compressive stress fields has
been developed. The method is based on modifications to an explicit discrete element/finite element code.
The model considers a failure envelope consisting of the Mohr–Coulomb failure surface in compression and
three independent anisotropic rotating crack models in tension. To model fracture in the extension
mechanism observed in compressive stress fields, the plastic strain induced by compressive failure is coupled
to anisotropic rotating crack models. Once the tensile strength has been depleted, discrete fractures are
inserted into the finite element mesh. An explicit Lagrangian contact algorithm is used to enforce non-
penetration of the newly created fracture surfaces.
The model is applied to triaxial and plane strain tests as well as punch tests and borehole breakouts.
These models exhibit fracturing in a range of stress states, sometimes with different modes in a single
sample. The model is able to quantitatively predict the stresses at the initiation of fracturing, the defor-
 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056 3055

mations associated with fracturing and the appropriate load–displacement response of the structure. The
evolution of discrete fracturing is predicted and the fracture patterns are very similar to those observed in
the experiments. Future work will concentrate on three-dimensional models to supplement the two-
dimensional applications presented here.

References

[1] D.K. Hallbauer, H. Wagner, N.G.W. Cook, Some observations concerning the microscopic and mechanical behaviour of
quartzite specimens in stiff, triaxial compression tests, Int. J. Mech. Min. Sci. Geomech. Abstr. 10 (1973) 713–726.
[2] P. Tapponnier, W.F. Brace, Development of stress-induced microcracks in westerly granite, Int. J. Mech. Min. Sci. Geomech.
Abstr. 13 (1976) 103–112.
[3] J.A. Franklin, Triaxial strength of rock materials, Rock Mech. 13 (3) (1971) 86–98.
[4] W. Janach, Failure of granite under compression, Int. J. Rock Mech. Sci. Geomech. 14 (1977) 209–215.
[5] K. Mogi, Some precise measurements of fracture strength of rocks under uniform compressive stress, Rock Mech. Engrg. Geol.
4 (4) (1966) 41–55.
[6] K. Mogi, Effect of the intermediate principal stress on rock failure, J. Geophys. Res. 72 (20) (1967) 5117–5131.
[7] K. Mogi, Fracture and flow of rocks, in: A.R. Ritsema (Ed.), The Upper Mantle, Tectonophysics, vol. 13 (1–4), 1971, pp. 541–568.
[8] M. Yumlu, M.U. Ozbay, A study of the behaviour of brittle rocks under plane strain and triaxial loading conditions, Int. J. Rock
Mech. Min. Sci. Geomech. Abstr. 32 (7) (1995) 725–733.
[9] E. Hoek, Z.T. Bieniawski, Brittle fracture propagation in rock under compression, Int. J. Fract. 1 (1965) 137–155.
[10] W.F. Brace, B.W. Paulding, C. Scholz, Dilatancy in the fracture of crystalline rocks, J.Geophys. Res. 71 (16) (1966) 3939–3953.
[11] E.Z. Wang, N.G. Shrive, Brittle fracture in compression: mechanisms, models and criteria, Engrg. Fract. Mech. 52 (6) (1995)
1107–1126.
[12] J. Ryder, CSIR Rock Engineering Textbook, C.S.I.R. (Miningtek), Johannesburg, South Africa, 1999.
[13] A.A. Griffith, Theory of rupture, 1st International Congress Appl. Mech., Delft, 1924, p. 55.
[14] W.F. Brace, E.G. Bombolakis, A note on brittle crack growth in compression, J. Geophys. Res. 68 (12) (1963) 3709.
[15] A.V. Dyskin, L.N. Germanovich, H. Joer, J.S. Krasinski, K.K. Lee, J.-C. Roegiers, E. Sahouryeh, K.B. Ustinov, Some results on
three-dimensional crack propagation in compression, in: H.P. Rossmanith (Ed.), Mechanics of Jointed and Faulted Rock,
Balkema, Rotterdam, 1995, pp. 91–96. ISBN 90-5410-5410.
[16] E.G. Bombolakis, Photoelastic investigation of brittle crack growth within a field of uniaxial compression, Tectonophysics 1 (4)
(1964) 343–351.
[17] H. Horii, S. Nemat-Nasser, Compression induced microcrack growth in brittle solids: axial splitting and shear failure, J. Geophys.
Res. 90 (B4) (1985) 3105–3125.
[18] E. Sahouryeh, A.V. Dyskin, The mechanism of 3-D crack growth in uniaxial compression: the role of non-singular stresses, Int. J.
Rock Mech. Min. Sci. 34 (3–4), 1997, Paper No. 071.
[19] C.H. Scholz, Microfracturing and the inelastic deformation of rock in compression, J. Geophys. Res. 73 (4) (1968) 1417–1432.
[20] T. Engelder, Joints and shear fractures in rock, in: B.K. Atkinson (Ed.), Fracture Mechanics of Rock, Academic Press Inc, 1987.
[21] W.D. Ortlepp, Rock fracture and rockbursts an illustrative study, The South African Institute of Mining and Metallurgy,
Johannesburg, SA, 1997. ISBN 1-874832-67-6.
[22] E. Sellers, F. Scheele, Prediction of anisotropic damage in experiments simulating mining excavations in Witwatersrand Quartzite
blocks, Int. J. Rock Mech. Min. Sci. 33 (7) (1996) 659–670.
[23] A. Munjiza, D.R.J. Owen, N. Bicanic, A combined finite-discrete element method in transient dynamics of fracturing solids,
Engrg. Comput. 12 (1995) 145–174.
[24] P. Klerck, The finite element modelling of discrete fracture in quasi-brittle materials. Ph.D. Thesis. University of Wales, Swansea,
2000.
[25] R. De Borst, Numerical aspects of cohesive-zone models, Eng. Frac. Mech. V 70 (2003) 1743–1757.
[26] P.H. Feenstra, Computational aspects of biaxial stress in plain and reinforced concrete, Ph.D., Delft University, 1993.
[27] J. Yu, A contact interaction framework for numerical simulation of multi-body problems and aspects of damage and fracture for
brittle materials, Ph.D. dissertation, University of Wales, Swansea, 1999.
[28] T. Belytschko, An overview of semidiscretisation and time integration procedures, in: T. Belytschko, T.J.R. Hughes (Eds.),
Computational methods for transient analysis, vol. 1, Elsevier Sci. Publ., 1983, pp. 1–65.
[29] E.A. De Souza Neto, D. Peric, D.R.J. Owen, Continuum modelling and numerical simulation of material damage at finite strain,
Arch. Comput. Meth. Engrg. 5 (November) (1998) 311–384.
[30] L. Szabo, M. Balla, Comparison of some stress rates, Int. J. Solids Struct. 25 (3) (1989) 279–297.
[31] J. Bonet, R.D. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge University Press, 1997.
3056 P.A. Klerck et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 3035–3056

[32] W.T. Koiter, Stress–strain relations, uniqueness and variational theorems for elasticplastic materials with a singular yield surface,
Quart. Appl. Math. 11 (3) (1953) 350–354.
[33] R. De Borst, P. Nauta, Non-orthogonal cracks in a smeared finite element model, Engrg. Comput. 2 (1985) 35–46.
[34] Z.P. Bazant, J. Planas, Fracture and Size Effect in Concrete and Other Quasibrittle Materials, CRC Press LLC, 1998.
[35] Z.P. Bazant, B.-H. Oh, Crack band theory for fracture of concrete, Mater. Struct. 16 (1983) 155–177.
[36] K.J. Willam, E. Pramonno, S. Sture, Fundamental issues of smeared crack models, in: S.P. Shah, S.E. Schwartz (Eds.), Proc.
SEM/RILEM Int Conf. on Fracture of Concrete and Rock, Springer Verlag, New York, Houston, 1989, pp. 142–157.
[37] J.G. Rots, R. De Borst, Analysis of mixed-mode fracture in concrete, J. Engrg. Mech. 113 (11) (1987) 1739–1758.
[38] Z.P. Bazant, L. Cedolin, Blunt crack band propagation in finite element analysis, J. Engrg. Mech. (ASCE) (105) (1979) 297–315.
[39] G.I. Barenblatt, The mathematical theory of equilibrium of cracks in brittle fracture, Adv. Appl. Mech. 7 (1962) 55–129.
[40] D.S. Dugdale, Yielding of steel sheets containing slits, J. Mech. Phys. Solids 8 (1960) 100–108.
[41] A. Hillerborg, M. Modeer, P.-E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture
mechanics and finite elements, Cem. Concr. Res. 6 (1976) 773–782.
[42] J.C. Simo, J.G. Kennedy, S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions
and numerical algorithms, Int. J. Numer. Methods Engrg. 26 (1988) 2161–2185.
[43] D. Peric, D.R.J. Owen, Computational model for 3-D contact problems with friction based on the penalty method, Int. J. Numer.
Methods Engrg. 35 (1992) 1289–1309.
[44] D. Peric, E.A. De Souza Neto, A new computational model for Tresca plasticity at finite strains with an optimal parameterisation
in the principal space, Comput. Methods Appl. Mech. Engrg. (171) (1999) 463–489.
[45] N.J. Carpenter, R.L. Taylor, M.G. Katona, Lagrange constraints for transient finite element surface contact, Int. J. Numer.
Methods Engrg. 32 (1991) 103–128.
[46] Z. Zhong, L. Nilsson, Lagrange multiplier approach for evaluation of friction in explicit finite-element analysis, Comm. Numer.
Methods Engrg. 10 (1994) 249–255.
[47] T. Dede, Fracture onset and propagation in layered media, M.Sc. dissertation, University of the Witwaterstrand, Johannesburg,
South Africa, 1997.
[48] E. Hoek, E.T. Brown, Underground Excavations in Rock, Inst. of Mining and Metallurgy, London, 1980.
[49] B.J. Carter, E.Z. Lajtai, A. Petukhov, Primary and remote fracture around underground cavities, Int. J. Numer. Anal. Meth.
Geomech. 15 (1991) 21–40.
[50] B.J. Carter, E.Z. Lajtai, Y. Yuan, Tensile fracture from circular cavities loaded in compression, Int. J. Fract. (57) (1992) 221–236.
[51] E.J. Sellers, P.A. Klerck, Modelling of the effect of discontinuities on the extent of the fracture zone surrounding deep tunnels
under pressure, Tunnelling Underground Space Technol. 15 (2001) 463–469.
[52] M. Lee, B. Haimson, Laboratory study of borehole breakouts in Lac du Bonnet granite: a case of extensile failure mechanism, Int.
J. Mech. Min. Sci. Geomech. Abstr. 30 (7) (1993) 1039–1045.
[53] R.T. Ewy, N.G.W. Cook, Deformation and fracture around cylindrical openings in rock––II. Initiation, growth and interaction of
fractures, Int. J. Mech. Min. Sci. Geomech. Abstr. 27 (5) (1990) 409–427.
[54] Z. Zheng, J. Kemeny, N.G.W. Cook, Analysis of borehole breakouts, J. Geophys. Res. 94 (B6) (1989) 7171–7182.
[55] B. Haimson, I. Song, Laboratory study of borehole breakouts in cordova cream: a case of shear failure mechanism, Int. J. Mech.
Min. Sci. Geomech. Abstr. 30 (7) (1993) 1047–1056.
[56] B.J. Carter, Size and stress gradient effects on fracture around cavities, Rock Mech. Rock Engrg. 25 (3) (1992) 167–186.
[57] E.Z. Lajtai, Effect of tensile stress gradient on brittle fracture initiation, Int. J. Rock Mech. Min. Sci. 9 (1972) 569–578.
[58] C.D. Martin, J.B. Martino, E.J. Dzik, Comparison of borehole breakouts from laboratory and field tests, Eurock ’94, Balkema,
Rotterdam, 1994, pp. 183–190. ISBN 90-5410-502-X.
[59] C.D. Martin, Seventeenth Canadian Geotechnical Colloquium: the effect of cohesion loss and stress path on brittle rock strength,
Can. Geotech. J. (34) (1997) 698–725.
[60] R.T. Ewy, N.G.W. Cook, Deformation and fracture around cylindrical openings in rock––I. Observations and analysis of
deformations, Int. J. Mech. Min. Sci. Geomech. Abstr. 27 (5) (1990) 387–407.
[61] L.P. Gonano, Stress gradient and size effect phenomena in brittle materials, Ph.D. dissertation, James Cook University,
Townsville, Australia, 1974.
[62] F.J. Santarelli, E.T. Brown, Failure of three sedimentary rocks in triaxial and hollow cylinder compression tests, Int. J. Rock
Mech. Min. Sci. Geomech. Abstr. 26 (1989) 401–413.

View publication stats