An application of a cohesive fracture model combining compression, tension and shear in soft rocks

Computers and Geotechnics 66 (2015) 142–157

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

An application of a cohesive fracture model combining compression, tension and shear in soft rocks Y. Gui, H.H. Bui, J. Kodikara ⇑ Department of Civil Engineering, Monash University, Melbourne 3800, Australia

a r t i c l e

i n f o

Article history: Received 22 September 2014 Received in revised form 27 December 2014 Accepted 24 January 2015

Keywords: Cohesive fracture model Mixed-mode Tensile failure UDEC Uniaxial test Brazilian disc test

a b s t r a c t A cohesive fracture model considering tension, compression and shear material behaviour is implemented into the hybrid continuum–discrete element method, i.e. Universal Distinct Element Code (UDEC), to simulate possible multiple fracture in soft rocks. The fracture model considers both elastic and inelastic (decomposed to fracture and plastic) displacements. The norm of the effective inelastic displacement is used to control the fracture behaviour. Three numerical examples, including Mode-I, ModeII and Mixed-mode tests, are conducted to verify the model and its implementation in UDEC. The model is subsequently applied to simulate uniaxial compression and Brazilian disc tests on soft rocks and the results are compared with experimental results. The results indicate that the cohesive fracture model is capable of realistically simulating the combined tensile, shear and mixed-mode failure behaviour applicable to geomaterials. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction It is generally accepted that the linear elastic fracture mechanics (LEFM) is a useful approach for addressing fracture problems involving cracks where the non-linear process zone in front of the crack is sufficiently small [1]. In addition, LEFM assumes that the intact material behaviour is linear elastic. However, for many geomaterials, such as soils, rocks, concrete, cement-stabilised rock aggregates and soil, it may be unrealistic to consider the size of the non-linear zone to be negligible and intact materials to be linear elastic [1,2]. To overcome these shortcomings, a cohesive fracture model (also referred to as a fictitious crack model) which was first proposed by Dugdale [3], has been advanced [1,2,4–7]. In cohesive fracture model, both strength and fracture criteria are combined together offering advantages over models based on LEFM [8]. In Mode-I fracture, the cohesive fracture model and its constitutive behaviour may be described as shown in Fig. 1. The fracture is considered to consist of two components: a real crack and a fictitious crack (also known as the process zone), each of which is associated with a crack tip, namely the real tip and the virtual crack tip, respectively. The process zone is then defined as the zone between the real and fictitious crack tips, and comprises the material that is partially damaged but is still able to transfer load across the fracture [4]. The crack opening behaviour is governed by the value ⇑ Corresponding author. Tel.: +61 3 9905 4963. E-mail address: [email protected] (J. Kodikara). http://dx.doi.org/10.1016/j.compgeo.2015.01.018 0266-352X/Ó 2015 Elsevier Ltd. All rights reserved.

of the opening displacement and the strength of the material. The initial hardening behaviour is assumed to be linear elastic when the opening displacement is smaller than a critical value (i.e. wc1 in Fig. 1b). The crack surface traction at the critical opening displacement is the material tensile strength (i.e. rt in Fig. 1b). For openings larger than wc1 , the bridging stress across the fracture will decrease featuring softening behaviour, and will become zero at a limiting displacement (i.e. wc2 in Fig. 1b). Generally, the initial hardening response is relatively small in comparison to the softening response, thus softening response has therefore received more attention in the past [1]. To explain the softening behaviour, several softening laws have been proposed, including mono-linear, bi-linear and other softening laws as shown in Fig. 1. Although there have been some advancements in cohesive fracture modelling, for example [1,2,4–7], only few models have considered Mixed-mode fracture and most are still focused on Mode-I case in geomaterials. However, some experimental and numerical evidences have shown that cohesion-softening (or sometimes referred to as decohesion) is necessary when considering the plastic/frictional dissipation in soft rocks. For example, Vermeer and de Borst [9] indicated that cohesion-softening is to be expected due to micro-fracturing and the degradation of bond between grains when material undergoing plastic deformation. Edelbro [10] also demonstrated the applicability of cohesion-softening in simulating rock mass behaviours. Therefore, it is rational to take Mix-mode cohesive fracture model in describing fracturing behaviour of soft rocks.

143

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

Most existing cohesive models in the literature are based on damage theory and do not take into account the plastic/frictional dissipation [8]. In addition, it has been realised that although pure damage model is more efficient in terms of computation than plasticity/damage coupled model, coupled plasticity/damage cohesive model can have true predictive capability as friction and fracture are both dissipative mechanisms which contribute to interface damage process [8]. In tensile loading, the plastic deformation during loading–unloading is contributed from two aspects: (1) the rock fracture cannot be fully healed (i.e. there is an increase of porosity locally) due to the mismatch of the interface and the loss of rock grains at the rock fracture, especially, the rock fracture is under shearing. (2) There is cohesive zone in front of real crack tip which contributes to the plastic deformation. This plastic deformation cannot be reversed even if the load on the fracture has been removed. In fact, some experiments on concrete (e.g. [11]), soft rock (e.g. [12]), and soil (e.g. [13]) have demonstrated the existence of the plasticity during loading–unloading through bending tests. Kazerani et al. [5] proposed a cohesive fracture model considering the cohesive behaviour of tensile, compressive-shear in rocks. First, the hardening stage of the load–displacement response was simulated using an exponential law instead of a linear law adopted in most cohesive models. In fact, there has been no strong experimental evidence that can show the initial elastic behaviour is exponential. In addition, the plasticity dissipation was not considered in this model. The unloading/reloading path was treated as elastic, where the unloading/reloading path went through the origin of stress-displacement space (i.e. state with zero displacement and zero stress), assuming that the cohesive process zone behaves fully elastically. Moreover, there is no cohesive effect actually included in the compressive-shear behaviour as the cohesion decreases to a constant residual value immediately once the initial strength has been reached. In this paper, a cohesive fracture model originally presented in Ref. [17] that takes into account tensile, shear and compressive behaviour combined with an evolutionary failure model applicable to general Mixed-mode rock fractures is investigated using a hybrid continuum–discrete element method. Implementation of

b

σt

wc2

real crack tip

In this section, the hybrid continuum–discrete element method is presented. The method can model the interaction of discrete materials, the solid continuum within discrete bodies and the fracturing process. A distinct element numerical scheme needs to allow for finite displacements and rotations of the discrete bodies,

c

σt

σt

traction, σ

process zone intact zone

2. Hybrid continuum–discrete element modelling

traction, σ

a

the model in the hybrid method has the advantage of handling multiple fracture and deformation problems in materials over other methods, such as finite element method (FEM) and discrete element method (DEM). More specifically, FEM has problem of handling multiple fractures. In conventional DEM, the block deformation is not possible to be considered. The original work on this model was developed for applications in concrete primarily for single crack simulations. However, our implementations catered for multiple crack formation as needed for rock and soil fragmentation where both shear and tension modes can be prevalent. In this model, elastic–plastic-damage is coupled together. Therefore, the dissipative mechanism (i.e. plasticity during unloading–reloading) has been taken into account in the model. In addition, the definition of elastic stiffness is based on fracture energy which eases the selection of parameters. The applicability of the model is tested through the verification of the implementation and then applying to real soft rock experimental data. For other cohesive fracture models, the implementation can be time-consuming. LEFM cannot easily simulate the development of multiple interacting cracks and disintegration of elements. Therefore, in this paper, the verification of implementation is only tested through the simulation result with the data obtained from model equations. The paper is organised as follows: in Section 2, the 2-dimensional hybrid continuum–discrete element model is introduced; the cohesive model framework is illlustrated in Section 3. In Section 4, the cohesive model is implemented and verified using uniaxial tension and compressive/shear tests. In addition, the model is applied to simulate geomechanical test examples, including uniaxial compression and Brazilian disc tests.

virtual crack tip wc1

wc1

wc 2

crack opening w

e

f

σt

σt

traction, σ

traction, σ

σt

traction, σ

d

wc1

wc 2

crack opening w

wc 2

crack opening w

wc1

wc 2

crack opening w Fig. 1. Mode-I cohesive fracture model and its consititutive behaviour.

wc1

wc 2

crack opening w

144

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

including the possibility of complete detachment and regaining new contacts automatically as the deformation progresses [14]. The material mass is represented by an assemblage of discrete blocks [15], each with their own continuum behaviour. The fractures or potential fractures in the materials are simulated by interface elements between blocks. The contact forces and displacements across these interfaces are calculated and the movements of the blocks are traced. The features of the hybrid continuum–discrete element method, i.e. Universal Distinct Element Code (UDEC) are elaborated in the following five sub-sections.

R €i ¼ u

rij nj ds þ F i

s

m

þ gi

ð7Þ

where rij is the element stress, s is the surface enclosing the mass m, which is considered to concentrate at the node, and nj is the unit normal to s. F i and g i are the resultants of all external forces applied to the node and gravitational acceleration, respectively. The nodal forces are calculated by summing all external applied forces, contact forces (applying only to the nodes along the block boundary) and the contribution from internal stresses in elements adjacent to the node, as given below:

2.1. Equations of block motion

F i ¼ F li þ F ci þ F iz

The equations of motion describing the translation and rotation of the block about its centroid are developed [15]. In one dimension, Newton’s second law of motion can be written in the following form:

where F li is the external applied forces, F ci the contact forces and F iz the contribution from internal stress in the elements adjacent to the node (Fig. 2b). The contribution from internal stresses in the elements adjacent to the node is calculated as:

du_ F t ¼ þg dt m

F iz ¼

ð1Þ

_ t and m are the velocity, time and mass, respectively. g is where u; the gravitational acceleration. The central difference scheme for the left-hand side of Eq. (1) can be written as:

du_ u_ tþDt=2  u_ tDt=2 ¼ dt Dt

ð2Þ

Substituting Eq. (1) into Eq. (2) yields:

u_ ðtþDt=2Þ ¼ u_ ðtDt=2Þ þ



 Ft þ g Dt m

ð3Þ

With velocities stored at the half-time step, it is possible to express displacement as:

uðtþDtÞ ¼ uðtÞ þ u_ ðtþDt=2Þ Dt

ð4Þ

As the force depends on the displacement, the force/displacement calculation is performed at one time-step. In two dimensions, blocks can be loaded by multi-directional forces including gravity. Therefore, the velocity equations become:

8  t  RF > < u_ iðtþDt=2Þ ¼ u_ iðtDt=2Þ þ mi þ g i Dt  t > : h_ ðtþDt=2Þ ¼ h_ ðtDt=2Þ þ RMi Dt i i I

ðtþDtÞ

xi

ðtþDtÞ

hi

ðtÞ ðtþDt=2Þ ¼ xi þ u_ i Dt ðtÞ ðtþDt=2Þ ¼ h þ h_ Dt i

rij nj ds

ð9Þ

c

where nj is the unit outward normal to the contour c. The contour c follows the closed polygonal line defined by the straight segments that bisect the element edges converging on the node under consideration. The stress in Eq. (9) is calculated using a constitutive relation adopted for deformable blocks. Although it is not essential, in this paper, blocks are assumed to be linear elastic and their constitutive response can be described using the following stress–strain incremental form so that implementation on nonlinear problems can be accomplished easily:

Dreij ¼ kl Dev dij þ 2ll Dei;j

ð10Þ

where kl ; ll are the Lamé constants, Dr Dei;j and Dev are the incremental forms of elastic stress, strains and volumetric strain. Dev ¼ De11 þ De22 þ De33 (for plane strain De33 ¼ 0Þ. dij is the wellknown Kronecker delta function. During each time–step, the strain and rotation are calculated based on nodal displacements as: e ij ;

 1 u_ i;j þ u_ j;i 2  1 x_ ij ¼ u_ i;j  u_ j;i 2

e_ ij ¼ ð5Þ

where h_ i ; I; RMti ; u_ i and g i are respectively the angular velocity of the block about centriod, the moment of inertia of the block, the total moment acting on the block, the velocity components of block centroid and the component of gravitational acceleration (body force). Finally, based on the velocity equations, the block location can be determined according to the following equations:

(

Z

ð8Þ

ð6Þ

i

where hi and xi are respectively the rotation angle of a particle about the centroid and its coordinates. 2.2. Constitutive model of blocks Deformable blocks are discretised into finite-difference triangular elements (Fig. 2a). The use of these triangular elements is to calculate block deformation rather than fracture propagation. UDEC uses mixed discretization technique to ensure that plastic collapse and flow accurately modelled under Lagrangian analysis [16]. The computational accuracy depends on the number of elements contained in each block. For each triangular element, there are three nodes and the equation of motion for each node is described as:

ð11Þ

For large deformation treatment allowing for rigid body rotation, Jaumann stress rate was adopted to Eq. (10). This stress rate _ jk  x _ ik  rkj with r _ Jaumann ^ ij ¼ r_ ij  rik  x ^ ; r; x is calculated as r stress rate, total stress and spin-rate, respectively. Once all the forces are defined, each node of the element is accelerated according to the finite-difference form of Newton’s second law of motion. ðtþDt=2Þ ðtDt=2Þ ¼ u_ i þ u_ i

RF iðtÞ m

! þ g i Dt

ð12Þ

2.3. Constitutive law of contact In the current numerical method, the crack or flaw in the material (or potential crack or flaw) is represented numerically using a contact surface (i.e. interface element). The constitutive law of contact is based on the techniques of contact detection during calculation. Fig. 3(a) shows the contact logic used in UDEC. Two blocks (i.e. blocks 1 and 2) are considered. The contact interfaces are made up of three pairs of contact points (i.e. A1  A2 ; B1  B2 ; C 1  C 2 Þ and their related block boundary lengths. For each pair of contact logic points, there are two spring bonds which are responsible for normal and shear displacements at the contact. The normal and shear

145

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

I

a

H

J

G 13 12

16

14 15

P

F z11

b

1

Fl

L

17

O

2

F

K 18

9

10

8

M

11

m

Fc

F z2

7 N

L

E 2

3

4

6

3

1

F z1

5

4

A

D

A

s

mg

B

B

C

Fig. 2. Schematic illustration of block: (a) block discretisation and (b) forces applied on node.

a

b block 1 A1

block 1 A1 k nA

k sA

A2

B1 k nB

C1

k sB

B2

k nC

B1

C1

ΔunA Δ u sA

k sC

A2

C2

B2

C2

block 2

block 2

Fig. 3. Schematic illustration of contact logic in UDEC (a) the mechanical model of contact logic and (b) the relative displacement detection logic. DunA and DusA are relative normal and shear dispalcement, respectively.

stiffness of the contact logic are kn and ks . Any form of displacement in one of the blocks can cause deformation of springs at contact, and thus force will be generated in the springs. Meanwhile, the force on the springs will be exerted on the neighbouring block through the contact point pair. The calculation of relative displacement between two blocks is illustrated in Fig. 3(b). Next, the stress calculation based on the relative displacement is formulated. In the normal direction of the contact interface, the incremental form of the stress-displacement relation is assumed to be linear and governed by normal stiffness (kn Þ as:

Drn ¼ kn Dun

Ds ¼

Dtn ¼ 2 min

ð15Þ

where Dues is the elastic component of the incremental shear displacement.

 1=2 mi ki

ð16Þ

where mi and ki are respectively the mass associated with block node i and the stiffness of the elements surrounding the node. The stiffness is calculated, considering both interface and the continuum of intact block, as:

  ki ¼ R kzi þ kji

ð17Þ

In Eq. (17) Rkzi represents the contribution of the stiffness of all the zones connected to node i. The kzi is obtained using P-wave modulus which is defined as the ratio of axial stress to axial strain [15]:

ð14Þ

If the shear stress is smaller than the ultimate value (i.e. smax Þ, the incremental form of shear stress-shear displacement relations can be written down as:

ks Dues

The solution scheme used for the distinct element method is conditionally stable. Therefore, it is necessary to determine a limiting (psuedo) time-step that satisfies the stability criterion for the calculation of internal block deformation as well as the inter-block relative displacement. The limiting time-step for the stability of block deformation computation is estimated according to [15]:

ð13Þ

There is a limiting tensile strength, rt for a contact. In any cases when the tensile strength is exceeded by the tensile stress (i.e. if rn > rt Þ, then its value is set to zero (i.e. rn ¼ 0Þ. The contact may undergo compression if two blocks overlap each other. For the case of shear, similarly, the response is controlled by a shear stiffness, ks . The shear stress, s, is limited by a combination of cohesive and frictional strength through the contact cohesion, c, and friction angle, u, respectively. Thus, the ultimate shear stress is calculated as:

smax ¼ c þ rn tan u

2.4. Time-step determination

kzi ¼

  2 8 4 b K þ G max 3 3 hmin

ð18Þ

with K and G representing respectively the bulk and shear moduli of the block material, while bmax and hmin stand for the largest zone size and the minimum height of the triangular zone. The joint stiffness term, kji , exists only for nodes located on the block boundary, and is taken as the product of the normal or shear joint stiffnesses

146

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

(whichever is larger) and the sum of the lengths of the two block edge segments adjacent to node i. As for the calculation of inter-block relative displacement, the limiting time-step is calculated by analogy to a single degree-offreedom system with the mass of the smallest block in the system (i.e. M min Þ and the maximum contact stiffness (i.e. kmax Þ as:

Dtb ¼ 2ðfracÞ

(6) under compressive normal stresses, neither the shear nor the normal stiffness decreases to zero and; (7) dilation angle decreases with sliding or opening displacement. In addition, irreversible relative displacements are caused by broken segments of the interface material and by friction between the two fracture surfaces. 3.2. The constitutive model of the cohesive fracture model

 1=2 Mmin kmax

ð19Þ

where frac is a user-supplied value that accounts for the fact that a single block may be in contact with several blocks simultaneously. A typical value of frac is 0.1. Based on the above two limiting time-steps, the chosen timestep is the smaller one between the two time-steps, i.e.

Dt ¼ min ðDt n ; Dt b Þ

ð20Þ

When a material is undergoing loading, its displacement (i.e. u) can be partitioned into elastic (i.e. ue Þ and inelastic parts (i.e. ui Þ as:

u ¼ ue þ ui

ð22Þ

The inelastic displacement can be further decomposed into plastic displacement (i.e. up Þ, which is irreversible, and fracture displacement (i.e. u f Þ as:

ui ¼ up þ u f

ð23Þ

The norm of the inelastic displacement is defined as:

2.5. Damping Mechanical damping is introduced in the distinct element method to make the model tolerant of numerical instability. In this simulation, the damping is proportional to the magnitude of the unbalanced force. The damping force’s direction is such that energy is always dissipated. For deformable blocks, Eq. (12) is replaced by the following equation incorporating local damping:

  n   o Dt  ðtÞ  ðtÞ u_ ðtþDt=2Þ ¼ u_ ðtDt=2Þ þ RF i  HRF i sign u_ ðtDt=2Þ mn

ð21Þ

where H is a constant set to be 0.8 in UDEC [15] and found to be suitable for the following numerical validations and applications and mn is the nodal mass.

uieff ¼ jjui jj ¼ jjup þ u f jj ¼

rt



8   < r 1  uieff  uieff < wr t0 wr ieff u ¼ : 0 uieff P wr

The fracture interfaces are idealised with zero thickness and the major features of the model are as follows: (1) shear strength is affected by the normal stress; (2) softening occurs both in shear and tension; (3) there is residual shear strength due to the friction along the interface; (4) strength degradation is caused by fracture propagation; (5) there is zero normal and shear stiffness when the interface is fully destroyed under tension or full fracture opening;

aσ

2GIf

ð26Þ

rt0

In Eqs. (25) and (26), wr ; rt0 and GIf are the ultimate norm of the inelastic displacement corresponding to zero tensile strength, the initial tensile strength and the Mode-I fracture energy, respectively. The ultimate norm of the inelastic displacement corresponds to the threshold condition where the crack is fully developed and the material is no longer capable of transferring stress. The initial tensile strength is the stress at which the cohesive zone starts to develop and the crack starts to undergo softening. The tensile strength evolution during crack opening is also illustrated in Fig. 4(a). In order to cater for stiffness degradation during softening, a micro-damage variable D is introduced as the percentage of fracture surface (i.e. Af Þ to the overall surface area of an interface (i.e. A0 Þ. D can be calculated as:

D¼

Af kns ¼1 A0 kn0

ð27Þ

b

c

t

σt0

c0 I

G IIa f

Gf

0

ð25Þ

and

3. Cohesive fracture model

3.1. General features of the model

wσ

u ieff

ð24Þ

In tensile loading, the governing variables are the tensile strength (rt Þ and the norm of the inelastic displacement (uieff Þ. The tensile strength is a linear function of the norm of the inelastic displacement given as:

wr ¼

The cohesive fracture model presented herein is based on the framework of the interface elastic–plastic-damage constitutive model presented in [17] that takes into account the cohesive effect on both tension and shear. The hardening stage is assumed to be linear elastic, and the emphasis is placed on the treatment of the softening stage.

r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 2 uin þ uis

0

Fig. 4. (a) Tensile and (b) shear softening law of the cohesive fracture model.

wc

uieff

147

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

where kns and kn0 are respectively the degraded and initial normal stiffness. kns can be computed as:

kns ¼ ¼

rn un  upn

¼

rt ðuieff Þ uen þ upn þ unf  upn

rt ðuieff Þ rt ðuieff Þ=kn0 þ ð1  gÞuieff

ð28Þ

where g is the ratio of irreversible inelastic normal displacement to the total value of inelastic displacement (i.e. g ¼ up =ui with ui the norm of inelastic displacement before unloading) and can be determined experimentaly using the pure Mode-I test. Substituting Eq. (28) into Eq. (27), the micro-damage variable D is expressed based on the norm of inelastic displacement as:

D¼1

rt

ðuieff Þ

rt ðuieff Þ þ ð1  gÞuieff kn0

ð29Þ

Accordingly, the normal stress-displacement relationship, i.e. the relationship between rn and ðun  upn Þ of the interface, is presented as:









3.3. Failure function A general failure function is needed to cater for the combined influence of normal and shear stresses. The normal stresses can be either tensile or compressive. As described earlier, with the variation of the norm of the inelastic displacement, the material tensile strength and cohesion are changed. In this paper, a modified non-linear form of the Mohr–Coulomb failure criterion is adopted to represent the initial failure states, as shown in Fig. 5. As the fracture develops, the tensile strength and cohesion will evolve, shrinking the failure surface to the final failure surface, which is the conventional Mohr–Coulomb failure surface with zero tensile strength and zero cohesion. Mathematically, the failure surface function can be expressed as:

  f ¼ s2  2c tan ðuÞðrt  rÞ  tan2 ðuÞ r2  r2t ¼ 0

ð36Þ

where u is the friction angle which is kept unchanged during the calculation, while rt and c are evolved as per Eqs. (25) and (32), respectively.

ð30Þ

3.4. Numerical implementation of the cohesive fracture model in UDEC

where parameter a is the integrity parameter defining the relative active area of the fracture and is related to the micro-damage variable D through the following relation:

For simplicity of description, the superscripts N; N þ 1 in the following illustration indicate variables at step N, and N þ 1, respectively. Superscript zero indicates the initial variables used in the simulation. The subscipts n and s indicate variables corresponding to the direction perpendicular and along the cohesive fracture considered. At step N, the relative normal and shear displacement are uNn

rn ¼ kns un  upn ¼ akno un  upn

a¼1

jrn j þ rn D 2jrn j

ð31Þ

The intergrity parameter is used to consider tensile unloading– reloading behaviour in the cohesive fracture model. It can be seen from Eq. (31) that the activation of the micro-damage variable is controlled through the fraction normal stress, which is activated in tension (i.e. rn > 0Þ and deactivated in compression (i.e. rn < 0Þ. Thus, both the normal and shear stiffness in tension can be degraded, while they are kept unchanged in compression. In a similar fashion, the governing variables for the shear loading are the cohesion c (the contribution of normal stress to shear strength can be neglected if the friction angle is taken to be zero) and the norm of the inelastic displacement uieff . As the crack propagates, the cohesion is degraded and can be expressed as a linear function of the norm of the inelastic displacement uieff as:

8   < c 1  uieff  uieff < wc 0 wc c uieff ¼ : 0 uieff P wc 

ð32Þ

and

2GIIa f wc ¼ c0

ð33Þ

In Eqs. (32) and (33), wc is the ultimate norm of inelastic displacement corresponding to zero cohesion (see Fig. 4b), c0 is the initial cohesion and GIIa is the Mode-II fracture energy dissipated during f shear at high confining normal stress (i.e. without the influence of the tensile loading regime). The ultimate norm of inelastic displacement corresponding to zero cohesion gives the threshold condition at which the shear crack is fully developed and the material is no longer capable of transfering cohesion. Similar to the softening treatment under tensile loading, the degraded shear stiffness (i.e. kss Þ can be written as:

kss ¼ aks0

ð34Þ

The shear stress (sÞ is computed by:







s ¼ kss us  ups ¼ akso us  ups



rn with a ¼ 1  jr2njjþ rn j D, and kso the initial shear stiffness.

ð35Þ

N

N

and uNS , the normal and shear stiffness are kn and ks , the micro N

damage variable and integrity parameter are D and aN , the norm of inelastic displacement, tensile strength and cohesion, respectively, are uNieff ; rNt and cN . The normal and shear stress is rNn and

sN . The dilation angle is uNd . At step N þ 1, the relative normal and shear displacement are uNþ1 and uNþ1 , thus the incremental form of the relative normal n S and shear displacement can be expressed as:

Dun ¼ uNþ1  uNn n

ð37Þ

Dus ¼ uNþ1  uNs s

where Dun and Dus are the increment of relative normal and shear displacement. Based on the incremental form of the relation between stress and displacement in UDEC, the new stress at step N þ 1 can be calculated as: N rNþ1 ¼ rNn þ kn Dun n N sNþ1 ¼ sN þ ks Dus

ð38Þ

Substituting Eq. (38) into the failure function (i.e. Eq. (36)), produces: 2

F ¼ ðsNþ1 Þ  2cN tanðuÞ







2



rNt  rNþ1  rNt  tan2 ðuÞ rNþ1 n n

2  ð39Þ

If F 6 0, the fracture surface behaviour is under the elastic region and no intervention needs to be imposed into UDEC. However, in the case of F > 0, the fracture surface has entered the plastic region and any departure from the failure surface is corrected by a return-mapping algorithm. An elastic predictor is used and the normal and shear stresses at step N þ 1 are recalculated using the elastic relation and initial stiffnesses as:



e

0 rNþ1 ¼ rNn þ kn Dun n e 0 ðsNþ1 Þ ¼ sN þ ks Dus

ð40Þ

148

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

 e e where rNþ1 and ðsNþ1 Þ are elastic predictors. As the elastic pren dictors at step N þ 1 in Eq. (40) are in the plastic range, a return step is needed to correct the stress back to the failure surface as:



e

0 rNþ1 ¼ rNþ1  kkn mn n n e sNþ1 ¼ ðsNþ1 Þ  kk0s ms

ð41Þ

Therefore, the new expresssion of the failure surface after stress returning is:

   2 F ¼ ðsNþ1  kms Þ  2cN tanðu0 Þ rNt  rNþ1  kmn n  2  2   kmn  rNt  tan2 ðu0 Þ rNþ1 n

ð42Þ

In Eqs. (41) and (42), k is the inelastic multiplier and m defines the direction of the inelastic displacement. Determination of the direction of the inelastic displacement is based on the stress state defined in Eq. (40) and the failure surface as shown in Fig. 6. The stress space is divided into two domains, segmented by the line OP 1 . The return direction in space rOP 1 is along the line OS1 conand sNþ1 (i.e. point S1 Þ necting the stress space position of rNþ1 n and the origin O. The return direction in the rest of the space is defined by the line P2 S2 which is orthoganal to the current plastic potential line. The plastic potential is defined as:

!

N

kn

Q ¼s

N

ks

tan u

N d

r¼0

ð43Þ

Therefore, the direction of the inelastic displacements can be computed as:

( mn ¼ (

0 =kn

r in space rOP1 and sNþ1 tan uNd in rest of the space Nþ1

sNþ1 =k0s in space rOP1 ms ¼ sNþ1 in rest of the space

Based on Eqs. (44) and (42) and making Eq. (42) equal zero, the inelastic multiplier (kÞ can be determined. From the inelastic multiplier and direction of the inelastic displacement, the norm of inelastic displacement is updated as:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N uNþ1 m2n þ m2s ieff ¼ uieff þ k

ð45Þ

Based on the updated norm of the inelastic displacement, the tensile strength, cohesion, and dilation angle at step N þ 1 are updated as:

rNþ1 ¼ r0t 1  t Nþ1 d

u

¼u

0 d

uNþ1 ieff

1

! ;

wt uNþ1 ieff

cNþ1 ¼ c0 1 

!

uNþ1 ieff wc

! ; ð46Þ

wd

The normal and shear stresses are updated for step N þ 1 using Eq. (40) together with the direction of inelastic displacement and inelastic multiplier calculated from Eqs. (42) and (44) as:



e

0 rNþ1 ¼ rNþ1  kkn mn n n e 0 sNþ1 ¼ ðsNþ1 Þ  kks ms

ð47Þ

The micro-damage variable and the integrity paremeter are expressed as:

DNþ1 ¼ 1 

rNþ1 t

; 0 rNþ1 þ ð1  kÞuNþ1 t ieff kn

aNþ1 ¼ 1 

 Nþ1  r  þ rNþ1 n  n DNþ1 2rNþ1  n

ð48Þ Therefore, the degraded stiffness are calculated as: Nþ1

kn ð44Þ

0

Nþ1

¼ aNþ1 kn ; ks

0

¼ aNþ1 ks

ð49Þ

The flow chart of the implementation is shown in Fig. 7. 4. Model validation and applications

τ

c

The numerical implementation of the cohesive fracture model was programmed using FISH language in UDEC which was introduced in Section 2. In the following sections, three numerical examples, including Mode-I, Mode-II and Mixed-mode tests, are first conducted to verify the cohesive fracture model. Then, the proposed model is applied to simulate rock uniaxial compression and Brazilian disc tests as the preliminary applications of the cohesive fracture model.

initial failure surface

final failure surface

0 compression range

4.1. Mode-I test

σ

σt

The purpose of this test is to verify the implementation of the tensile and compressive behaviour of the cohesive fracture model in UDEC. Transjurane sandstone, previously utilised by Kazerani et al. [18,19], is used as the material for verification. The input parameters for the cohesive fracture model are the Young’s modulus, Poisson’s ratio, initial cohesion, intial tensile strength, friction angle and initial normal and shear stiffness, as well as initial dilation angle, which are summarised in Table 1. The normal and shear stiffness are calculated from the fracture energy which is related to the fracture toughness. In fracture mechanics, the fracture energy is expressed as:

tension range

Fig. 5. The evolution of failure surface.

τ P1 S2

S1

c

GIf ¼

Q=0

P2

o

σt

Fig. 6. Schematic illustration of the inelastic return direction.

σ

K 2IC E0

and GIIf ¼

K 2IIC E0

ð50Þ

where GIf is the Mode-I fracture energy, K IC the Mode-I fracture toughness, GIIf the Mode-II fracture energy and K IIC the Mode-II fracture toughness. E0 is calculated as:

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

E0 ¼



E

plane-stress

E=ð1  t2 Þ plane-strain

ð51Þ

where E the Young’s modulus and t the Poisson’s ratio. The ultimate norm of inelastic displacement corresponding to zero tensile strengh and zero cohesion, i.e. wr and wc can then be calculated from Fig. 4 as:

wr ¼

2GIf

rt0

and wc ¼

2GIIa f c0

ð52Þ

149

Table 1 Mechnical properties of Transjurane sandstone. Young’s modulus, E (GPa) Poisson’s ratio, t pffiffiffiffiffi Fracture toughness in Mode I, K IC ðMPa= mÞ pffiffiffiffiffi Fracture toughness in Mode II, K IIC ðMPa= mÞ Cohesion, c0 (MPa) Tensile strength, rt0 (MPa) Friction angle, u (degree) Dilation angle, ud0 (degree)

Fig. 7. Flow chart of the algorithm implemented in UDEC to calculate contact stress.

12.5 0.3 0.7 0.8 8.5 2.8 41 10

150

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

It should be noted that GIIa in Eq. (52) is not the pure Mode-II f fracture energy

GIIf .

However, due to the fact that the pure Mode-

II fracture energy is close to GIIa f , the pure Mode-II fracture energy is used here to approximately calculate wc .

vy

a

vy

b

The contact initial stiffness was evaluated on the basis of material mechanical properties and thickness of cohesive zone by Kazerani et al. [5] from material strength in molecular mechanics and material fracturing. They are calculated as:

8 2 > E rt0 E < 4K 2IC kn ¼ ¼ 2 l > : E rt0

plane-stress and plane-strain

4ð1t2 ÞK 2IC

8 EGr t0 G < 4K 2IIC ks ¼ ¼ : EGrt0 2 l 2

plane-stress

100 mm

100 mm

Interface

Interface

where, l is the cohesive zone thickness perpendicular to the orientation of the fracture propagation, and G is the elastic shear modulus which can be calculated as:

G¼

E 2ð1 þ tÞ

ð54Þ

The cohesive zone thickness may be estimated using material mechanical properties based on the theory of material strength in molecular mechanics and material fracturing as [5]:

l¼

8 I > < 4Gf

Pure Mode I

II > : 4Gf

Pure Mode II

rt0 rt0

50 mm

ð53Þ

plane-strain

4ð1t ÞK IIC

ð55Þ

50 mm 3

Normal/Shear stress, MPa

3.5

Numerical implementaion Model equation

3 2.5

2.5

Tensile strength, MPa

Fig. 8. The numerical models and the boundary condition in Mode-I test: (a) uniaxial tension, (b) uniaxial compression.

Numerical implementation 2

Model equation

1.5

1

0.5

2 1.5

0 0

0.000005

1

0.00001

0.000015

0.00002

0.000025

0.00003

Norm of inelastic displacement, m

0.5

Fig. 10. Tensile strength with norm of the inelastic displacement in Mode-I tension case.

0 0

0.00001

0.00002

0.00003

0.00004

Micro- damage variable D and integrity parameter α

Normal/shear relative displacement, m

Normal/Shear stress, MPa

3 2.5 2

Numerical implementaion Model equation

1.5 1 0.5

1.2 1 0.8

0.4 0.2 0 - 0.2

-0.000005

0 0

0.000001

0.000002

0.000003

0.000004

0.000005

α (numerical implementation) D (numerical implementation) α (model equation) D (model equation)

0.6

0

0.000005 0.00001 0.000015 0.00002 0.000025 0.00003

Norm of inelastic displacement, m

Normal/shear relative displacement, m Fig. 9. Tensile stress with opening displacement in Mode-I tension case.

Fig. 11. Micro-damage variable (DÞ and integrity parameter (a) with norm of the inelastic displacement in Mode-I tension case.

151

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

2.4E+08 1.E+09 1.E+08

Kn, MPa/m

2.0E+08

Kn, MPa/m

1.6E+08

1.E+07 1.E+06 1.E+05 1.E+04 1.E+03

1.2E+08

0

0.00001 0.00002 Norm of inelastic displacement, m

0.00003

8.0E+07 Model equation Numerical implementation

4.0E+07

0.0E+00 -0.00001 -0.000005

0

0.000005

0.00001

0.000015

0.00002

0.000025

0.00003

0.000035

Norm of inelastic displacement, m Fig. 12. Normal stiffness (kn Þ with norm of the inelastic displacement in Mode-I tension case.

3.5E+05 3.0E+05

1

2.5E+05

Kn, MPa/m

0.8 2.0E+05 0.6 1.5E+05 0.4 1.0E+05 0.2

5.0E+04 0.0E+00

Micro damage variable, D

1.2

Normal stiffness (model equationl) Normal stiffness (numerical implementation) Micro damge variable (model equation) Micro damge variable (numerical implementation)

0

0.000005

0.00001

0.000015

0.00002

0.000025

0 0.00003

Norm of inelastic displacement, m Fig. 13. Normal stiffness (kn Þ and micro-damage variable (DÞ with norm of the inelastic displacement in Mode-I tension case (initial normal stiffness is 3  1011 Pa/m).

Fy 9 8

Interface

50 mm

vx

Numerical implementation

Shear stress, MPa

vx

7

Model equation

6 5 4 3 2 1

50 mm 0 0

0.000002 0.000004 0.000006 0.000008 0.00001 0.000012 0.000014

Shear relative displacement, m Fig. 14. The numerical models and the boundary conditions in Mode-II test.

Fig. 15. Shear stress with relative shear displacement in Mode-II test.

152

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

order to simulate the behaviour of the fracture under tension, a tensile load is applied on the top of the sample, while the bottom is fixed along vertical direction. However, the behaviour of the fracture under compression is also studied. In this case, a compressive load is imposed on the top. In the actual simulation, the force is accomplished by applying constant velocity at the top boundary. The simulation is terminated once the normal stress on the interface decreases to zero. The models are set as the plane-stress condition. During modelling, the normal stress (rn Þ along the interface and relative displacement (un Þ at the interface, tensile strength (rt Þ, micro-damage variable (DÞ, integrity parameter (a), and normal stiffness (kn Þ are tracked. The modelling of opening crack failure is based on a mono-linear softening model, as shown in Fig. 4(a). The simulated results are shown in Figs. 9–12. Fig. 9 gives the comparison of the simulation result from the computer code implementation of the analytical model and the calculation result from our anlytical model for the opening displacement and the tensile stress on the crack. The tensile stress initially increases with opening until it reaches the peak stress (2.80 MPa). Then the softening stage starts and the tensile stress on the fracture surface decreases with open-

9 8 7

Model equation

5 4 3 2 1 0 0

0.000005

0.00001

0.000015

Norm of inelastic displacement, m

Micro damage variable D and integrity parameter α

Fig. 16. Cohesion with the norm of inelastic displacement in Mode-II test

1.2

ing until zero when the opening is around 2:81  105 m. The simulated tensile strength with the norm of inelastic displacement is shown in Fig. 10. Fig. 11 shows the evolution of the micro-damage variable and the integrity parameter with the norm of the inelastic displacement. The changes of the micro-damage variable and integrity parameter are opposite and the micro-damage variable increases from 0 to 1, while the integrity parameter decreases from 1 to 0. Fig. 12 demonstrates the degradation of the normal stiffness with the norm of the inelastic displacement. Both linear and logarithmic scale are presented. To demonstrate the effect of initial

1 0.8 0.6 0.4 0.2

stiffness, an initial normal stiffness value of 3  1011 Pa=m is used in the simulation and its degradation is presented in Fig. 13. In this simulation, other parameters are the same as the previous one. As can be seen that all the simulation results are in good agreement with the values directly calculated from model equations, demonstrating the validity of the numerical implementaion of the cohesive fracture model. Then the Mode-I compression is studied. As under compressive loading, the normal stress is negative and never hits the failure surface. Then all relative displacement on the fracture is considered as elastic displacement. Since the fracture is intact and there is no fractured area (i.e. Af ¼ 0Þ, the micro-damage variable is zero and the integrity parameter is 1. Accordingly, there is no change of tensile strength and normal stiffness.

α (numerical implementation) D (numerical implementation) α (model equation) D (model equation)

0

-0.2 -0.000005

0

0.000005

0.00001

0.000015

Norm of inelastic displacement, m Fig. 17. Micro-damage variable (D) and integrity parameter (a) with norm of inelastic displacement in Mode-II test.

The Mode-I test is performed on a rectangular sample with dimensions of 100  50 mm as shown in Fig. 8. Fig. 8(a) and (b) are respectively for the uniaxial tension and uniaxial compression tests. An interface is imbedded in the middle of the samples. In

vy

vy

a

b

10 mm

vx

Interface

10 mm

vx

10 mm

Coheison, MPa

Numerical implementation 6

vx

Interface

10 mm

Fig. 18. The numerical models and the boundary conditions used in mixed-mode tests (a) tensile-shear loading and (b) compressive-shear loading.

153

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

7 Initial failure surface

Shear stress, MPa

6

Evolved failure surface Numerical implementation

5

Model equation

4 3 2 1 0

0

0.5

1

1.5

2

2.5

3

Tensile stress, MPa Fig. 19. The loading path of tensile-shearing simulation (the evolved failure surface is one of the intermediate surfaces between initial and final failure surfaces).

4.3. Mixed-mode test

10

8 Shear stress, (MPa)

to avoid the contribution of normal stress on the fracture to shear strength, the friction angle adopted for the shearing test is zero. The modelling of shear crack failure is based on the mono-linear softening model, as shown in Fig. 4(b). Figs. 15–17 present the simulation results using the same material properties given in Table 1. Fig. 15 shows the shear stress with shear displacement. The shear stress initially increases to the initial cohesion value and then decreases linearly. The cohesion decreases with the norm of the inelastic displacement linearly and reduces to zero when the critical norm of the inelastic displacement (i.e. wc Þ is reached, as shown in Fig. 16. The integrity parameter a in the shearing test is constant because the normal stress on the interface is compressive due to the applied load on the top boundary; thus, the microdamage variable D cannot be activated for the stiffness (see Fig. 17). Accordingly, both normal and shear stiffnesses are kept unchanged and equal to their initial values. All the results are in good agreement with the analytical results, demonstrating the validity of the numerical implementation of the model.

Initial failure surface Evolved failure surface Numerical

6

4

2

0 -5

-4

-3 -2 Compressive stress, σn (MPa)

-1

0

Fig. 20. The loading path of compressive-shearing simulation.

For the mixed-mode test, the Transjurane sandstone was used again. Two cases were studied including tensile-shearing and compressive-shearing loading. The model’s geometry and boundary conditions are shown in Fig. 18. The dimension of the model is 10  10 mm. A constant velocity (i.e. v x in Fig. 18) is applied at the two side boundaries of the top part of the model. Another constant velocity (i.e. v y in Fig. 18) is exerted on the top boundary. The parameters are listed in Table 1. During the simulation, the average stresses including normal and shear stress along the interface are tracked and the results are plotted in Figs. 19 and 20 for tensile and compressive shearing, respectively, with the consideration of the failure surface change. Good agreement between simulation and analytical results from model equation calculation is achieved, suggesting that the proposed model and its numerical implementation in UDEC can be applied to simulate laboratory tests of soft rocks.

4.2. Mode-II test 4.4. Uniaxial compression test The Mode-II shearing test was also conducted using Transjurane sandstone. The model geometry and boundary conditions are shown in Fig. 14. The dimension of the model is 50  50 mm. A constant velocity (i.e. v x in Fig. 14) is applied at the two side boundaries of the top part of the model. In order to keep the fracture interface in contact during shearing, a constant vertical load (i.e. F y in Fig. 14) is applied on the top of the model. Meanwhile,

The laboratory testing of Transjurane sandstone subjected to uniaxial compression [19] was simulated to verify the proposed cohesive fracture model. The experimental specimen had dimensions of 80 mm in diameter and 160 mm in height, respectively. Accordingly, a numerical model which has same dimensions as the experimental model was built, as shown in Fig. 21(a). To

vy

a

b

c

d

e

45 ϕf 2

or

ϕd 2

A B

Fig. 21. Uniaxial test (a) model discretisation and boundary conditions, (b) final fracture pattern from experiment, (c) final fracture pattern from numerical simulation using the current model, (d) final fracture pattern from numerical simualtion of Kazerani’s model [16] and (e) Coulomb shear failure band.

154

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

45

b

Axial stress, MPa

40 35

c

30

a

25

Experiment

d

Numerical (current model)

20 15

e

Numerical (Kazerani)

10 5 0

0

0.2

0.4

0.6

0.8

1

Axial strain, % Fig. 22. The comparison of experimental and simulation results of the strain–stress relation in uniaxial compression tests of Transjurane sandstone.

capture the randomness of fracture patterns in real uniaxial compression tests, the model was discretised using Voronoi tessellation [18], which can create randomly-sized polygonal blocks (particles) and remove mesh bias [20]. The shared boundaries of each pair of neighbouring blocks are treated as the cohesive fracture interface. The fracture interface is assigned with the cohesive fracture model, as descibed in Section 3. ‘‘Fracturing’’ can only occur along the contact of block when the fracture strength between Voronoi blocks is exceeded. The geometry and boundary conditions of the numerical model are shown in Fig. 21(a). The size

of the particle for the uniaxial compression test is 5 mm. A vertical loading is simulated by applying vertical downward velocity on the top, while the bottom is held fixed using rollers. According to the specimen’s geometry and loading conditions, a plane-strain analysis is used for the modelling [19]. Fig. 21(b) and (c) show the comparison between the experiment and simulation for the final failure pattern using the proposed cohesive fracture model. Fig. 21(d) presents the fracture pattern simulated by Kazerani’s simulation (i.e. from [19]). It can be seen that the current numerical model captures well the fracture pattern observed in the experiment. Comparing the fracture pattern of the two numerical simulations, a closer to the experimental pattern is obtained by the current cohesive fracture model. Fig. 22 shows the comparison of the laboratory test and the numerical simulation for axial strain–stress curves. The laborotory result has shifted left a little by truncating the initial curve before the peak to remove the initial effect of compaction. The simulated peak load is 39.7 MPa and the axial strain at the peak load is 0.31%, which are consistent with the results of the laborotory test. Comparing the two simulation result, the current model has better prediction ability for the peak stress. In terms of post peak behaviour, Kazerani’s model has better capacity to capture the post peak ‘‘brittle’’ behaviour. Nevertheless, the post peak behaviour from experiment may not be correctly obtained due to the possible relatively low stiffness of testing machine. According to the Coulomb shear band theory, which is based on the critical shear stress surface and correlates the shear band orientation with friction angle, the maximum shear stress occurs

Fig. 23. Transjurane sandstone fracture progression under unaxial compression testing. (a) to (e) correspond to the positions in Fig. 22.

155

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

vy

a

b

c

Fig. 24. Brazilian disc test (a) model discretization and boundary conditions, (b) final fracture pattern from experiment and (c) final fracture pattern from numerical simulation.

along a direction of 45° from the major principle stress direction (i.e. A in Fig. 21e), but the actual shear failure band forms along a direction (hm ¼ 45  u=2Þ from the direction of the maximum load (i.e. B in Fig. 21e). However, the Roscoe shear band theory [21], which is based on the critical kinematic surface and relates shear band orientation with dilation angle, proposes that the shear bands must be parallel to the no-extension lines. The no-extension line is necessary for a continuous displacement field [22]. The orientation of the no-extension line is (hr ¼ 45  ud =2Þ with the major principle stress. Due to the fact that the dilatancy angle is usually much smaller than the friction angle, the difference between the Coulomb and Roscoe shear band theories is significant. Authur et al. [23,24] proposed an intermediate shear band theory by considering both the Coulomb and Roscoe shear band theories and experimental evidence. In this theory, the orientation of the shear band is (ha ¼ 45  ðud þ uÞ=4Þ from the major principle stress orientation. For the Transjurane sandstone, the shear failure band and the maximum load directions should form angles of 24.5°, 40° and 32.25° based on the Coulomb and Roscoe shear band theories, respectively. Laboratory testing demonstrated that the Coulomb shear band theory is more applicable for the prediction of shear band orientation for the rock and this may explain the reason for the wide application of the Coulomb strength criterion in rock mechanics. Fig. 23 gives the modelling prediction of the initiation and progression of the fractures in numerical simulation. The stress– strains of (a) to (e) are marked in Fig. 22. The fracture initiates in the pre-peak stage and progresses with compression in random locations. In particular, no crack is observed in the simulated sample until point (a). From that point on (i.e. (a) to (b)), the fracture keeps initiating and propagating. The axial stress keeps ascending until a continuous through-running failure band appears, i.e. the failure band at the left bottom of the model (see Fig. 21b). This stage is normally named the failure stage for a typical uniaxial compression test curve in rock and rock experiences accelerating expansion in its volume in the stage. It is known that during the deformation of rock, the rock testing machine accumulates strain energy due to its limited stiffness. After point (b), the rock sample goes into the post-peak failure stage. In this stage, the rock failure cannot control the release of the strain energy stored in the machine, resulting in the dramatic destruction of the rock sample. This is the reason for the drastic decrease of compressive stress in the laboratory test, as shown in Fig. 22. Two major failure bands finally form in the simulation and this agrees well with the experimental result shown in Fig. 21(b).

4.5. Brazilian disc test The simulation for the Brazilian disc test is based on the authors’ experimental result on Gosford sandstone which has not simulated before. The purpose of the Brazilian disc test is to apply the mixed-mode cohesive fracture model in the soft rock and demonstrate the applicability of the cohesive fracture model. The implementation of other cohesive fracture model can be time-consuming. Thus, only the simulation result of the current cohesive fracture model is presented. Gosford sandstone is a medium-grained sandstone with medium particles, and it is a poorly cemented, immature sandstone with minor quartz overgrowths and contains 20 to 30% feldspar and clay minerals and iron staining [25,26]. The intact sandstone samples were cut into discs with 50 mm in diameter and 25 mm in height. Prior to testing, all sandstone specimens were Table 2 Parameters used in Brazilian disc simulation.

q E

t rt0 c0

2600 kg/m3 4.0 GPa 0.25 4.0 MPa 10 MPa

kn0 ks0

u wr wc

6.0 GPa/mm 3.0 GPa/mm 40° 0.01 mm 0.15 mm

Fig. 25. Comparison of experimental and simulation results of the strain–stress relation in Brazilian disc test.

156

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157

a

b

d

e

c

Fig. 26. Gosford sandstone fracture progression under Brazilian disc testing (a) to (e) are correspond to the positions in Fig. 25.

conditioned by placing them for a week in a climate chamber with a constant relative humidity and temperature. The geometry and boundary conditions are shown in Fig. 24(a). A plane-stress loading analysis was adopted for the modelling. The simulation results gradually converged when decreasing the particle size. Balancing the computation time and result resolution, the particle size used in the simulation was 1 mm. A vertical loading was simulated by applying a vertical downward velocity (i.e. v y Þ on the top while the bottom was held fixed. The detailed parameters used in the simulation are presented in Table 2. The simulated fracture pattern is shown in Fig. 24(c), as compared to the experimental final pattern shown in Fig. 24(b). The simulated load–displacement curve is shown in Fig. 25. The laborotory result has shifted left a little by truncating the initial curve before the peak to remove the initial effect of compaction. Fig. 26 gives the modelling prediction of the initiation and propagation of the fractures in numerical simulation. The stress–strains of (a) to (e) are marked in Fig. 25. The fracture initiates in the pre-peak stage and progresses with compression mostly at locations close to the top and bottom boundaries. In particular, no crack is observed in the simulated sample until point (a). At point (a) only a few randomly positioned non-continuous fractures appear. From that point on (i.e. (a) to (c)), fractures keep initiating and propagating. The compressive load keeps increasing until a continious through-running failure band appears in the middle of the sample (i.e. Fig. 26c) and the peak compressive load is reached. After that, the load decreases with compressive displacement. From the simulation, it is obvious that fracture initiation does not mean rock specimen collapse. Rock collapse only occurs when a continuous failure band appears. The failure band is not straight due to the discrete nature of the model used and the real rock specimen. 5. Conclusion In this paper, a cohesive fracture model considering compression, tension and shear is investigated. The fracture model uses

an effective inelastic displacement to relate tensile strength, cohesion, and stiffness. A micro-damage variable is used to degrade both normal and shear stiffness under tension. The cohesive fracture model is successfully implemented in coupled continuum– discrete element code (i.e. UDEC) using FISH and verified using Mode-I, Mode-II and Mixed-mode tests. The cohesive fracture model is applied to uniaxial compression and Brazilian disc tests of soft rocks, which involve complex inter-block interactions. Good agreements between numerical simulations and experimental results have been achieved. Although linear elastic behaviour was assumed within the material matrix, it is possible to incorporate elasto-plastic models to simulate the material matrix. Acknowledgement The funding received by the second and third authors through ARC Linkage (LP130100884) and by the third author through ARC Discovery Project (DP110104808) is gratefully acknowledged. The support of funding from the Monash Engineering Seed Funding Scheme 2014 is gratefully acknowledged. The fund received by the first author through open fund SKLGDUEK1401 of State Key Laboratory of Deep GeoMechanics and Underground Engineering is also highly acknowledged. References [1] Elices M, Guinea GV, Gomez J, Planas J. The cohesive zone model: advantages, limitations and challenges. Eng Fract Mech 2002;69:137–63. [2] Ruiz G, Pandolfi A, Ortiz M. Three-dimensional cohesive modelling of dynamic mixed-mode fracture. Int J Numer Anal Meth Geomech 2001;52:97–120. [3] Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–82. [4] Carpinteri A, Cornetti P, Valente BP. Cohesive crack model description of ductile to brittle size-scale transition: dimensional analysis vs. renormalization group theory. Eng Fract Mech 2003;70:1809–39. [5] Kazerani T, Yang ZY, Zhao J. A discrete element model for predicting shear strength and degradation of rock joints by using compressive and tensile test data. Rock Mech Rock Eng 2012;45:695–709.

Y. Gui et al. / Computers and Geotechnics 66 (2015) 142–157 [6] Amarasiri AL, Kodikara JK. Numerical modelling of desiccation cracking using the cohesive crack method. Int J Geomech 2013;13(3):213–21. [7] Amarasiri A, Kodikara JK. Determination of cohesive properties for mode I fracture from beams of soft rock. Int J Rock Mech Min Sci 2011;48:336–40. [8] Guiamatsia I, Nguyen GD. A generic approach to constitutive modelling of composite delamination under mixed-mode loading conditions. Compos Sci Technol 2012;72:269–77. [9] Vermeer PA, Borst R. Non-associated plasticity for soils, concrete and rock. Heron 1984;29(3):1–64. [10] Edelbro C. Numerical modelling of observed fallouts in hard rock masses using an instantaneous cohesion-softening friction-hardening model. Tunn Uudergr Sp Tech 2009;24:398–409. [11] Kolluru SV, O’Neil EF, Popovics JS, Shah SP. Crack propagation in flexural fatigue of concrete. J Eng Mech 2000;126(9):891–8. [12] Ferreira LET, Bittencourt TN, Sousa JLAO, Gettu R. R-curve behavior in notched beam tests of rocks. Eng Fract Mech 2002;69:1845–52. [13] Shannon BM. Fracture propagation of cohesive soils under tensile loading and desiccation. PhD thesis, Monash University, Australia; 2013. [14] Cundall PA, Hart RD. Numerical modelling of discontinua. Eng Comput 1992;9(2):101–13. [15] Itasca Consulting Group, Inc. UDEC manual, USA; 2008. [16] Marti J, Cundall P. Mixed discretization procedure for accurate modelling of plastic collapse. Int J Num Anal Mech Gomech 1982;6:129–39. [17] Galvez JC, Cervenka J, Cendon DA, Saouma V. A discrete crack approach to normal/shear cracking of concrete. Cem Concr Res 2002;32:1567–85.

157

[18] Kazerani T, Zhao J. Micromechanical parameters in bonded particle method for modelling of brittle material failure. Int J Numer Anal Meth Geomech 2010;12:1877–95. [19] Kazerani T. Effect of micromechanical parameters of microstructure on compressive and tensile failure process of rock. Int J Rock Mech Min Sci 2013;64:44–55. [20] Wang DS, Du Q. Mesh optimization based on the centroidal voronoi tessellation. Int J Numer Anal Model 2005;2:100–13. [21] Roscoe KH. The influence of strains in soil mechanics. Geotechnique 1970;20:129–70. [22] Bui HH, Fukagawa R. An improved SPH method for saturated soils and its application to investigate the mechanism of embankment failure: case of hydrostatic pore-water pressure. Int J Numer Anal Meth Geomech 2013;37:31–50. [23] Authur JR, Dunstan T, Al-Ane QAJ, Assadi A. Plastic deformation and failure of granular media. Geotechnique 1977;27:53–74. [24] Vermeer PA. The orientation of shear bands in biaxial tests. Geotechnique 1990;40:223–36. [25] Ord A, Vardoulakis I, Kajewski R. Shear band formation in gosford sandstone. Int J Rock Mech Min Sci 1991;28(5):397–409. [26] Gui Y, Zhao G, Khalili N. Experimental and numerical modelling of sandstone bending using modified three-point test. In: Zhao J, Li JC, editers. Rock dynamics and applications – state of the art; May 2013. p. 289–394.