A coupled elastoplastic damage model for brittle rocks and its application in modelling underground excavation

International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

Contents lists available at ScienceDirect

International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms

A coupled elastoplastic damage model for brittle rocks and its application in modelling underground excavation J.C. Zhang a, W.Y. Xu b, H.L. Wang c, R.B. Wang b,n, Q.X. Meng b, S.W. Du d a

College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China Research Institute of Geotechnical Engineering, Hohai University, Nanjing 210098, China c College of Harbour, Coastal and Offshore Engineering, Hohai University, Nanjing 210098, China d Changjiang Institute of Survey, Planning, Design and Research, Wuhan 430010, China b

art ic l e i nf o

a b s t r a c t

Article history: Received 9 February 2015 Received in revised form 29 October 2015 Accepted 30 November 2015

This study develops a coupled elastoplastic damage model for brittle rocks within the framework of irreversible thermodynamics. In the model, degradations in the deformation and strength of rocks induced by damage are considered in the elastic stiffness matrix and plastic yield criterion. The thermodynamic force associated with damage is expressed as a state function of elastic strain and the internal plastic variable. The damage potential and criterion, as functions of the thermodynamic force, are constructed to describe the damage evolution of rocks, and to determine whether they are damaged, respectively. Constitutive relationships, and their numerical formulations of a return mapping algorithm for finite element method implementation, are deduced under three different loading conditions: damage, plastic and coupled plastic damage. For engineering application, the developed model is implemented in the finite element programme Abaqus as a user-defined material model (UMAT). The coupled model is then used for simulating the re-distributed stress fields, displacement and EDZs of the sidewall rock mass around the underground powerhouse of the Wudongde hydropower plant. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Rock damage Coupled elastoplastic damage model Numerical implementation Underground excavation

1. Introduction Brittle rocks (like, granite, shale, marble, and limestone) are always chosen as geological barriers for radioactive waste disposal1,2 and as host rocks for underground caverns.3,4 These rocks, however, are prone to damage during excavation. The integrity of the rocks surrounding these underground openings is compromised and excavation damaged zones (EDZs) are formed (Fig. 1). These zones give rise to considerable degradation of the mechanical, physical, hydraulic, and acoustic properties of the rock.1 EDZs may produce adverse deformation and unstable conditions around the underground caverns, so engineers need to evaluate the EDZ distributions in the surrounding rocks. Knowledge concerning the mechanical behaviours of the brittle rocks would help engineers to understand EDZ developed around an underground opening. Many tests provide this information.5–7 They indicate that the brittle rocks undergo deformation in a number of stages (Fig. 1). These stages are closely associated with four cracking states: initial crack closure, new crack initiation, crack propagation, and crack coalescence, which are represented n

Corresponding author. E-mail address: [email protected] (R.B. Wang).

http://dx.doi.org/10.1016/j.ijrmms.2015.11.011 1365-1609/& 2015 Elsevier Ltd. All rights reserved.

by characteristic stress thresholds: crack closure σcc , crack initiation σci , crack damage σcd , and peak strength σp .2,8 In general, when the stress level exceeds σcd , volumetric strain turns from compression to dilation and the number of acoustic emission (AE) events starts to rise.2,8,9 This accompanies a synchronous development of damage. As the stress level approaches σp , brittle rocks fail and rupture, and this rapidly degrades their mechanical performance.6,10,11 Subsequently, mechanical behaviour of brittle rocks goes into the post-peak regime where bearing capacity stresses drop sharply at first and then decline more slowly and approach the residual strength, σr . Experiments show that in a narrow range from pre- to post-peak, AE events and damage induced by crack propagation and coalescence increase rapidly.6 With this fundamental knowledge about damage properties, the degree and extent of the EDZs around an underground opening can be evaluated.2,5,8 Crack growth is the key factor that drives the development of damage and irreversible deformation in both pre- and post-peak regimes. In order to describe this crack-induced plasticity and damage mechanism, micromechanical models have been developed.10,12,13 In these models, crack opening, closure, friction, and interaction are all taken into account. Thus, micromechanical models can reproduce the physical mechanisms that cause the

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

131

Fig. 1. (a) Stress–strain diagram illustrating the mechanical behaviour of brittle rocks (after Martin et al.,2 Goodman37 and Cai et al.8). (b) Diagram illustrating the distribution of excavation disturbed and damaged zones around a tunnel opening (after Read38).

mechanical behaviour of the brittle rocks.14 However, excessive preoccupation with crack behaviour can result in extremely complex micromechanical models that greatly limit their applicability to real-world rock engineering problems. On the other hand, phenomenological models generated from plasticity and damage mechanics provide an alternative. Although phenomenological models cannot be related to crack-induced plasticity and damage mechanisms, they are easily implemented and applied to real engineering situations.14 In this field of study, which is mainly based on continuum mechanics, researchers introduce damage concepts into their inelastic models to describe the degradation of the mechanical properties of brittle rock-like materials.15,16 These phenomenological models usually consider the effect of damage on stress, strain, elastic stiffness, and strength,17–19 where interactions between damage and plasticity are paid less attention. However, further experimental analyses on brittle rocks6,11 has shown that damage and plasticity are coupled as shown schematically in Fig. 1. Thus, a number of coupled elastoplastic damage models20–22 have been proposed. Shao et al.,22 Chiarelli et al.23 and Mohamad-Hussein et al.24 formulated the basic framework and discussed the mechanics of plastic flow and damage evolution in coupled elastoplastic damage models. Most of the coupled elastoplastic damage models proposed describe the mechanical behaviour of semi-brittle rocks like argillite25,26 and claystone.23 These rocks do not show brittle failure and strength softening, behaviours that are salient features of brittle rocks. In order to describe the mechanical behaviour of the brittle rocks, Chen et al.6 and Zhou et al.11 developed suitable coupled elastoplastic damage models for simulating the plastic hardening effect in the pre-peak regime and rapid damage degradations in elastic stiffness and yield strength in the post-peak regime. These models have achieved considerable success in modelling laboratory tests of brittle rocks. In this study, the authors attempt to advance a coupled elastoplastic damage model to a state where it is possible to apply the model to a practical problem like evaluating EDZs around an underground excavation. This paper presents three main studies. First, a comprehensive analysis is conducted on building a coupled elastoplastic model within a thermodynamics framework. Second, based on return mapping algorithms, the computations are implemented numerically by the finite element method (FEM) in order to apply the coupled elastoplastic model to realistic rock engineering. Following basic principles and formulations, a coupled elastoplastic model based on a modified nonlinear Mohr–Coulomb yield criterion is proposed. Third, as a realistic validation of the model and to demonstrate its application, the proposed model is applied to one of

the large underground excavations for the Wudongde hydropower project, southern China. The numerical modelling results would help rock engineers to design excavation and reinforcement. The work presented in this paper is significant for the application of the coupled elastoplastic damage models in rock engineering.

2. General formulations of the coupled models 2.1. Definition of damage variable Many laboratory tests on granites and other rocks27,28 indicate that stress-induced damage is anisotropic. This changes the structure and affects the mechanical behaviour of the rocks. Thus, a tensor-type damage variable should be introduced into a constitutive model to accurately describe the anisotropic damage.23 For practical engineering application, however, this study uses a scalar isotropic damage variable to build the coupled elastoplastic damage model. This approach may be weaker for modelling the anisotropic degradation and mechanical performance of brittle rocks. Following the tradition of continuum damage mechanics, the degradation ratio of the elastic modulus during stress loading is commonly used to define the scalar isotropic damage variable ω . That is

ω=1−

E , E0

(1)

where E0 is the initial undamaged elastic modulus, and E = (1 − ω) E0 the damaged elastic modulus. From Eq. (1), cyclic tests are useful for measuring the degraded modulus and evaluating the extent of damage to the test materials. 2.2. Thermodynamics framework In an isothermal process of the continuum of damage to rock materials, the Helmholtz free energy is taken as a function of three state variables: elastic strain εe , internal plastic variable γ¯ p and damage variable ω . To simplify its mathematical form, the Helmholtz free energy Ψ (εe, γ¯ p, ω) is generally decomposed into elastic part Ψ e (εe, ω) and plastic part Ψ p (γ¯ p, ω), expressed as follows:

Ψ = Ψ (εe , γ¯ p, ω) = Ψ e (εe , ω) + Ψ p (γ¯ p, ω) .

(2)

The second law of thermodynamics dictates that the inelastic dissipation must be positive in terms of the following Clausius– Duhem inequality:

132

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

σT dε p − dΨ ≥ 0,

(3)

where dε p denotes the incremental plastic strain. Here, σ and ε are the stress and strain vectors, as the simplified versions of corresponding tensors. The superscript “ T ” means the transpose operation of a vector or a matrix. By substituting the differential of the Helmholtz free energy dΨ into the above inequality, we have

σT dε p −

∂Ψ dω ≥ 0, ∂ω

(4)

where dε p is used to denote the increment of plastic strain. The inequality must hold for any state variable results in that the thermodynamic conjugate force for damage Y can be determined by the following partial derivatives of the Helmholtz free energy e

Y=−

p

∂Ψ ∂Ψ − = Y e + Y p, ∂ω ∂ω

(5)

where Y e is the thermodynamic conjugate force for damage caused by elastic strains, which is always taken to be the following expression:

Ye = −

1 e T ∂D e ( ε ) ∂ω ε , 2

(6)

with D being the damaged elastic stiffness matrix as a function of damage variable ω. Y p in Eq. (5) denotes a contribution of the plastic flow to the thermodynamic conjugate force for damage, which is usually regarded as a function of internal plastic variable, with an abstract formulation of Y p = Y p (γ¯ p).6,22

2.3.1. General stress–strain relationship Generally, the assumption of small strain is suitable for the brittle rocks. Thus, the total strain ε can be decomposed into an elastic part εe and a plastic part ε p , as follows:

(7)

For rocks that have been damaged and incur plastic deformation, the stress at each point is determined by the elastic strain εe and the damaged elastic stiffness matrix, as follows: p

e

σ = D (ε − ε ) = (1 − ω) D0 ε ,

(8)

or

εe = C σ =

C0 σ, 1−ω

(9)

where D( D0 ) is the damaged (undamaged) elastic stiffness matrix and its inverse is the damaged (undamaged) elastic compliance C (C 0), that is C = D−1 ( C 0 = D−0 1). For the stiffness matrix, D = (1 − ω) D0 , and for the compliance matrix, C = C 0 / (1 − ω). From Eq. (8), the stress increment can be written in the following differential form:

dσ = D (dε − dε p) − D0 (ε − ε p) dω.

dω = dμ

∂Gω = dμgYω , ∂Y

(10)

Eq. (10) shows the stress increment d σ is composed of two parts: the “positive” part, caused by the incremental elastic strain dεe = dε − dε p on the current damaged elastic stiffness matrix D, and the “negative” part, produced by the incremental damage variable dω on the current elastic strain state εe = ε − ε p . Using this definition, the coupled elastoplastic damage model provides a more universal formulation in which damage and plasticity are just two special cases. They are described below. 2.3.2. Descriptions of damage In terms of describing damage evolution, a damage potential G ω = G ω (Y , ω) is used following a non-normal dissipation law.

(11)

where dμ is the non-negative incremental multiplier for damage. The onset of damage evolution depends on the damage loading condition, which is expressed by

F ω (Y , ω) = 0, dμ ≥ 0 , and F ω (Y , ω) dμ = 0,

(12)

where F ω (Y , ω) = 0 denotes the damage criterion. As the elastic damage loading state is active and no plastic flow is involved, the damage consistency condition can be expressed as follows:

dF ω = fYω dY + fωω dω = 0, where

(13)

ω

f Yω

f ωω

= ∂F /∂Y and

ω

= ∂F /∂ω . By substituting Eq. (5) into

Eq. (13), and noting that dε p = 0 , dμ can be determined by T

( ) dε . dμ = ∂Y

∂εe fωω gYω

(14)

In this condition, the incremental stress–strain relationship of Eq. (10) becomes

dσ = D eωdε,

(15)

where Deω is the tangent elastic damage stiffness matrix, which is defined as follows:

D eω = D −

2.3. Framework of the coupled elastoplastic damage model

ε = εe + ε p , and dε = dεe + dε p.

Thus, the damage variable increment dω can be determined by:

fYω fωω gYω

D0 εe (εe)T D0 .

(16)

2.3.3. Descriptions of plasticity In terms of rock materials, non-associated flow rule is generally used. This means that the plastic strain increment dε p can be determined by a plastic potential G p = G p (σ, γ¯ p, ω) with the following expression:

d ε p = dλ

∂G p = dλ g σp , ∂σ

(17)

where dλ is the non-negative plasticity multiplier. The onset of plastic flow is dependent on the plastic loading condition, which is expressed by

F p (σ, γ¯ p, ω) = 0, dλ ≥ 0, and F p (σ, γ¯ p, ω) dλ = 0. p

(18)

p

where F (σ, γ¯ , ω) = 0 is the plastic yield criterion. As the plasticity loading state is active and no damage evolution is involved, the plastic consistency condition can be expressed as follows: T

dF p = ( f σp) dσ + f γpp dγ¯ p + f ωp dω = 0, ¯

(19)

where f σp = ∂F p/∂σ, f γ¯pp = ∂F p/∂γ¯ p and f ωp = ∂F p/∂ω . By substituting of Eqs. (10) and (17) into Eq. (19), and noting dω = 0, dλ is determined with

( f σp)T Ddε

dλ =

( f σp)T Dgσp − f γ¯ p

∂γ¯ p ∂ε p

T

( )

. g σp

(20)

Thus, in this condition, the incremental stress–strain relationship of Eq. (10) is of elastoplasticity and simplified as follows:

dσ = D epdε,

(21)

where Dep is the tangent elastoplastic stiffness matrix, which is defined as follows:

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

133

T

Dg σp ( f σp) D

D ep = D −

(f σ)T Dg σp − f γ¯ p

T

∂γ¯ p

( ) ∂ε p

. g σp

(22)

2.3.4. Descriptions of coupled elastoplastic damage As damage and plasticity loading states are usually active, damage evolution and plastic flow occur in a coupled way. For this situation, dε p and dω are both commonly determined by solving a set of linear equations composed of damage and plastic consistency conditions, Eqs. (13) and (19), which is

⎡ A11 A12 ⎤ ⎧ dμ ⎫ ⎧ B1⎫ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ , ⎣ A21 A22 ⎦ ⎩ dλ ⎭ ⎩ B2 ⎭

(23)

with dμ and dλ as the unknown variables. The solution of Eq. (23) is

⎡ A22 −A21⎤ ⎧ B1⎫ ⎧ dμ ⎫ 1 ⎨ ⎬ = ⎢ ⎥⎨ ⎬ . A11A22 − A12 A21 ⎣ −A12 A11 ⎦ ⎩ B2 ⎭ ⎩ dλ ⎭

(24)

In Eqs. (23) and (24), A11, A12, A21, A22, B1 and B2 are defined by T

A12 = f Y εe D 0 g σp + f Y

ω

ω

A11 = f ω gYω ,

ω

( )

( )

T

p ⎞T

⎛ ∂¯γ p T A22 = f σ Dg σp − fγ¯ p ⎜ p ⎟ g σp, ⎝ ∂ε ⎠

⎛ p T p⎞ A21 = ⎜ f σ D 0 εe − fω ⎟ gYω , ⎝ ⎠

B1 = fYω ( εe) D0 dε,

T ∂Y ⎛ ∂¯γ p ⎞ p ⎜ ⎟ g ∂¯γ p ⎝ ∂ε p ⎠ σ

( )

(25)

(27)

where Depω is the coupled elastoplastic damage stiffness matrix, which is defined as follows:

=D−

T

T

Dg σp A11 ( f σp) D − A12 ( εe) D0 A11A22 − A12 A21

(29)

(26)

dσ = D epωdε,

(

calculation for the fully elastic trial stress σtrial and the trial thermodynamic conjugate force Y trial using purely elastic relationships, which are

σtrial = σn + DΔε,

T

B2 = ( f σp) Ddε.

In Eqs. (23) and (24), the element components A12 and A21 indicate the coupled influences of plasticity and damage. In the coupled loading state, the incremental stress–strain relationship of Eq. (10) becomes

D epω

Fig. 2. Diagram illustrating the principle of return mapping for stress updating (after 30, 31).

).

(28)

3. Numerical implementation of the coupled model The return mapping algorithms proposed by Simo and Hughes29 are now in general use in FEM. Generally, in a return mapping algorithm, two calculation steps are usually involved, which are elastic prediction and plastic correction,30,31 as schematically illustrated in Fig. 2. However, the algorithms currently used are mainly used for updating stress and plastic strain in a solitary physical process. On the basic principles, the formulations presented in the following section provides a computation procedure, which not only updates stress and plastic strain, but also updates the thermodynamic conjugate force and damage variable in the coupled plasticity and damage processes. 3.1. Elastic prediction Numerical implementation of the coupled elastoplastic damage model is schematically illustrated with Fig. 2. It is assumed that the main finite element programme returns the state variables: σn , ε en , ε np , Yn , ωn and Dn at time tn , and then provides the total strain increment Δεn + 1 at the integrating points in each element at time tn + 1. The constitutive integrating formulations perform the

Y trial = −

1 e ( ε n + Δε)T D0 ( ε en + Δε) + Ynp . 2

(30)

Because the following calculation process is carried out at tn þ 1, the subscript “ n + 1” is omitted. According to the trial results, F p = F p (σtrial , γ¯np , ωn ) and F ω = F ω (Y trial, ωn ) may be greater or smaller than zero. These two functions exhibit three types of combinations, corresponding to three inelastic loading states, with their related numerical formulations deduced as follows.

3.2. Stress and damage update 3.2.1. Coupled plastic damage loading state If the trial results make F ω > error and F p > error (where error is a very tiny value set for numerical computation), this point is in the coupled plastic damage loading state, with both plastic flow and damage evolution. Therefore, a part of the total strain increment Δε will translate into plastic strain to reduce the trial stress, and during this process the elastic part Y e of the thermodynamic conjugate force for damage will also reduce and translate into the plastic part Y p . At the same time, the damage variable increases. In the coupled plastic damage loading state, from time tn to tn þ 1 the stress and the thermodynamic conjugate force for the damage update as

σ = σn + Δσ = σn + DΔεe − D0 ε en Δω, Y = Yn + ΔY = −

(31)

T 1 e ( εn + Δεe) D0 ( εen + Δεe) + Y p γ¯np + Δγ¯ p . 2

(

)

(32)

As illustrated in Fig. 2, the plastic damage coupling corrector stress Δσpω from σtrial to σ can be written as

Δσpω = σ − σtrial = − ΔλDg σp − ΔωD0 ε en ,

(33)

and the plastic damage coupling corrector thermodynamic conjugate force ΔYpω from Y trial to Y is expressed as

134

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

T

ΔYpω = Y − Y trial = ( ε en ) D0 ( Δε − Δεe) + Y p γ¯np + Δγ¯ p

(

)

T 1⎡ ⎤ + ⎢⎣ (Δε)T D0 Δε − ( Δεe) D0 Δεe⎦⎥ − Y p γ¯np , 2

( )

eT

T

(34)

e

where (Δε) D0 Δε and (Δε ) D0 Δε are, in some sense, infinitesimally high-order variables compared to (ε en)T D0 (Δε − Δεe), and can be omitted. Additionally, a differential is used to represent the increment of ΔY p from γ¯np to γ¯np + Δγ¯ p . Then Eq. (34) becomes T

ΔYpω = Y − Y trial = ( ε en ) D0 Δε p +

∂Y p p Δγ¯ . ∂γ¯ p

(35)

However, the real route from the elastic prediction states (σtrial , γ¯np , Y trial, ωn ) to the newly updating states (σ, γ¯ p, Y , ω) is usually finished after several iterations, instead of being directly finished in one step, as described by Eqs. (33) and (35). At every iteration, according to Huang and Griffiths,31 the plastic yield function and damage criterion can be expanded as two Taylor series at the current state (σ(k ) , γ¯ p (k ) , Y (k ) , ω(k ) ), that is

F

ω (k + 1)

=F

ω (k )

+ fYω (k) δY p(kω) + fωω (k) δω (k) ,

(36)

T

(

)

F p (k + 1) = F p (k) + f σp (k) δ σ (pkω) + f γ(kp) δγ¯ p (k) + f ω(k) δω (k) . ¯

(37)

where k indicates iteration time. To simplify Eqs. (36) and (37), the symbols δσ(pkω) , δεp (k ) , δγ¯ p (k ), δY p(kω), and δω(k ) are used to represent the increments of stress, plastic strain, internal plastic variable, thermodynamic conjugate force, and damage variable from k to k + 1, respectively. Here, the symbol ‘ δ ’ is used to represent the small difference between two iterations, as distinct from the total difference ‘ Δ’ between the starting and ending points (as shown in Fig. 2). They are

δ σ (pkω) = σ (k + 1) − σ (k) = − δλ (k) D (k) g σp (k) − δω (k) D0 ε en ,

(38)

δ ε p (k) = ε p (k + 1) − ε p (k) = − δλ (k) g σp (k) ,

(39)

⎛ ∂γ¯ p ⎞(k) T δγ¯ p (k) = γ¯ p (k + 1) − γ¯ p (k) = δλ (k)⎜ p ⎟ g σp (k) , ⎝ ∂ε ⎠

(40)

T ⎤(k) ⎡ ∂Y p ⎛ ∂γ¯ p ⎞ ⎥ T δY p(kω) = Y (k + 1) − Y (k) = δλ (k)⎢ ( ε en ) D0 + g σp (k) , ⎜ ⎟ ⎢⎣ ∂γ¯ p ⎝ ∂ε p ⎠ ⎥⎦

(41)

δω (k) = ω (k + 1) − ω (k) = δμ(k) gYω (k) .

(42)

(k ) , A11

(k ) , A12

(k ) , A21

(43)

(k ) A22

where are the same as those listed in Eq. (25). However, in the numerical implementation, it should be noted that B (1k ) and B (2k ) are different from that ( B1 and B2) listed in Eq. (26). They are

B (1k)

=

F ω (k ) ,

B (2k)

=

F p (k ) ,

(44)

With δμ(k ) and δλ(k ) solved from Eq. (43), the increments: δσ(pkω) ,

δγ¯ p (k ) , δY p(kω), and δω(k ) can be calculated after which the state variables σ(k + 1) , εp (k + 1) , εe (k + 1) , γ¯ p (k + 1) , Y (k + 1) , and ω(k + 1) can be updated.

δεp (k ) ,

3.2.2. Damage loading state If the trial results make F ω > error and F p ≤ error , this point is in the damage loading state. Because no plastic flow occurs, the total strain increment Δε is entirely an elastic one, with no incremental plastic part included, that is Δεe = Δε and Δε p = 0 . The variables updating formulations can be simplified from Eqs. (38), (41) and (42). 3.2.3. Plastic loading state If the trial results make F ω ≤ error and F p > error , this point is in the plastic loading state. In this case, the cutting plane algorithm29 is adopted in this study, as shown in Fig. 2. Because no damage evolution ( Δω = 0) occurs, the stress, plastic and elastic strains updating formulations can be exactly the same as the classical return mapping algorithm. For the three loading states, at the end of each iteration, convergence is evaluated. The values of the damage and plastic functions are calculated. If F p (σ(k + 1) , γ¯ p (k + 1) , ωn ) ≤ error and

F ω (Y (k + 1) , ω(k + 1) ) ≤ error are both satisfied, the iteration ends. 3.3. Programming UMAT of the coupled model in Abaqus Aiming at engineering application, the numerical formulations mentioned above are programmed as a UMAT subroutine in Abaqus/Standard.32 It is worth noting that the strain vector and the material Jacobian matrix in Abaqus are engineering versions, following the traditions of engineering mechanics, not the tensor versions used previously. Table 1 summarises the computational procedures required for the coupled elastoplastic damage model. The table includes the procedures that are not described in Section 3.3. The specific functions used for programming the UMAT are discussed in Section 4.

4. Specific functions used for the numerical study

It is worth noting that F ω (k + 1) and F p (k + 1) are much closer to zero than F ω (k ) and F p (k ) , during iterations. Assuming F ω (k + 1) = 0 and F p (k + 1) = 0,31 an equation system similar to Eq. (23) for solving the damage multiplier increment δμ(k ) and plastic multiplier increment δλ(k ) , can be obtained as follows:

⎡ A ( k ) − A (k ) ⎤ ⎧ (k ) ⎫ ⎧ δμ(k) ⎫ 1 21 ⎥ ⎪ B 1 ⎪ ⎢ 22 ⎨ ⎬ = ⎨ ⎬, (k ) (k ) (k ) (k ) ⎢ (k ) (k ) ⎥ ⎪ (k ) ⎪ ⎩ δλ (k) ⎭ A11 A22 − A12 A21 ⎣ −A12 A11 ⎦⎩ B 2 ⎭

As described above, in addition to the coupled loading states, the coupled elastoplastic damage model has another two inelastic loading states: damage and plasticity. Their numerical computation procedures are a little different because there is no coupled (k ) (k ) component A12 or A21 involved as there is for solving δμ(k ) or δλ(k ) .

In the coupled model framework and with the numerical formulations mentioned above, specific plastic and damage functions are selected to build a coupled elastoplastic damage model, which is used for modelling laboratory tests and simulation of damage around underground excavations. 4.1. Plastic functions The description of rock plasticity requires a combination of a plastic yield criterion, a plastic hardening function, and a plastic potential. In this study, a modified hyperbolic Mohr–Coulomb criterion is adopted to describe the plastic behaviour of brittle rocks. This provides a factor for plastic hardening and damage degradation of strength, which is

Fp =

(qR (θ ))2 + (ξc cos φ)2 + p sin φ − h (1 − ω) c cos φ = 0,

(45)

where c and φ are the cohesion and frictional angle strength parameters of the Mohr–Coulomb materials. The variables p , q , and θ are, respectively, the mean stress, generalised deviatoric stresses, and Lode’s angle. These can be calculated using stress tensors:

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

135

Table 1 Computation procedure of the coupled elastoplastic damage model. 1

Input Vε at the n + 1-th load step; Input the converged values of ε en , ε np , σn , Yn and ωn at the n -th load step.

2

Make elastic predictions: σ trial = σn + (1 − ωn ) D 0Δε and Y trial = −

3

Initialisation: set k = 0 ,

4

IF F ω (Y , ωn ) ≤ error AND F (σ

5

εp (0)

=

ε np , εe (0) trial , γ¯ p , n

=

ε en

+ Vε , Δλ

(0)

= 0,

Δμ (0)

1 e (ε 2 n

+ Δε)T D0 (ε en + Δε) + Ynp

= 0 , σ (0) = σ trial , Y

(0)

=Y

trial

ωn ) ≤ error , THEN pure elastic loading state, GO TO 7. ELSE IF THEN GO TO 5.

5-1. F ω (Y (k ) , ω(k ) ) > error AND

5-2. F ω (Y (k ) , ω(k ) ) ≤ error AND

5-2. F ω (Y (k ) , ω(k ) ) > error AND F (σ(k ) , γ¯ p (k ) , ω(k ) ) > error

F (σ(k ) , γ¯ p (k ) , ω(k ) ) ≤ error [1] Calculate:

F (σ(k ) , γ¯ p (k ) , ω(k ) ) > error [1] Calculate:

[1]

Solve Eq. (43) to calculate: δλ (k ) ; δμ (k ) .

[2]

Calculate: δω(k ) = δμ(k ) gYω (k ) ;

δ σ(pkω)

= − δλ(k ) D(k ) g σp (k ) − δω(k ) D0 ε en ;

δY p(kω)

= δλ(k )[(ε en)T D0 +

[3]

Update variables:

δω(k ) = (F ωgYω /A11)(k ) ; δ σ(ωk ) = − δω (k ) D0 ε en . [2]

Update variables:

δλ [2]

(k )

= (F

σ

=

σ( k )

+

δ σ(ωk ) ;

εp (k + 1) = εp (k ) ;

δ σ(pk ) = − δλ(k ) D(k ) g σp (k ) .

Update variables:

ω(k + 1) = ω(k ) + δω (k ) ; (k + 1)

p/ A )(k ) ; 22

εp (k + 1) = εp (k ) − Cδ σ(pk ); σ (k + 1) = σ( k ) +

σ (k + 1) = σ( k ) + δ σ(pkω) ;

δ σ(pk ) ,

Y (k + 1) = Y (k ) ; k = k + 1.

Y (k + 1)

=

⎡ ⎢ ⎢Y ⎢⎣

+

∂Y p ⎛ ∂γ¯ p ⎞ δλ p ⎜⎜ p ⎟⎟ ∂γ¯ ⎝ ∂ε ⎠

T

εp (k + 1) = εp (k ) − Cδ σ(pkω) ;

ω(k + 1) = ω(k ) + δω (k ) ;

γ¯ p (k + 1) = [γ¯ p + δλ (∂γ¯ p/∂ε p)T g σp](k ) ;

εe (k + 1) = εe (k ) ;

∂Y p ∂γ¯ p T (k ) p (k ) ( ) ] gσ ∂γ¯ p ∂ε p

εe (k + 1) = εe (k ) + Cδ σ(pkω) ; γ¯ p (k + 1) = [γ¯ p + δλ (∂γ¯ p/∂ε p)T g σp](k ) ;

⎤ (k ) ⎥ g σp⎥ ; ⎥⎦

Y (k + 1) = Y (k ) + δλ(k )[(ε en)T D0 +

∂Y p ∂γ¯ p T (k ) p (k ) ( ) ] gσ ; ∂γ¯ p ∂ε p

k = k + 1.

k = k + 1. 6

IF F ω (Y (k ) , ω(k ) ) ≤ error AND F (σ(k ) , γ¯ p (k ) , ω(k ) ) ≤ error , Iteration is converged. THEN GO TO 7. ELSE IF, THEN GO TO 5.

7

IF variables come from 4, THEN σ = σ trial , ε p = ε np , εe = ε en + Δε and the tangent stiffness matrix D = (1 − ω(k ) ) D0 is returned to the FEM main programme. IF variables come from 5, THEN σ = σ (k ) , ε p = ε np (k ) , εe = ε en + Vε − εp (k ) , the tangent stiffness matrices: Deω (k ) in Eq. (16) for 5-1; Dep (k ) in Eq. (22) for 5-2; and Depω (k ) in Eq. (28) for 5-3 are returned to the FEM main programme.

p=

1 σii , q = 3

1.5Sij Sij and θ =

1 −1⎛ 27J3 ⎞ sin ⎜ − 3 ⎟ 3 ⎝ 2q ⎠

(46)

where Sij denotes the deviatoric stress tensor ( Sij = σij − pδ ij ), and δ ij denotes the Kronecker tensor. J3 is the third invariant of the deviatoric stress tensor, which is the determinant of Sij , namely J3 ¼det Sij . In Eq. (43), R (θ ) is a function representing the dependency of rock strength on Lode's angle; this function determines the shape of the modified Mohr–Coulomb criterion in the deviatoric plane,33 as follows:

R (θ ) =

3 − sin φ , K= 2 2 ∘ 3 + sin φ K + (1 − K ) sin (45 − 1.5θ ) K

(47)

where K represents the ratio of the triaxial extension (at θ = − 30°) to the triaxial compression (at θ = 30°). In Eq. (45), ξ is the eccentricity ratio of the hyperbolic Mohr– Coulomb yield surface in the p– q meridian plane. The eccentricity ratio rounds the apex of the cone and reduces the tensile strength σt to a reasonable value. As the prediction of Griffith’s criterion and empirical observations in rock engineering show,2,8 the tensile strength of rocks and rock masses σt is about 1/8 to 1/20 times compressive strength σc . On average, σt ¼ 1/15 σc . Thus, it is reasonable to take ξ ¼ 0.8 in this study as a rough estimation of σt . In Eq. (45), h denotes the plastic hardening function, which describes that evolution laws from its initial yielding of the subsequent yield surface. The multiplication of h and (1  ω) makes the modelling of pre-peak plastic hardening and post-peak damage softening deformation possible. Inspired by Chen et al.,6 h is expressed as follows:

h = h0 + (h1 − h0 )

γ¯ p , b + γ¯ p

(48)

where the parameters h0 and h1 represent the initial and peak yielding thresholds. b is a parameter to describe plastic hardening rate. For building a concrete model, the internal plastic variable γ¯ p

is traditionally defined as the equivalent plastic shear strain,22,34 that is

γ¯ p =

∫ dγ¯ p = ∫

2 (de p)T de p/3 ,

(49)

p

where de is the deviatoric plastic strain and is defined by de p = dε p − δ i dεvp /3 ( δ i = 1, as i=1, 2, 3 and δ i = 0, as i=4, 5, 6). Here, dεvp is the incremental plastic volumetric strain. To describe the plastic flow, a plastic potential G p with a mathematic expression similar to the yield function is used, that is

Gp =

(qR (θ ))2 + (ξc cos ψ )2 + p sin ψ ,

(50)

where ψ is the dilation angle, with its formulations proposed by Vermeer et al.35 as follows:

sin ψ = where

dεvp , −2dε1p + dεvp

dε1p

ψ = arc sin

dεvp , −2ε1p + εvp

(51)

is the increment of the major principal plastic strain.

4.2. Damage functions Following the work of Chen et al.,6 a similar function is used for the damage criterion, and is very suitable for describing the evolution of the damage variable, that is

F ω (Y , ω) = ωc − ωc exp (αY ) − ω = 0,

(52)

where ωc represents the asymptotic damage value of ω . α is a parameter that describes the evolution rate of the damage variable. For simplicity, the damage evolution is regarded as being associative in this study. This means that the damage potential is the same as the damage criterion, namely G ω = F ω . For simplicity, the plastic component Y p of the thermodynamic conjugate force for damage Y used in this study is specifically defined as a linear function of the equivalent plastic shear strain γ¯ p . That is

136

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

Y p = − κγ¯ p

(53)

where κ is a material parameter for the stress with dimensions of [F/L2].

5. Application of the model: modelling an underground excavation 5.1. Wudongde hydroelectric project overview The Wudongde reservoir, the fourth largest hydropower plant in China, has being built on the downstream of the Jinsha River. The total installed capacity of the power plant will be 10,200 MW to be generated by 12 hydroelectric turbines evenly spaced in two underground cavern complexes being excavated in both abutments of the dam. There will be three main underground caverns in each complex: the powerhouse, designed to be 333 m long, 32.5 m wide, and 89.3 m high, the transformer house, 270  18  35 m3 high, and the tailrace surge chamber, 215.7  40  111.5 m3 high. The underground cavern complexes are being excavated in the third Section (Pt2l3) of the Mesoproterozoic Luoxue Formation (Pt2l). In this area, the strata form a monocline striking generally 260–280° (essentially EW) and dipping steeply at 60–80°. At the Wudongde reservoir, the Luoxue Formation is mainly thick-bedded limestone (Fig. 3). The rock is quite competent. 5.2. Triaxial tests and modelling of the Wudongde limestone In order to understand the mechanical behaviour of the limestone, triaxial compression tests were conducted on three cylinder specimens of Wudongde limestone. The cores were obtained from an exploration adit near the powerhouse excavation site. The confining pressures applied to the specimens were 4, 8, and 12 MPa, which were estimated from the in situ stresses. A servocontrolled triaxial rheological test system36 was used to conduct the triaxial compression tests. The displacement stroke control mode was used to apply deviatoric stress to the ends of the specimen until the specimen fractured. The displacement loading rate was 0.02 mm/min (about 2  10  4 ε /min). Fig. 4 shows the triaxial test results. Note that the limestones show strength enhancement in the pre-peak deformation stage when the deviatoric stresses exceed the crack damage stress threshold σcd . This is followed by a sharp ‘stress drop’ in the postpeak deformation stage as the samples are pushed beyond their

peak strength. This ‘stress-drop’ is closely related to the brittle failures. The three photos in Fig. 4 show the brittle failure modes of the Wudongde limestone samples. For validating the coupled elastoplastic damage model, the UMAT for this material was adopted to simulate the limestone triaxial compression tests within Abaqus. The mechanical parameters of the limestone for the developed model are listed in Table 2. They were determined as follows: (a) Elastic parameters: the initial undamaged elastic modulus E0 and Poisson's ratio ν were determined from the test data in the linear elastic deformation stage. (b) Plastic parameters: Taking the test data for peak strengths under different confining pressures, Mohr–Coulomb criterion is used to determine cohesion c and frictional angle φ . Using the method suggested by Vermeer et al.,35 dilation angle can be derived with Eq. (51). (c) Hardening parameters: h0 is set as h0 = σcd /σp (Cai et al.8 studied this ratio). The values of h1 and b are calibrated from the relationship between the subsequent yield stresses and their corresponding equivalent shear strain γ¯ p , from σcd to σp . (d) Damage parameters: the observed residual strengths σr are approximately 0.38–0.66 times the peak strengths σp . This means that, very roughly, there is a reduction of about 62–34% from peak strengths to residual strengths. According to Eq. (52), it is reasonable to set ωc equal to 0.6. α and κ are difficult to calibrate directly from the test data, so in order to obtain a better simulation, they are estimated by repeated calculation. With the mechanical parameters given in Table 2, the Abaqus software can simulate the elastoplastic deformations and damage behaviour of the limestone. The numerical modelling results are also presented in Fig. 4. It is observed that the modelling stress– strain curves show a good agreement with the test data. Additionally, the modelling damage curves reveal the damage characteristics of the brittle limestone. In the unstable cracking growth stage and post-peak stage, limestone damage increases sharply. The shape of the damage curves are similar to the curves of cumulative AE counts reported by Chen et al.6 5.3. Excavation and EDZs at the Wudongde powerhouse site Fig. 5a shows a section shape of the underground powerhouse at chainage 0 þ100 in the right abutment. The section is sliced through the alignment of the hydroelectric turbine No. 12. The

Fig. 3. Photo of the thick-bedded limestone in the Luoxue Formation.

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

137

Fig. 4. Three plots comparing the test results and the numerical modelling results for three specimens of Wudongde limestone subjected to triaxial compression. The photographs show the brittle failure modes of the three limestone specimens after at the confining pressures of 4, 8, and 12 MPa. Table 2 Mechanical parameters used for modelling the triaxial compression tests of the Wudongde limestone. Elasticity

Plasticity

Hardening

Damage

E0 (MPa)

ν

c (MPa)

φ (deg)

ψ (deg)

h0

h1

b

ωc

α (MPa-1)

κ (MPa)

28,530 ( σ 3 ¼ 4 MPa) 35,530 ( σ 3 ¼ 8, 12 MPa)

0.3

17.85

54.3

54.4

0.77

1.21

0.00018

0.60

0.095

750

Fig. 5. (a) Section of the underground powerhouse, along the alignment of the hydroelectric turbine No. 12, showing the excavation sequences and the boreholes used for the acoustic (P-wave) wave velocity tests. (b) Results of P-wave velocity tests from the four boreholes.

138

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

excavation of the underground powerhouse has been designed to be done in a series of 10 sequences. Generally, P-waves are used to evaluate EDZs4 by comparing the relative changes in the acoustic wave velocities in damaged rocks to wave velocities in the undamaged rock. The following simple equation shows the relationship:

ω=1−

V2 E =1− , E0 V02

(54)

where V and V0 are the P-wave velocity of the damaged and undamaged rocks, respectively. After excavation sequences I and II in Fig. 5a were completed in April 2014, four boreholes: BH1 to BH4 (Fig. 5a) were drilled in the arch. Each borehole was about 9.2 m long. To determine the the EDZs in the surrounding rocks, acoustic P-wave velocities were measured from these boreholes. P-wave test results are shown in Fig. 5b. It was found that the average P-wave velocity in the undamaged zones was about 5142 m/s, whereas the P-wave velocities in the EDZs were in the range of 3774–4878 m/s. The EDZ extended to distances of about 1.2–2.0 m into the surrounding rocks. According to Eq. (54), the maximum damage variable of the surrounding rocks is about ω = 1 − V2/V02 ¼1  37742/ 51422 ¼0.46. 5.4. Numerical modelling of excavation 5.4.1. Numerical model For modelling the underground powerhouse excavation, a finite element model with the dimensions x  y  z ¼444  153  650 m was constructed in the right abutment. The x -axis is roughly perpendicular to the river channel oriented east–west and the z -axis is vertical. For simplicity, the model only encompasses the powerhouse and does not take other underground chambers and water conveyance tunnels into account. Hexahedral meshing (that is, 8-noded elements) were used to construct the model. The elements

near the opening surface are about 2.0–3.0 m in size. The numerical model is shown in Fig. 6. There are 72,830 model elements and 78,439 nodes in the model. 5.4.2. Boundary and in situ stress conditions The model has five planar boundaries and one irregular boundary, the irregular boundary being the ground surface (the mountain under which the chambers are being excavated). Boundary conditions are the five planar boundaries are fixed and the mountain surface boundary is free. The top of the powerhouse chamber ranges from about 140 to 470 m below the ground surface and the side walls near the river valley are about 100–120 m inside the canyon wall. Field tests indicate that the in situ stresses caused by vertical gravity are dominant. Thus, the initial in situ stress field applied to the model is simply gravitational along the z -axis. The two lateral pressure coefficients are about 0.45 along the x -axis and 0.65 along the y -axis. 5.4.3. Mechanical parameters The surrounding rocks are of rock masses, in which many macro joints and fracture are contained. This is rather different from the intact rock specimens tested in the laboratory that had no obvious cracks. Thus, most of the mechanical parameters of the surrounding rocks have not been measured directly but have instead been estimated or back-analysed based on the results of in situ monitoring. For this study, the engineers suggested the density, elastic and plastic parameters (Table 3) of the Wudongde limestone rock masses in the underground excavation site, according to the in situ geological investigations. The plastic hardening and damage parameters were estimated, taking the following points into consideration: (a) According to the work presented by Cai et al.,8 for massive to moderately jointed rock masses, the hardening parameter

Fig. 6. Numerical model: (a) a view of the whole model; (b) elements of the powerhouse showing the excavation sequences (as viewed from upstream side).

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

139

Table 3 Mechanical parameters of the surrounding rock mass used for modelling underground excavation. Density

Elasticity

Plasticity

Hardening

ρ (kg/m3)

E0 (MPa)

ν

c (MPa)

φ (deg)

ψ (deg)

h0

h1

b

ωc

α (MPa  1)

κ (MPa)

2600

15,000

0.3

2.0

43

33

0.80

1.20

0.0003

0.6

0.3

400

h0 = σcd /σcm are in the range of 0.8–0.9 (where σcm denotes the uniaxial compressive strength of the rock masses). In this study, h0 ¼ 0.8 was selected. (b) Compared with the intact rock specimens tested in the laboratory, the limestone surrounding rocks in the construction site are relatively “softer” in elastic stiffness and “weaker” in strength (compare Table 3 and Table 2). This means that the surrounding rocks undergo much more plastic deformations than the intact rock specimens at a same stress level. A bigger b parameter can describe this difference in deformation. b ¼ 0.0003 and h1 ¼1.20 produce relatively smooth stress– strain curves whose peak strengths arrive at the Mohr– Coulomb envelope. (c) It is assumed that the limestone surrounding rocks have similar damage properties like the intact rock specimens. As their deviatoric stress levels reach peak strengths, their damage variables rise to about 0.10 (see Fig. 4). According to this assumption, the asymptotic damage parameters ωc ¼0.6, α ¼0.3, and κ ¼400 can simulate the underground excavation well. The mechanical parameters used for the limestones surrounding the powerhouse chamber are listed in Table 3. 5.4.4. Simulation results and analysis The excavation of the underground powerhouse was numerically simulated with the mechanical parameters given in Table 3.

Damage

The redistributed stress contours after excavation are shown in Fig. 7a and b. The values of the major principal stresses σ1 in the stress concentration areas are about 19–23 MPa. Near the openings in the sidewalls, redistributed minor principal stresses σ 3 reduce to tensile stress states. This indicates that the stress states and deformation of the surrounding rocks move into the postpeak regime. Fig. 7c shows the displacements of the surrounding rocks. The maximum displacements are about 35 mm, which are mainly observed at the tops of the central piers between foundation pits for the turbines. As shown in Fig. 5a, 2 multi-point displacement extensometers: DM_1 and DM_2 were installed in the arch dome (EL 855.0 m) and upstream sidewall (EL 838.3 m), respectively. They were used to measure the deformations of the surrounding rocks. The convergence displacements monitored by DM_1 and DM_2 in November 2014 were about 12.8 mm and 16.5 mm. The modelled displacements (see Fig. 7c) at these two positions are about 22 mm and 18 mm, respectively. Considering that any displacement that took place during excavation could not be observed, the predicted displacements are considered acceptable. Figs. 7d and 8 show approximations of the EDZ distributions in the surrounding rocks. Most of the surrounding rocks do not undergo serious damage, because their damage values are mainly around 0.15–0.20. However, the maximum damage values in the central piers can be greater than 0.5. This indicates that the surrounding rocks in the central piers can be more seriously damaged because they are unloaded from three directions. Overall, the

Fig. 7. Simulation results of (a) major principal stress σ1 (MPa), (b) minor principal stress σ 3 (MPa), (c) Displacements (m) and (d) EDZs in the surrounding limestones. All views are from the upstream side.

140

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

damage variables from the model agree well with the variables estimated by acoustic P-wave velocities.

6. Conclusions A coupled elastoplastic damage model for brittle rocks was developed on the basis of irreversible thermodynamics. Based on non-negative inelastic energy dissipation, the thermodynamic conjugate force for the damage is expressed as the sum of elastic and plastic parts, which are state functions for the elastic strain and internal plastic variable, respectively. This thermodynamic conjugate force impels the development of rock damage. Damage degradation in deformation and strength properties was considered. This resulted in coupled influences between damage and elastoplasticity. Constitutive stress–strain relationships for the coupled model were deduced under damage, plastic and coupled plastic damage loading conditions. Based on the principles of the return mapping algorithm, general constitutive integrating formulations of the coupled elastoplastic damage model carried out in FEM were deduced. For engineering application, a plastic yielding criterion on the basis of the modified nonlinear Mohr–Coulomb model and a damage criterion as a negative exponential function of the thermodynamic conjugate force for damage were introduced to code a specific UMAT on Abaqus. For validating the developed model and its engineering applicability, it is firstly used for modelling the triaxial compression tests of a type of limestone in Wudongde hydropower plant. A comparison between numerical prediction and test data indicates that the developed elastoplastic damage coupled model can describe plastic hardening (pre-peak) and damage softening behaviours (post-peak), respectively. The developed model is then used for modelling the re-distributed stress fields, displacement and EDZs after the excavation of the powerhouse. This work is very useful for providing a good prediction for the engineers before their rock excavations are started.

Acknowledgements This study was financially supported by the National Natural Science Foundation of China (Grant nos. 51479049, 51209075 and 11272113) and the Natural Science Foundation of Jiangsu Province (Grant no. BK20130846).

References

Fig. 8. Sections through the underground powerhouse showing EDZs in the surrounding rocks: (a) Section along the alignment of the foundation pit for hydroelectric turbine No. 12; and (b) Section along the alignment of the central pier between the chambers for turbine No. 11 and turbine No. 12. The left side is the upstream sidewall.

1. Tsang CF, Bernier F, Davies C. Geohydromechanical processes in the Excavation Damaged Zone in crystalline rock, rock salt, and indurated and plastic clays—in the context of radioactive waste disposal. Int J Rock Mech Min Sci. 2005;42:109– 125. 2. Martin CD. Seventeenth Canadian Geotechnical Colloquium: the effect of cohesion loss and stress path on brittle rock strength. Can Geotech J. 1997;34:698– 725. 3. Liu HB, Xiao M, Chen JT. Parametric modeling on spatial effect of excavationdamaged zone of underground cavern. J Cent South Univ. 2013;20:1085–1093. 4. Yan P, Lu WB, Chen M, Hu YG, Zhou CB, Wu XX. Contributions of in-situ stress transient redistribution to blasting excavation damage zone of deep tunnels. Rock Mech Rock Eng. 2014:1–12. 5. Chang SH, Lee CI, Lee YK. An experimental damage model and its application to the evaluation of the excavation damage zone. Rock Mech Rock Eng. 2007;40:245–285. 6. Chen L, Wang CP, Liu JF, et al. Damage and plastic deformation modeling of Beishan granite under compressive stress conditions. Rock Mech Rock Eng. 2015;48:1623–1633. 7. Qiu SL, Feng XT, Xiao JQ, Zhang CQ. An experimental study on the pre-peak unloading damage evolution of marble. Rock Mech Rock Eng. 2014;47:401–419. 8. Cai M, Kaiser PK, Tasaka Y, Maejima T, Morioka H, Minami M. Generalized crack

J.C. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 84 (2016) 130–141

9.

10. 11. 12. 13. 14. 15. 16. 17. 18.

19.

20. 21.

22.

23.

initiation and crack damage stress thresholds of brittle rock masses near underground excavations. Int J Rock Mech Min Sci. 2004;41:833–847. Alkan H, Cinar Y, Pusch G. Rock salt dilatancy boundary from combined acoustic emission and triaxial compression tests. Int J Rock Mech Min Sci. 2007;44:108–119. Xie N, Zhu QZ, Xu LH, Shao JF. A micromechanics-based elastoplastic damage model for quasi-brittle rocks. Comput Geotech. 2011;38:970–977. Zhou H, Zhang K, Feng XT. A coupled elasto-plastic-damage mechanical model for marble. Sci China Technol Sci. 2011;54:228–234. Hu DW, Zhu QZ, Zhou H, Shao JF. A discrete approach for anisotropic plasticity and damage in semi-brittle rocks. Comput Geotech. 2010;37:658–666. Pensée V, Kondo D, Dormieux L. Micromechanical analysis of anisotropic damage in brittle materials. J Eng Mech. 2002;128:889–897. Molladavoodi H, Mortazavi A. A damage-based numerical analysis of brittle rocks failure mechanism. Finite Elem Anal Des. 2011;47:991–1003. Dragon A, Mroz Z. A continuum model for plastic–brittle behaviour of rock and concrete. Int J Eng Sci. 1979;17:121–137. Mazars J, Pijaudier-Cabot G. Continuum damage theory – application to concrete. J Eng Mech. 1989;115:345–365. Contrafatto L, Cuomo M. A framework of elastic–plastic damaging model for concrete under multiaxial stress states. Int J Plast. 2006;22:2272–2300. Hansen NR, Schreyer HL. A thermodynamically consistent framework for theories of elastoplasticity coupled with damage. Int J Solids Struct. 1994;31:359– 389. Carol I, Rizzi E, Willam K. On the formulation of anisotropic elastic degradation. I. Theory based on a pseudo-logarithmic damage tensor rate. Int J Solids Struct. 2001;38:491–518. Ju JW. On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects. Int J Solids Struct. 1989;25:803–833. Menzel A, Ekh M, Runesson K, Steinmann P. A framework for multiplicative elastoplasticity with kinematic hardening coupled to anisotropic damage. Int J Plast. 2005;21:397–434. Shao JF, Jia Y, Kondo D, Chiarelli AS. A coupled elastoplastic damage model for semi-brittle materials and extension to unsaturated conditions. Mech Mater. 2006;38:218–232. Chiarelli AS, Shao JF, Hoteit N. Modeling of elastoplastic damage behavior of a

141

claystone. Int J Plast. 2003;19:23–45. 24. Mohamad-Hussein A, Shao JF. Modelling of elastoplastic behaviour with nonlocal damage in concrete under compression. Comput Struct. 2007;85:1757– 1768. 25. Jia Y, Song XC, Duveau G, Su K, Shao JF. Elastoplastic damage modelling of argillite in partially saturated condition and application. Phys Chem Earth, Parts A/ B/C. 2007;32:656–666. 26. Chen L, Shao JF, Huang HW. Coupled elastoplastic damage modeling of anisotropic rocks. Comput Geotech. 2010;37:187–194. 27. Wong TF. Micromechanics of faulting in westerly granite. Int J Rock Mech Min Sci Geomech Abstr. 1982;19:49–64. 28. Wu XY, Baud P, Wong TF. Micromechanics of compressive failure and spatial evolution of anisotropic damage in Darley Dale sandstone. Int J Rock Mech Min Sci. 2000;37:143–160. 29. Simo JC, Hughes TJR. Computational Inelasticity. New York: Springer Science Business Media, LLC; 1998. 30. Clausen J, Damkilde L, Andersen L. An efficient return algorithm for non-associated plasticity with linear yield criteria in principal stress space. Comput Struct. 2007;85:1795–1807. 31. Huang J, Griffiths DV. Observations on return mapping algorithms for piecewise linear yield criteria. Int J Geomech. 2008;8:253–265. 32. Hibbitte K. ABAQUS User Subroutines Reference Manual. HKS INC.; 2005. 33. Aubertin M, Li L, Simon R. A multiaxial stress criterion for short- and long-term strength of isotropic rock media. Int J Rock Mech Min Sci. 2000;37:1169–1193. 34. Salari MR, Saeb S, Willam KJ, Patchet SJ, Carrasco RC. A coupled elastoplastic damage model for geomaterials. Comput Method Appl Mater. 2004;193:2625– 2643. 35. Vermeer PA, de Borst R. Non-associated plasticity for soils, concrete and rock. Heron. 1984;29:3–64. 36. Xu WY, Wang RB, Wang W, Zhang ZL, Zhang JC. Creep properties and permeability evolution in triaxial rheological tests of hard rock in dam foundation. J Cent South Univ Technol. 2012;19:252–261. 37. Goodman RE. Introduction to Rock Mechanics. 2nd ed.,John Wiley; 1989:289. 38. Read RS. 20 Years of excavation response studies at AECL Underground Research Laboratory. Int J Rock Mech Min Sci. 2004;41:1251–1275.