Modeling time-dependent deformation behavior of brittle rock using grain-based stress corrosion method

Modeling time-dependent deformation behavior of brittle rock using grain-based stress corrosion method

Computers and Geotechnics 118 (2020) 103323 Contents lists available at ScienceDirect Computers and Geotechnics journal homepage: www.elsevier.com/l...

10MB Sizes 0 Downloads 62 Views

Computers and Geotechnics 118 (2020) 103323

Contents lists available at ScienceDirect

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

Research Paper

Modeling time-dependent deformation behavior of brittle rock using grainbased stress corrosion method Guang Liua,b,d, Ming Caib,c,

T



a

School of Civil Engineering, Hefei University of Technology, Hefei, Anhui 230009, China Geomechanics Research Centre, MIRARCO, Laurentian University, Sudbury, Ontario P3E 2C6, Canada c Key Laboratory of Ministry of Education for Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110004, China d Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering, Ministry of Education, Wuhan University, Wuhan 430072, China b

ARTICLE INFO

ABSTRACT

Keywords: Brittle rock Grain-based model Time-dependent deformation Stress corrosion Static-fatigue tests

This research presents a grain-based stress corrosion model to simulate time-dependent failure of brittle rock. Time-dependent deformation of rock is considered as the manifestation of stress corrosion that may occur inside grains and on grain boundaries. A damage-rate law is adopted and modified to describe continuous degradations of parallel-bond diameters and contact properties in the stress corrosion process. A detailed numerical scheme with self-adaptive determination of time step is used to simulate the time-dependent deformation behavior of rock. The results show that both the short-term mechanical properties and the fatigue test results of Lac du Bonnet granite can be reproduced by the proposed method. In addition, compared with the original parallelbond stress corrosion method, the modified grain-based stress corrosion method predicts the failure time of rock better under low driving-stress ratios. With the increase of the driving-stress ratio, the failure time of the numerical specimen decreases, and the creep curves tend to have a stair-stepping shape. The numerical simulation results show that an increase of 15 MPa confining pressure under a driving-stress ratio of 0.6 can lengthen the failure time of the rock from 139 days to 11 years.

1. Introduction Knowledge of time-dependent deformation behavior of brittle rock is important for predicting the long-term stability of rock structures. The interest in the time-dependent deformation behavior of brittle rock has increased gradually over the last 40 years with the growing demands for resource exploitation and underground construction at depth. Numerous studies have confirmed that long-term stability of rock engineering structures such as rock slopes and tunnels is controlled not only by the stress field but also by time-dependent deformation and damage of rock [1–4]. The macroscopic manifestation that a rock can deform and fail at a constant applied stress is often termed as brittle creep [5]. Creep test has been widely used to study time-dependent deformation of brittle rock [6–10]. In general, the idealized creep behavioral curve can be divided into three regimes (Fig. 1): (1) primary or transient, decelerating creep, (2) secondary, constant strain-rate creep, and (3) tertiary, accelerating creep [11]. Heap et al. [12] presented experimental data of brittle creep of an igneous rock under triaxial stress loading. Their test results indicate that the creep strain rates are highly dependent on the



level of applied stress. With a 20% increase in stress, there are three orders of magnitude increase in the creep strain rate. Therefore, the time needed to conduct creep tests under low driving-stress ratios (i.e., the ratio of the applied stress to the short-term compressive strength) is quite long, which brings a challenge to laboratory testing. Some theoretical models [9,14–18] and numerical simulations [19–22] have thrown light on the mechanism of time-dependent deformation behavior of rock. Debernardi and Barla [23] proposed a nonlinear visco-elasto-plastic creep model that considers creep threshold and long-term strength by connecting an instantaneous elastic Hooke body, a visco-elasto-plastic Schiffman body, and a nonlinear visco-plastic body in series. More recently, an exponential creep model on the basis of material property degradation was proposed by Xu et al. [24] to simulate the rheological behavior of inhomogeneous rock and time-dependent deformation of tunnels. The solid rheology (e.g. lattice distortion, dislocation slip, van der Vaal’s bonds) may be affected by new or pre-existing crack propagation in the time-dependent deformation process [25]. Yang et al. [26] reported that with the increase in damage extent both the cohesion and the internal friction angle of the intact mudstone specimens display an attenuation trend.

Corresponding author at: Geomechanics Research Centre, MIRARCO, Laurentian University, Sudbury, Ontario P3E 2C6, Canada. E-mail address: [email protected] (M. Cai).

https://doi.org/10.1016/j.compgeo.2019.103323 Received 2 May 2019; Received in revised form 9 September 2019; Accepted 23 October 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

Nomenclature α β1, β2 β3, β4 up , low γp γs s

a

1,

c

f

¯ ¯' ¯s ¯n ¯s0 t tip

¯ ¯' ¯s

a

min

3

A c c0 ¯, D ¯' D D¯a D¯ min E F¯ n , F¯ s F¯ n' , F¯ s'

the constant of proportionality between the corrosion rate and the reaction rate two constants in damage-rate law for parallel-bond two constants in damage-rate law for smooth-joint contact the upper boundary and lower boundary of β1 the reducing rate of parallel-bond diameters the damage rate of the degradation factor for smooth-joint contacts rotation at contacts ¯ ¯ ' to initial D the ratio of current D ¯ the ratio of D¯a to D

G I k¯ n , k¯ s kv ¯s M ¯ s' M m ms mi nc p (t ) p0 R¯ , R¯ ' rs rmid

value of when a contact bond breaks in tensile failure value of when a contact bond breaks in shear failure axial and confining stress in triaxial tests uniaxial compressive strength peak axial stress microscopic tensile stress at contacts new microscopic tensile stress after the change of parallelbond diameters tensile strength at contacts microscopic normal stress at contacts initial tensile strength tensile strength crack-tip stress microscopic shear stress at contacts new microscopic shear stress after the change of parallelbond diameters shear strength at contacts internal friction angle degradation factor for smooth-joint contacts the allowable lower boundary of the degradation factor the degradation factor when a smooth-joint contact breaks in tensile failure the degradation factor when a smooth-joint contact breaks in shear failure the minimum degradation factor

T t

t tp ts tf U n, U s V0 v+

Based on this view, it is considered that numerical models that can describe time-dependent deformation of rock better should consider micro-mechanical failure mechanisms such as evolution of microstructures and microcracks [27–29]. Hence, modeling time-dependent

the cross-sectional area of the parallel bond cohesion initial cohesion initial and current bonding diameter the allowable lower boundary of parallel-bond diameters the minimum value of parallel-bond diameters apparent activation energy normal and shear forces at contacts new normal and shear forces after reducing the parallelbond diameter the universal gas constant moment of inertia for parallel-bond normal and shear stiffness per unit area the decreasing rate from up to low moment at contacts new moment after reducing the parallel-bond diameter a parameter expressing the stress level on the contacts micro-activation stress level H–B model parameter for intact rock the number of cycles until the first bonded contact breaks. contact property at time t initial property of the contacts initial and current bonding radius driving-stress ratio the middle driving-stress ratio which corresponds to a half value of up low absolute temperature time increment of time the elapsed time for the failure of next parallel-bond contact the elapsed time for the failure of next smooth-joint contact the elapsed time for the failure of next contact normal and shear displacements at contacts experimental constant activation volume

deformation behavior at the grain-scale level is needed to better understand the long-term strength of rock. There are many mechanisms such as temperature effect, chemical reaction and environmental changes that may influence time-dependent deformations of brittle rock. However, the stress-corrosion process is considered the most likely factor responsible for the time-dependent precursory cracking, deformation and accelerated seismic activities that commonly precede shear ruptures that cause earthquakes [30–32]. Stress corrosion expresses a thermally activated surface reaction that is preferential between a chemically active pore fluid, most commonly water, and the strained atomic bonds close to crack tips [33]. It has been demonstrated that the tensile stress intensity factor plays an important role in determining the crack propagation rate [34,35]. Because stress corrosion is a key factor leading to time-dependent deformation of rock, a method to simulate time-dependent deformation must consider the property degradation induced by stress corrosion. By introducing a damage-rate law to the parallel-bond radius in the Particle Flow Code (PFC), Park [36] and Potyondy [37] employed a parallel-bond stress corrosion (PSC) model to study time-dependent deformation behavior of rock. The PSC model can predict the failure time of rock in a certain range of the driving-stress ratios (about 0.7–0.9). However, the PSC model is less efficient when a confining pressure is applied to a specimen or fatigue in the tensile mode is considered, because this stress corrosion method is based on the parallel-bond contact model which cannot capture well the strength

Fig. 1. Schematic diagram for different phases of time-dependent deformation of rock. Three stages in an idealized creep curve are represented by I, II, and III. Reproduced from Damjanac and Fairhurst [13]. 2

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

envelope and the ratio of compressive to tensile strengths [38]. On the contrary, grain-based model (GBM) is a more advanced model compared with the parallel-bond model and it can reproduce most shortterm macro-mechanical behaviors of rock [39,40]. Unfortunately, few researchers simulate time-dependent deformation behavior using the GBM modeling approach because no explicit numerical schemes in GBM on time-dependent deformation process are readily available in PFC. The purpose of this research is to develop a numerical method to simulate time-dependent deformation of rock using the GBM modeling approach, through stress-induced corrosions occurring in the granule interior and on the grain boundaries. A grain-based stress corrosion method (GSC for short) is proposed in GBM to simulate time-dependent deformation of brittle rock.

and the ratio of uniaxial compressive to tensile strengths [41]. Although the flat joint model can capture the mechanical behavior well with fewer model parameters, GBMs reflect the grain structure better. The GBM modeling approach is particularly suitable for some rocks with polygonal grain structures. In order to study the inter- and intra-grain stress corrosion process, we adopted the GBM modeling approach in this research. Fig. 2a shows the grain structure of a granite [43]. It is seen that the shapes of the grains are quite similar to polygons. Hence, GBMs incorporate this micro-characteristic using polygonal grains (Fig. 2b), with the grain boundaries represented by smooth-joints. Note that the GBM in Fig. 2b is PFC-based. The balls located on the opposite sides of a smooth-joint contact are allowed to overlap and slide along the joint plane when the joint bond breaks. Each grain consists of several particles, which are cemented by parallel-bonds [38]. The GBM has been widely used to study mechanical behaviors of rock [40,44–48]. Steps to generate a GBM specimen in PFC (PFC-GBM) are summarized as follows (Fig. 3):

2. Model framework In this section, a grain-based stress corrosion (GSC) model, which is based on the parallel-bond stress corrosion (PSC) model [36,37], is presented. In the GSC model, the stress corrosion process can occur on the inter-grain and intra-grain contacts. The stress corrosion process on the intra-grain contacts is represented by a modified PSC model. A strength degradation model is developed for the stress corrosion process on the inter-grain contacts.

(a) Generation of a contact network from an initial large disc packing Based on the grain information such as grain sizes, mineral types, and contents of each mineral, a particle assembly is generated. The centers of adjacent discs are connected and a contact network is formed. The contact points are represented by the red dots in Fig. 3.

2.1. Review of grain-based model

(b) Determination of nodes of the grain structure

The grain-based model (GBM) in PFC was proposed by Potyondy [41] to investigate the failure process of brittle rock such as Äspö diorite. This method mimics a rock by a synthetic material that consists of a large number of deformable, breakable polygonal grains cemented along their adjoining sides. In a GBM, each grain includes several bonded particles and the contacts of the grains are depicted by smoothjoint contacts. The smooth-joint contacts are designed to describe the interface behavior regardless of the local particles’ contact orientations. When the smooth-joint contacts are applied, two contacting balls can move along a user-specified joint interface without causing local dilation. One major advantage of the GBM modeling approach is that it can reproduce most macro-mechanical properties of rock with the consideration of microstructures of rock. Because the macro-mechanical properties of rock are influenced by its microstructure [42], it is necessary to replicate the grain structure in numerical specimens to capture the macro-mechanical behaviors including the failure envelope

The centers for each triangle mesh are found and they are denoted as nodes of the grain structure. For polygon meshes, more than one node is found to ensure that every disc is surrounded by several nodes. The nodes are represented by the black dots in Fig. 3. (c) Generation of the grain structure The grain structure is created by linking the two nodes on the opposite sides of the contact points. This grain structure shows the same mean grain sizes as the real grains. However, the geometric topology of the grain structure is random, which has not been compared with the microstructure in real rocks. (d) Covering an assembly of small discs with the grain structure In this procedure, an assembly of small discs is generated in a box.

Fig. 2. Schematic diagram for a GBM. (a) The microscopic scanning image of a thin section of granite (1 – grain boundary cracks (red), 2 – intra-grain cracks (blue), Qtz – Quartz, KFs – K-feldspar, Ab – Albite) [43]; (b) GBM considering microstructures of rock. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 3

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

Fig. 3. A general procedure to generate GBM specimens. (a) Generation of a contact network from an initial large disc packing. The contact network is represented by the black meshes; (b) determination of nodes of the grain structure; (c) generation of the grain structure; (d) covering the assembly of small discs with the grain structure. The red lines indicate the smooth-joint contacts. The black meshes indicate the grain structure. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

initiation and propagation of microcracks on the grain level. The first stage in a typical crack-growth curve of brittle materials such as glasses and ceramics is controlled by a stress-corrosion reaction [49]. Wiederhorn and Bolz [50] provided a rate equation to describe crack propagation based on the static-fatigue theory of Hillig and Charles [51]. In the PSC model, the following equation is used to express the corrosion ¯ ) with the assumption that the rate of the parallel-bond diameter (D stress-corrosion reaction only affects the cement [37].

¯ dD = dt

E GT

V0 e

e

v+ tip GT

(1)

where α is the constant of proportionality between the corrosion rate and the reaction rate, V0 is an experimental constant, E is the apparent activation energy, v+ is the activation volume, tip is the crack-tip stress, G is the universal gas constant, and T is the absolute temperature. Based on the damage-rate relations (Fig. 4) proposed by Potyondy [37], a driving-stress ratio rs is introduced into the damage-rate relations to better simulate time-dependent deformation behaviors at various stress levels. The stress corrosion process is driven by the combined effects of the macroscopic driving-stress ratio and the microscopic tensile stresses at contacts. When the reaction-site stress is taken as ¯ , the damage-rate law for the PSC model can be written as

Fig. 4. Damage-rate relations in PSC model. (a) Schematic diagram for the damage of the parallel-bond diameter; (b) damage rate curve for the parallelbond diameter. Reproduced from Potyondy [37].

The assembly has the same grain size distribution as that in the real specimen. Then this assembly of small discs is covered by the grain structure generated in the previous step using the geometry information of grain structure. The contacts located on the grain boundary are set to the smooth-joint contacts. The discs located in each polygonal mesh are bonded by the parallel-bonds to form a assembly of mineral grains. 2.2. Time-dependent deformation behavior of intra-grain contact

¯ dD = dt

The PSC model considers time-dependent deformation behavior in the formulation of the bonded-particle model (BPM) by adding a damage-rate law to the parallel-bond diameter. Because the intra-grain contact is represented by the parallel-bond model in this research, the time-dependent deformation behavior of intra-grain contact is expressed by a modified PSC model. In essence, time-dependent deformation of rock is influenced by the

0, 1e

¯ < ¯a ¯ 2 ¯s ,

,

¯a ¯

¯ < ¯s ¯s

(2)

where ¯ is the microscopic tensile stress at the contact, ¯s is the parallelbond tensile strength, ¯a is the micro-activation stress, β1 and β2 are two material parameters. In our research, β1 is a function of the drivingstress ratio. To simulate the time-dependent deformation behavior of the 4

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

bonded particles, the parallel-bond diameter is assumed to reduce with ¯ ' ) after a time increment of t can be obtime. The new diameter (D tained by

¯' = D ¯ D

p

t

satisfy the following equations.

F¯ n' 0 0 0 F¯ s' = 0 0 0 3 ¯ s' M

(3)

where p is the reducing rate of the parallel-bond diameter. According to Eq. (2), p is written as p

0,

=

1e

where

¯a

¯ < ¯s

(4)

1

3

f

3

=

up

1 + e kv (rs

(5)

low rmid )

+

¯' D ¯ D

As stated and explained above, stress-induced damage is an essential factor causing time-dependent deformation of rock. As shown in Fig. 6, Martin [52] found that the strength and cohesion of Lac du Bonnet (LdB) granite decrease with the increase of damage. However, the change of the friction angle is relatively complex and it does not show a monotonic decreasing trend with the increase of damage in the LdB granite test data. The decreases of the strength and cohesion are large (more than 40%) and it seems that the curves can be fitted by exponential functions. Based on the above consideration, a property degradation method is proposed to model time-dependent deformation behaviors of inter-grain contact. In grain-based models, inter-grain contacts are defined by the smooth joint model. To model the time-dependent deformation of inter-

where 1 and 3 are the axial and confining stresses (maximum and minimum principal stresses) in the fatigue test, and f is the peak axial stress under a confining pressure of 3 . In our study material parameter β1 is considered as a function of the driving-stress ratio rs (see Fig. 5), which can be written as 1

=

2.3. Time-dependent deformation behavior of inter-grain contact

In the original damage-rate law [38], parameter β1 is a material constant given by users. That damage-rate law is only related to the microscopic tensile stresses at contacts when material parameters β1 and β2 are given, and the microscopic tensile stresses are unknown prior to the simulation. This brings a great challenge for one to match the t–rs curves in the model calibration process. The driving-stress ratio rs can be expressed as

rs =

(10)

is a multiplier describing the degradation degree of ¯ is the ¯ ' is the current bonding diameter, and D bonding diameter, D original bonding diameter.

¯ < ¯a

¯ 2 ¯s ,

F¯ n F¯ s ¯s M

(6)

low

where up and low are the upper and lower boundaries of β1, respectively, rmid is the middle driving-stress ratio that corresponds to the half value of up + low , and k v expresses the decreasing rate from up to low . By introducing the function shown in Eq. (6), the damage rate is now influenced by the driving-stress ratio. To understand the effect of parameters in Eq. (6), we plot a graph for this equation in Fig. 5. As illustrated in Fig. 5a, β1 is controlled by its upper boundary (βup), lower boundary (βlow), middle driving-stress ratio (rmid ) and decreasing rate (k v ). With the increase in the decreasing rate k v , the 1–rs curves show signs of leveling off (Fig. 5b). As a note, parameters βup, βlow and rmid , which have no actual physical meaning in reality, are just three factors that control the shapes of the functions shown in Fig. 5. The forces acting on a parallel-bond in two dimensions can be expressed by the stiffness and the relative bond displacements as

F¯ n F¯ s ¯s M

=

k¯ nA 0 0

0 k¯ sA 0

0 0 k¯ nI

Un Us s

(7)

where , , , , , and are normal force, shear force, moment in the out-of-plane direction, normal displacement, shear displacement, and rotation, respectively; k¯ n and k¯ s are normal and shear stiffness per k¯ nI in Eq. (7) are k¯ sA , unit area, respectively. Therefore, k¯ nA, equivalent to the effective stiffness of normal displacement, shear displacement and rotation, respectively. The area and the moment of inertia for the parallel-bond crosssection are given below (disk thickness )

F¯ n

F¯ s

¯s M

¯ = 2R¯ A = 2Rh I=

2 ¯3 2 R h = R¯ 3 3 3

Un

Us

s

(8) (9)

As the cross-sectional area ( A ) and moment of inertia (I ) are related to the parallel-bond diameter, reducing the diameter of the parallelbond would result in a reduction of the effective stiffness. It should be noted that the displacement components U n , U s and s are the same prior to and after the reduction of the diameter. To keep the displacements constant in Eq. (7), Potyondy [37] states that new forces should

Fig. 5. Variation of β1 with the driving-stress ratio. (a) A specification of parameters needed to calculate β1; (b) influence of decreasing rate (k v ) on β1, with up = 2 × 10 14 , low = 1 × 10 19 , rmid = 0.7 . 5

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

Normalized property

1.0

Mobilized strength Cohesion

0.8 0.6 0.4 0.2 0.0 0.0

0.2

0.4 0.6 Normalized damage

0.8

Fig. 7. Strength envelope for bonded smooth-joint contact.

1.0

level (m ) at the contact. Fig. 8 presents the property degradation curves given by the proposed degradation method when the stress level at a contact decreases linearly with time. The degradation curve calibrated by the proposed model displays the same trend of cohesion degradation shown in Fig. 6. This indicates that the macroscopic performance for the property degradation controlled by microscopic stress corrosion is acceptable. In this research, the properties of smooth-joint contacts, which are assigned by users, include tensile strength, normal and shear stiffness, cohesion, friction coefficient and internal friction angle. It is worth mentioning that the shear strength is not given by users directly. As can be seen from Fig. 7, the shape and the position of the strength envelope is controlled by tensile strength, cohesion, and internal friction angle. In this research, the tensile strength and cohesion are scaled down by the same property degradation factor. Therefore, the strength envelopes can keep the same shape in the stress-corrosion process.

Fig. 6. Variations of mobilized strength and cohesion with damage for LdB granite. The data were obtained from damage-controlled tests (confining pressure 3 = 0 ). Reproduced from Martin [52].

grain contacts, the mechanical properties (such as tensile strength and cohesion) of a smooth joint contact at time t can be expressed as

p (t ) =

(11)

(t ) p0

where p0 represents the original property of the contacts, and is a degradation factor which is a function of time. In this research, the property degradation method is mainly applied to the tensile strength and cohesion of smooth-joint contacts because the variation trend of the tensile strength and cohesion is relatively simple (Fig. 6) and these two parameters have a large influence on the strength of the bonded smooth-joint contacts (Fig. 7). Similar to the PSC model, a damage-rate law for the degradation factor is defined as

0, d dt

=

3e ,

m = max

(

¯ |¯| , ¯s ¯s

3. Numerical implementation

m < ms ms m < 1 m 1

4 m,

)

Potyondy [37] provided an elaborate scheme for the time-advancement procedure to the PSC modeling approach. In our research, the time-advancement procedure is extended to GSC models. In this section, the detailed numerical scheme for GSC models is given and the elapsed time for the next failure contact is derived to adaptively determine the time step.

(12)

where β3 and β4 are two model parameters, m is a parameter expressing the stress level on smooth-joint contacts, ¯s is the shear strength of the smooth joint contact, ¯ is the shear stress acting on the smooth-joint contact, and ms is the micro-activation stress level. 3 is also a function of the driving-stress ratio (similar to Eq. (6)). As shown in Fig. 7, both the shear and tensile stresses can cause the failure of contact bonds. It should be noted that contacts will not merely fail under compressive loading. Thus, stress corrosion only occurs when the tensile stress or shear stress is reaching their corresponding threshold. The stress level m in Eq. (12) incorporates the effect of both tensile and shear stresses. According to Potyondy [37], it is more reasonable to obtain the reaction-site stress by considering both tensile and shear stresses. In the modeling of the time-dependent deformation of smooth-joint contacts, the degradation factor after a time increment of t can be obtained by

(t + t ) =

(t )

s

t

1.0

Normalized property

0.8

0.6

(13)

0.4

where s is the damage rate of the degradation factor. According to Eq. (12), the degradation factor decreases with time, and the damage rate ( s ) depends on the stress level. s can be written as s

=

0, m < ms 4 m , ms e m<1 3

0.0

0.2

0.4 0.6 Normalized time

0.8

1.0

Fig. 8. Variation of property (e.g., tensile strength or cohesion) with time given by the proposed property degradation model. In this figure, rs = 0.8 . 3 is calculated by Eq. (6) with up = 1.5 × 10 14 , low = 1.0 × 10 19 , rmid = 0.7 , k v = 55 , 4 = 27 .

(14)

The property degradation curves may exhibit different trends at different contacts because the damage rate ( s ) depends on the stress 6

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

¯ = ¯=

F¯ n A |F¯ s| A

+

¯ s | R¯ |M I

(16)

¯ s and R¯ are the moment acting in the out-of-plane direction where M and the parallel-bond radius, respectively, and A and I are the area and the moment of inertia for the parallel-bond cross-section, respectively. As the parallel-bond diameter continuously reduces, tensile ( ¯ ) and shear ( ¯ ) stresses change in the stress-corrosion process. Assume that a new parallel-bond diameter can be expressed by D¯ ' = D¯ , we have R¯ ' = R¯ . At time t, the new tensile ( ¯ ' ) and shear ( ¯ ' ) stresses are obtained by substituting Eqs. (8) and (9) into Eq. (16). F¯ n 2R¯

¯' =

(17)

Obviously, the tensile and shear stresses vary with the reduction of the bonding radius. For an intact contact bond, the new tensile and shear stresses should be lower than the tensile ( ¯s ) and shear ( ¯s ) strengths of the contact. Let ¯ ' = ¯s , the ultimate diameter multiplier ( ) can be obtained when the parallel-bond breaks in tension. Similarly, let ¯ ' = ¯s , the ultimate diameter multiplier ( ) can be obtained when the parallel-bond breaks in shear. and can be expressed as [37]:

Fig. 9. Flow chart for modeling time-dependent deformation of brittle rock.

3.1. Time-advancement procedure

F¯ n ±

=

In the rock stress-corrosion calculation, the failure time is extremely long compared with the critical time step in PFC. Consequently, it would be impractical to execute PFC models in the fully dynamic mode, in which each PFC time step would correspond with an identical stresscorrosion time step [37]. The time-advancement procedure maintains two separate accumulated times, one for particle assembly and the other for the stress-corrosion process. Fig. 9 shows a flow chart for modeling time-dependent deformation of brittle rock. At first, the numerical specimen is run to a state of static equilibrium. The condition of static equilibrium is controlled by the user-specified equilibrium ratio limit. Once the model reaches the state of static equilibrium, the stress-corrosion conditions for parallel-bond contacts and smooth-joint contacts are checked. For the parallel-bond contacts, stress-corrosion activates when the maximum tensile stress ( ¯ ) is higher than the micro-activation stress ( ¯a ) but lower than the parallel-bond tensile strength ( ¯s ) ( ¯a ¯ < ¯s ). The activation condition for the smooth-joint contacts, as shown in Eq. (12), is ms m < 1. The time step can be obtained by a self-adaptive procedure or given by users. For the self-adaptive procedure, the time step decreases during the stage of major stress-induced damage formation and increases during the stage of relative quiescence of damage. The time step can be estimated by the following equation:

=

¯ s| ¯s (F¯ n )2 + 24 |M 4R¯ ¯s

|F¯ s| 2R¯¯s

(18) (19)

As mentioned above, the allowable parallel-bond diameter (D¯a ) is expressed as:

D¯a =

(20)

¯ aD

where a is an allowable lower boundary of , which is introduced to prevent contact diameters to reduce to impractical values (e.g., too small). The value of a is assigned by users. For this research, we use the same value ( a = 0.01) suggested by Potyondy [37]. Due to the reduction of the parallel-bond diameter, the maximum tensile ( ¯t' ) and shear ( ¯ ' ) stresses keep changing in the stress-corrosion process, which is shown by Eq. (17). For this reason, a parallel-bond contact may break before its diameter is reduced to D¯a . The critical parallel-bond diameters at tensile and shear failures are given using the multipliers shown in Eqs. (18) and (19). For a parallel-bond contact, its minimum contact diameter denoted by D¯ min can be obtained by Eq. (21). The parallel-bond will break in tension or shear when the contact diameter falls below D¯ min .

D¯ min = max { a ,

tf nc

¯ s| 3|M 2 2R¯ 2

|F¯ s| 2 R¯

¯' =

t=

+

(15)

,

(21)

} D¯

Eq. (21) gives the minimum diameter of a parallel bond. The actual damage rate depends on the local stress and the driving-stress ratio. To obtain the approximate failure time of a parallel bond, it is assumed that the damage rate for the parallel-bond diameter remains constant when the bond diameter is reduced from D¯ to D¯ min . Consequently, the elapsed time (tp ) for the parallel-bond failure can be estimated by

where t f is the estimated time for the next failure contact, and nc is the number of cycles until the bonded contact breaks. 3.2. Estimation of failure time As mentioned above, the elapsed time for the next contact failure is needed for the estimation of the self-adaptive time step. In the stresscorrosion process, a bonded contact may fail in three ways: (1) tensile stress reaches the tensile strength of the contact; (2) shear stress reaches the shear strength of the contact; (3) the bond diameter falls below an allowable value D¯a = a D¯ (for parallel-bond contact) or the strength degrades to an allowable value (for smooth-joint contact). For a parallel-bond contact in two dimensions, the tensile ( ¯ ) and shear ( ¯ ) stresses between the two adjacent particles can be calculated using the beam theory [53].

tp =

, D¯

D¯ min 1

e

¯ < ¯a ¯ 2 ¯s ,

¯a

¯ < ¯s

(22)

On the other hand, the strength of the smooth-joint contact is also reduced in the stress-corrosion process. In this case, the strength degradation method is applied to the cohesion and the tensile strength of smooth-joint contacts. At time t, the cohesion c (t ) and the tensile strength ¯s (t ) are given by

c (t ) = 7

(t ) c0

(23)

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

¯s (t ) =

(24)

(t ) ¯s0

Table 1 Macro-mechanical properties of common minerals and mineral content of granite [40,56,57]

where c0 and ¯s0 are the initial cohesion and tensile strength, respectively, and (t ) is a degradation factor. Similar to the allowable lower boundary of , an allowable lower boundary for the degradation factor (t ) is needed to prevent the strength from reducing to an unreasonable value. The allowable lower boundary of is set to a = 0.01 by Potyondy [37], which means that stress corrosion can occur as long as the condition of micro activation is satisfied. Therefore, the allowable lower boundary for the degradation factor ( a) is also set to 0.01 in our study. According to the strength envelope shown in Fig. 7 (compression is positive), the shear strength ¯s (t ) of a bonded smooth-joint contact can be calculated by its cohesion c (t ) and friction angle .

F¯ n tan A

¯s (t ) = c (t ) +

¯s (t ) (26)

(t ) ¯s0 F¯ n tan A

(t ) c 0 +

(27)

When a smooth-joint contact breaks in tension or shear, the degradation factors are denoted as and , respectively. Solving Eq. (27), the critical values for the degradation factors can be obtained as

= =

F¯ n A ¯s 0 1 |F¯ s| ( c0 A

F¯ n tan A

)

(28)

The minimum degradation factor min

= max {

a,

,

min

is (29)

}

Quartz

Biotite

Elastic modulus (GPa) Poisson’s ratio Density (kg/m3) Content (%)

69.8 0.28 2560 45

88.1 0.26 2630 20

94.5 0.08 2650 30

33.8 0.36 3050 5

LdB granite consists of four minerals: K-Feldspar, Plagioclase, Quartz and Biotite [55]. The macro-mechanical properties of common minerals of granite and their component content are shown in Table 1. In our simulation, the component content of the numerical specimens and the densities of the minerals follow the numbers listed in Table 1. The elastic moduli given in Table 1 provide a valuable reference for the calibration of the parallel-bond and contact moduli of particles. As shown in Fig. 10, the synthetic rock specimen has a width of 50 mm and a height of 125 mm. The ratio of height to width is 2.5, which is in accordance with that of the experimental specimens reported by Martin [52]. This numerical specimen is used to perform uniaxial and triaxial compression tests. The grain structure is represented by the mesh in black. The contents of the minerals in the model are consistent with those shown in Table 1. Note that the topological and statistical properties of the grain structure have not been compared with those of LdB granite. The mean grain size in the model is 2.87 mm, which is the same as the mean grain size of different minerals in LdB granite [55]. In general, an acceptable calibration method of micro-property parameters in PFC models can be conducted using the trial-and-error process to match the macro-properties of the simulated material [58–64]. As a consequence, some micro-property parameters that can

Substituting Eqs. (23), (24) and (25) into Eq. (26) gives Eq. (27). F¯ n A |F¯ s| A

Plagioclase

4.1. Numerical model setup

(25)

¯s (t )

K-Feldspar

and long-term micro-property parameters are calibrated using the short-term static [52] and the static-fatigue test data of LdB granite [54], respectively.

For an intact bonded contact, its tensile and shear stresses are not allowed to go above their corresponding tensile and shear strengths. Therefore, the following conditions should be satisfied for an intact smooth-joint contact: F¯ n A |F¯ s| A

Mineral type

Assuming that the damage rate for the degradation factor remains constant when the degradation factor changes from to min , the elapsed time (ts ) for the bonded smooth-joint contact can be estimated from

ts =

, min 3

¯

m = max( ¯ , s

m < ms

e

4 m,

ms

|¯| ) ¯s

m<1 (30)

The stress corrosion process should adopt a consistent time step for the parallel-bond and the smooth-joint contacts. An optional method to determine a consistent time step by failure time (see Eq. (15)) is using the shorter failure time of tp and ts , i.e.,

t f = min {tp, ts}

(31)

In the simulation, a self-adaptive procedure is adopted to determine the time step, which is similar to the method given by Potyondy [37]. The estimation of failure time is an important procedure to determine the time step. As mentioned above, the time step for the stress-corrosion process can be calculated by Eq. (15) after obtaining the failure time t f . 4. Model validation

Fig. 10. A numerical rock specimen using GBM. Yellow, green, blue, and red areas indicate K-Feldspar, Plagioclase, Quartz, and Biotite, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In this section, a grain-based numerical specimen model of LdB granite is generated. As mentioned above, the model parameters consist of short-term and long-term micro-property parameters. The short-term 8

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

reproduce the macro-properties of a simulated material may not be unique. It should be noted that this has been a major issue in PFC modeling, which obstructs its application in engineering design analysis. Because this research mainly focuses on the modeling of timedependent deformation of rock, the limitation of PFC is not discussed in detail hereon. The determination of parameters in this research follows the widely used trial-and-error process. The main macro-properties that should be captured in the numerical modeling include Young’s modulus, Poisson's ratio, uniaxial compressive strength (UCS), tensile strength, and the ratio of uniaxial compressive to tension strengths. Uniaxial compression, biaxial compression, and direct tension tests are performed using a set of presumed micro-property parameters [42]. The micro-property parameters are considered acceptable when the calculated macro-mechanical properties match the corresponding properties of LdB granite. The main macro-mechanical properties of LdB granite from the laboratory test and numerical calculation are compared in Table 2, showing a good agreement between the two. The Young’s modulus E and the Poisson’s ratio are computed by assuming the plane-strain condition using the stress and strain increments occurring between the start of the test and the point at which one-half of the peak stress has been reached. As shown in Table 2, the macro-mechanical properties, especially the tensile strength, can be well captured by the grain-based model. The ratio of the unconfined compressive strength to the direct tensile strength is 26.6, which is only slightly lower than the mean experiment value of 27.0. This reveals that the grain-based model is capable of modeling the mechanical properties of the rock by replicating the granule microstructure. The calculated macro-mechanical properties listed in Table 2 are all within the range of the experiment data, indicating that the micro-property parameters used in the simulation are reasonable. A comparison of the failure envelopes from the experiment and the PFC-GBM simulation is presented in Fig. 11, showing a good agreement in the confinement range of 0–60 MPa. Because no complete stress–strain curves for the specimens from shallow depth (0–200 m) were given by Martin [52], this research exhibits the comparison of the main macro-mechanical properties listed in Table 2. The Hoek–Brown (H–B) and Mohr–Coulomb (M-C) model parameters given by the regression analysis [65] are listed in Table 3. It is seen that the numerical estimation is close to the experimental estimation. In addition, the Hoek–Brown model parameter mi is close to the one estimated using the UCS test data [66,67]. Potyondy and Cundall [38] found that micro-property parameters calibrated using uniaxial compression test data often gave low strengths under confined conditions because the friction angles of the BPMs are small. On the other hand, PFC-GBM models capture more accurate friction angles by introducing a large number of deformable, breakable polygonal grains. This is seen from the good agreement of the strength envelopes from the numerical simulation with the test data. Because of their advantages over BPMs, PFC-GBMs are used to simulate time-dependent deformation of brittle rocks in the next section.

4.2. Modeling time-dependent deformation of LdB granite Static-fatigue tests are used to study time-dependent deformation of rock [9,68,69]. For example, static-fatigue tests of LdB granite were conducted by Schmidtke and Lajtai [54] using specimens with a height to width ratio of 2.0. During the fatigue tests, environmental conditions such as moisture and temperature were kept constant. In this section, the static-fatigue tests for LdB granite are simulated by the proposed GSC model, using specimens with a height of 100 mm and a width of 50 mm. Fig. 12 shows the fatigue test results from the laboratory test and the numerical simulation. The simulation is performed by considering a combined effect of stress-induced corrosions occurring in granule interiors and on grain boundaries because this combined effect is closer to real stress corrosion. The stress-induced corrosion in granule interiors is described by the modified PSC model, and the stress corrosion on grain boundaries is expressed by the proposed strength degradation model. In the numerical fatigue test, four PFC walls are located on the top, bottom surfaces and the sides of a specimen. The confining pressure in the laboratory fatigue test is 0 MPa. In the numerical test, a very small confining pressure (0.01 MPa) is applied to the lateral loading walls to make the servo system in PFC work. Table 4 presents the long-term micro-property parameters calibrated using the fatigue test data. The values of up and low are small because the damage propagation induced by stress corrosion is quite slow. 1 can be calculated by Eq. (6). The parameters used in the PSC model also are helpful for the parameter determination in our simulation because the damage-rate law used in this research is modified based on the damage-rate law used in the PSC model [37]. In the simulation, the damage rate of grains is controlled by parameters 1 and 2 , and the damage rate of grain boundaries is controlled by parameters 3 and 4 . However, there is a lack of direct laboratory evidence to assist in the determination of the stress-corrosion damage rates inside the grains and on the grain boundaries. As mentioned above, the calibration of model parameters is a trialand-error process to match the macro-properties of the simulated material. In order to reduce the number of parameters and simplify the calibration process, 1 = 3 and 2 = 4 are adopted in this simulation. This assumption may have an influence on the damage rate on the grain level, but its influence on the long-term macro-properties is small because the long-term macro-properties can be captured after the model is calibrated. The calibration process for the parameters listed in Table 4 follows the following three steps. (1) Assume a set of parameters. The objective of calibration is to match the curve of failure time and the driving-stress ratio in the staticfatigue tests. Thus, the choice of these parameters should consider their effect on the shape of the curve. An important criterion is that the stress corrosion rate has a positive correlation with 1 and 2 . 1 decreases with the reduction of the driving-stress ratio. up , low , rmid , and k v are four parameters to control the curve of 1 (see Fig. 5 for details). (2) Perform the numerical static-fatigue test simulation. The static-fatigue tests are simulated under different driving-stress ratios from 1.0 to 0.6. The failure time for each driving-stress ratio is recorded; the failure time is determined at the start of the accelerating creep phase. (3) Plot the curves of failure time and driving-stress ratios. The curves from the numerical simulation and the laboratory test are compared. If the simulation result matches the laboratory data well, the calibration process is stopped. Otherwise, go to step (1) and repeat the calibration process.

Table 2 Comparison of macro-mechanical properties of LdB granite.a Macro-mechanical properties

Experiment

Simulation

Young's modulus, E (GPa) Poisson's ratio, Uniaxial Compressive Strength, c (MPa) Tensile strength b, t (MPa) Ratio of compressive to tensile strengths,

69 ± 5.8 0.26 ± 0.04 200 ± 22 7.4 ± 1.04 27.0

68.3 0.26 209.9 7.9 26.6

c/ t

a The experiment samples from the URL were obtained from shallow depth (0–200 m) [52]. b The tensile strength is the direct tensile strength, which is converted from the Brazilian tensile strength given by Martin [52].

It can be seen from Fig. 12 that the failure time (t fs ) of the specimens increases significantly with the decrease of the driving-stress ratio (rs ). 9

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

800 700

σ1 (MPa)

Table 4 Input parameters for modeling long-term strength of LdB granite.a

H-B model of intact rock 0.5 σ1= σ3+σc(miσ3/σc+1)

Parameters

600

Upper boundary,

500

Lower boundary, rmid Reducing rate, kv

400 300

-10

6.0 × 10

15

low

1.0 × 10 0.7 −55.0 25.8

20

Micro-activation stress (MPa), ¯a Micro-activation stress level, ms

Experimental data PFC-GBM result Test H-B curve Simulation H-B curve

100 0

up

2

200

0

10

20

30

40

50

[37] captures the failure time of LdB granite in fatigue tests well when driving-stress ratios are in the range of 0.7 to 0.9. When the drivingstress ratios are less than 0.7, the PSC model predicts a shorter failure time. It should be noted that the vertical coordinate in Fig. 12 is the logarithm of failure time. A small deviation of curves can cause a large prediction error of failure time when the driving-stress ratio is less than 0.7. For high driving-stress ratios (rs > 0.9 ), the curves from both the PSC and the GSC models approach minus infinity, which means that the time-to-failure is zero. This seems reasonable because the failure time should approach zero theoretically (logarithm of time approaches minus infinity) when the driving-stress ratios approach 1. However, the experimental data do not show this trend. One reason may be due to the fact that the time-to-failure is very short (less than 0.1 s) under very high-stress ratios and it is difficult to measure the time accurately in an experiment. In the experiment, 140 specimens were divided into 10 groups, and group 1 was used to determine UCS by uniaxial compression test [54]. The stress ratios in the fatigue test for other groups were calculated using the UCS determined from group 1. Due to the strength variability among different groups of the rock specimens, this may cause errors in the plot of the experimental data. This is obvious because several stress ratios in the laboratory test data are larger than 1 (see Fig. 12). The laboratory test data include a low-load series (80 data points) and a high-load series (32 data points, rs > 0.9). The failure time under the high-stress ratios is significantly shorter than that under the lowstress ratios. Schmidtke and Lajtai [54] found that the experimental data can be well fitted by eliminating the high-load series data. Fig. 13 compares the simulation results with the laboratory test data after eliminating the high-load series data. It is seen that the simulation

60

Fig. 11. Failure envelopes of LdB granite. The experimental data and fitted Hoek–Brown (H–B) curves are quoted from Martin [52]. H–B model parameters from the experimental data are: mi = 28.9, c = 210.0 MPa. H–B model parameters from the numerical simulation are: mi = 26.5, c = 209.9 MPa. Table 3 Comparison of model parameters of Hoek–Brown and Mohr–Coulomb failure criteria. Parameters

Experimental estimation

Numerical estimation

mi c (MPa) Cohesion c (MPa) Internal friction angle

28.9 210.0 44.1 49.9

26.5 209.9 49.6 47.7

(°)

8

Lab data Exp-fit Lab data GSC PSC

6

11.6 d

4

2.8 h

2

1.7 min

80.0 0.4

a 1 is not shown in this table because it is not an independent parameter. It can be calculated by Eq. (6).

σ3 (MPa)

Log(t) (s)

Values

7 0.32 y

1s

0.5

Lab data Exp-fit Lab data GSC PSC

6 11.6 d

0.6

0.7

0.8

Driving-stress ratio

0.9

5 27.8 h

1.0

Log(t) (s)

0

Fig. 12. Comparison of failure times under different driving-stress ratios from simulations with laboratory test data of LdB granite [54]. The laboratory test data include a low-load series and a high-load series. The result by the PSC model is obtained from Potyondy [37].

4 2.8 h 3 16.7 min 2 1.7 min 1 10 s

In the simulation, the specimens were loaded first to an axial stress according to the specified driving-stress ratio (rs ). Then the stress-corrosion process was started, with a constant axial stress at a confining pressure. As shown in Fig. 12, the laboratory test data can be fitted by an exponential function: log(t fs) = 1118.2e 8.2277rs + 0.70239 , with R2 = 0.98. As can be seen in Fig. 12, the results from our model (GSC) match the trend line of the laboratory data well. The PSC model by Potyondy

0 0.55

0.60

0.65

0.70

0.75

0.80

0.85

0.90

Driving-stress ratio Fig. 13. Comparison of simulation results with laboratory test data after eliminating the high-load series data (32 points). The Exp-fit function is log (tfs) = 129.24e−4.3298rs + 1.9108. 10

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

results by the GSC model match the laboratory test data well. The failure time predicted by the PSC model is significantly lower than the laboratory data when the driving-stress ratios are less than 0.65. The GSC model predicts the failure time well when the stress ratios are low. The prediction by the GSC model also has a good precision when the stress ratios are higher than 0.75 (Fig. 13). Fig. 14 presents the crack distributions in the failed specimens under different driving-stress ratios. It is found that most cracks are tensile cracks and the failure mode of the specimens is similar to axial splitting. It seems that the cracking paths under high driving-stress ratios show more branches. Under high driving-stress ratios, multiple fractures intersect in the specimens and the crack distribution is more scattered. The failure time is considerably reduced with the increase of the driving-stress ratio. The crack number is counted at the beginning of axial loading. The fatigue test is started when the required constant stress is reached. The initial crack number in the fatigue test is not zero because of the initial damage induced by axial loading. The creep curves in Fig. 15 resemble the idealized creep behavior shown in Fig. 1. With the increase of the driving-stress ratio, the creep curves tend to have stair-stepping shapes. Under a high driving-stress ratio, a specimen reaches the crack damage stress, which implies the start of unstable crack propagation [70]. Unstable propagation of cracks may cause intensive strain energy release. That may be a reason why the creep curves show stair-stepping shapes under high driving-stress ratios. A noteworthy phenomenon is that the crack number curves embody the same variation with the axial strain curves. This indicates that the propagation, connection and coalescence of microcracks are important micro-mechanisms of the creep deformation of brittle rock. With the reduction of the driving-stress ratio, the crack propagation becomes stable, leading to smooth creep curves.

With the increase of confining pressure, the failure time increases drastically. It should be noted that the Y-axis in Fig. 16a is the logarithm of failure time. Although the curves for different confinements are almost parallel to each other in Fig. 16a, Fig. 16b shows that there is a clear effect of confinement on the slope of the time-to-failure curve. Under a driving-stress ratio of 0.6, the failure time for a specimen under 0 MPa confining pressure is about 139 days, while the failure time leaps to 11 years when the confining pressure is 15 MPa. This reveals that the long-term stability of rock can be improved by imposing lateral restrictions. It should be noted that under low confinement the specimen will break in a few seconds once the driving-stress ratio exceeds 0.9. The difference in failure time under different confining pressures is small when the driving-stress ratio is greater than 0.9. The crack pattern under different confining pressures also has an influence on the failure time. This may be the reason why two data points for 5 MPa confining pressure are below the curve of 0 MPa in Fig. 16b. As shown in Fig. 17, the accelerating creep stage in the idealized creep behavioral curves is in a short time window under low confining pressures. However, the accelerating creep stage under a high confining pressure is longer because the lateral restraints increase the load capacity of the rock. The ultimate creep strains for confining pressures of 0, 5, 10, and 15 MPa are 0.03, 0.12, 0.12, and 0.21%, respectively, and the failure times of the specimens are 3.4, 5.0, 15.5 and 39.1 min (see Fig. 16), respectively. The ultimate creep strains can be calculated by subtracting the initial axial strain from the axial strain at the failure of the specimen. It is clear that the confining pressure increases the creep strain and failure time significantly. This reveals that the lateral restraint can improve the load carrying capacity of rock structures during long-term loading.

Fig. 14. Failure modes of specimens under different driving-stress ratios. The crack distributions under different driving-stress ratios 0.6, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, and 0.95 are shown in (a), (b), (c), (d), (e), (f), (g) and (h), respectively. The black and red line segments represent tensile and shear cracks, respectively. rs is the driving-stress ratio, and t is the failure time for numerical specimens. The crack width is denoted by the line width. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 11

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

Fig. 15. Creep curves under different driving-stress ratios under 0 MPa confinement. The driving-stress ratios are 0.6, 0.7, 0.8 and 0.9 for (a), (b), (c) and (d), respectively.

Fig. 16. Failure time under different driving-stress ratios and different confining pressures. (a) The logarithm of failure time and driving-stress ratios under different confining pressures; (b) linear failure time and driving-stress ratios under different confining pressures.

degradation method is widely used to model time-dependent deformation behavior [4], and it is much simpler to be implemented in numerical simulations. Therefore, for the sake of simplicity, this research adopts the bonding diameter degradation method for parallelbond contacts and proposed a property degradation method for the smooth-joint contacts in GBMs. The objective of this research is to provide a solution to modeling time-dependent deformation behavior of brittle rock with a consideration of material microstructure. Although many efforts have been devoted to this issue, the method presented above is just one of numerous possible solutions to the problem. Because there is a lack of experimental research on grain-level stress corrosion, further research is needed to explore which type of degradation method is better.

5. Discussion For a typical assembly using GBM in PFC, two contact models, i.e., parallel-bond and smooth-joint models, are incorporated. In this research, the stress corrosion process is described by different degradation methods for the two contact models: one is the bonding diameter degradation for parallel-bond contacts, and the other is the property degradation for smooth-joint contacts. Bonding diameter degradation, which was suggested by Potyondy [37], conforms to the physical mechanism of microcrack propagation on the grain level. However, this degradation method results in a change of effective stiffness on the parallel-bond cross section, and it involves the adjustment of microscopic contact force component (see Eq. (10)). In general, the property 12

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai

Axial strain

0.006

and the Open Fund of Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering, Ministry of Education, Wuhan University, China.

0 MPa 5 MPa 10 MPa 15 MPa

References

0.005

[1] Bieniawski ZT. Time-dependent behaviour of fractured rock. Rock Mech Rock Eng 1970;2(3):123–37. [2] Hashiba K, Fukui K. Index of loading-rate dependency of rock strength. Rock Mech Rock Eng 2014;48(2):859–65. [3] Sone H, Zoback MD. Time-dependent deformation of shale gas reservoir rocks and its long-term effect on the in situ state of stress. Int J Rock Mech Min Sci 2014;69:120–32. [4] Tran-Manh H, Sulem J, Subrin D. Progressive degradation of rock properties and time-dependent behavior of deep tunnels. Acta Geotech 2016;11(3):693–711. [5] Cruden DM. A theory of brittle creep in rock under uniaxial compression. J Geophys Res 1970;75(17):3431–42. [6] Özşen Hakan, Özkan İhsan, Şensöğüt Cem. Measurement and mathematical modelling of the creep behaviour of Tuzköy rock salt. Int J Rock Mech Min Sci 2014;66:128–35. https://doi.org/10.1016/j.ijrmms.2014.01.005. [7] Yang SQ, Jing HW, Cheng L. Influences of pore pressure on short-term and creep mechanical behavior of red sandstone. Eng Geol 2014;179:10–23. [8] Mishra Brijes, Verma Priyesh. Uniaxial and triaxial single and multistage creep tests on coal-measure shale rocks. Int J Coal Geol 2015;137:55–65. https://doi.org/10. 1016/j.coal.2014.11.005. [9] Fabre G, Pellet F. Creep and time-dependent damage in argillaceous rocks. Int J Rock Mech Min Sci 2006;43(6):950–60. [10] Yang SQ, Hu B. Creep and long-term permeability of a red sandstone subjected to cyclic loading after thermal treatments. Rock Mech Rock Eng 2018;51(10):2981–3004. [11] Jaeger JC, Cook NG. Fundamentals of rock mechanics. London: Methuen & Co. Ltd; 1969. p. 513. [12] Heap MJ, Baud P, Meredith PG, Vinciguerra S, Bell AF, Main IG. Brittle creep in basalt and its application to time-dependent volcano deformation. Earth Planet Sci Lett 2011;307(1–2):71–82. [13] Damjanac B, Fairhurst C. Evidence for a long-term strength threshold in crystalline rock. Rock Mech Rock Eng 2010;43(5):513–31. [14] Shao JF, Zhu QZ, Su K. Modeling of creep in rock materials in terms of material degradation. Comput Geotech 2003;30(7):549–55. [15] Okubo S, Fukui K. An analytical investigation of a variable-compliance-type constitutive equation. Rock Mech Rock Eng 2005;39(3):233–53. [16] Zhao LY, Zhu QZ, Xu WY, Dai F, Shao JF. A unified micromechanics-based damage model for instantaneous and time-dependent behaviors of brittle rocks. Int J Rock Mech Min Sci 2016;84:187–96. [17] Zheng Hong, Feng Xia-Ting, Hao Xian-jie. A creep model for weakly consolidated porous sandstone including volumetric creep. Int J Rock Mech Min Sci 2015;78:99–107. https://doi.org/10.1016/j.ijrmms.2015.04.021. [18] Abu Al-Rub RK, Darabi MK. A thermodynamic framework for constitutive modeling of time- and rate-dependent materials. Part I: theory. Int J Plasticity 2012;34:61–92. [19] Lu Yinlong, Elsworth Derek, Wang Lianguo. A dual-scale approach to model timedependent deformation, creep and fracturing of brittle rocks. Comput Geotech 2014;60:61–76. https://doi.org/10.1016/j.compgeo.2014.04.001. [20] Li WJ, Han YH, Wang T, Ma JW. DEM micromechanical modeling and laboratory experiment on creep behavior of salt rock. J Nat Gas Sci Eng 2017;46:38–46. [21] Sharifzadeh M, Tarifard A, Moridi MA. Time-dependent behavior of tunnel lining in weak rock mass based on displacement back analysis method. Tunn Undergr Space Technol 2013;38:348–56. [22] Li X, Konietzky H, Li XB. Numerical study on time dependent and time independent fracturing processes for brittle rocks. Eng Fract Mech 2016;163:89–107. [23] Debernardi D, Barla G. New viscoplastic model for design analysis of tunnels in squeezing conditions. Rock Mech Rock Eng 2009;42(2):259–88. [24] Xu T, Tang CA, Zhao J. Modeling of rheological deformation of inhomogeneous rock and associated time-dependent response of tunnels. Int J Geomech 2011;12(2):147–59. [25] Paraskevopoulou C, Perras M, Diederichs M, Loew S, Lam T, Jensen M. Time-dependent behaviour of brittle rocks based on static load laboratory tests. Geotech Geol Eng 2017;1–40. [26] Yang SQ, Tian WL, Jing HW, Huang YH, Yang XX, Meng B. Deformation and damage failure behavior of mudstone specimens under single-stage and multi-stage triaxial compression. Rock Mech Rock Eng 2018;52(3):673–89. [27] Bikong C, Hoxha D, Shao JF. A micro-macro model for time-dependent behavior of clayey rocks due to anisotropic propagation of microcracks. Int J Plast 2015;69:73–88. [28] Chen YF, Wei K, Liu W, Hu SH, Hu R, Zhou CB. Experimental characterization and micromechanical modelling of anisotropic slates. Rock Mech Rock Eng 2016;49(9):3541–57. [29] Li X, Huai Z, Konietzky H, Li XB, Wang Y. A numerical study of brittle failure in rocks with distinct microcrack characteristics. Int J Rock Mech Min Sci 2018;106:289–99. [30] Main IG, Meredith PG, Sammonds PR. Temporal variations in seismic event rate and b-values from stress corrosion constitutive laws. Tectonophysics 1992;211(1–4):233–46. [31] Main IG, Meredith PG. Stress corrosion constitutive laws as a possible mechanism of intermediate-term and short-term seismic quiescence. Geophys J Int

0.004

0.003 0.0

0.2

0.4

0.6

0.8

1.0

Normalized time Fig. 17. Creep curves under different confining pressures. The loading time is normalized by the failure time. The driving-stress ratio is 0.8.

6. Conclusions This research presents a grain-based stress corrosion (GSC) approach to simulate time-dependent deformation of brittle rock. The method considers the time-dependent deformation of rock as the manifestation of stress corrosion that may occur inside the grains and on the grain boundaries. By introducing a damage-rate law for parallelbond diameters and contact properties with the consideration of the driving-stress ratio, which is defined by the ratio of current deviatoric stress to the peak deviatoric stress, a detailed numerical scheme with a self-adaptive determination of time step is presented to simulate timedependent deformation of rock. The short-term static strength data of LdB granite were used to calibrate the grain-based numerical specimen and the static-fatigue tests of LdB granite were simulated by the proposed method. It is shown that the GSC method captures both the shortterm and the long-term mechanical behaviors of the rock well. The GSC model predicts the failure time of rock under both low and high drivingstress ratios. Both the numerical simulation and the laboratory test results indicate that the driving-stress ratio has a large influence on failure time. The failure time of specimens increases with the decrease of the driving-stress ratio. The creep behavior is influenced by the microcracking process in the fatigue tests. The specimens show unstable crack propagations under high driving-stress ratios, which results in stair-stepping creep curves. The failure mode of specimens under long-term loading is similar to axial split. Most cracks are caused by tensile failure and the spatial distribution of cracks seems to be more scattered under high driving-stress ratios. Although the difference of failure time under different confining pressures is small (a few seconds) when the drivingstress ratios are high, this difference becomes very large when the driving-stress ratios are low. The numerical simulation results show that an increase of confining pressure from 0 to 15 MPa under a drivingstress ratio of 0.6 can lengthen the failure time of the rock from 139 days to 11 years. This reveals that long-term stability of rock structures can be greatly improved by imposing lateral restrictions. The GSC model is a natural extension of the GBM for modeling time-dependent deformation of brittle rock. Future work will focus on modeling time-dependent failure of blocky rocks. Acknowledgements Financial supports from the Natural Sciences and Engineering Research Council (NSERC, CRDPJ 485174 - 15) of Canada, Rio Tinto, Helca, MIRARCO, and Northeastern University in China for this work are gratefully acknowledged. The first author also acknowledges the support by Anhui Natural Science Youth Fund (No. 1908085QE216) 13

Computers and Geotechnics 118 (2020) 103323

G. Liu and M. Cai 1991;107(2):363–72. [32] Heap MJ, Baud P, Meredith PG. Influence of temperature on brittle creep in sandstones. Geophys Res Lett 2009;36(19). [33] Brantut N, Heap MJ, Meredith PG, Baud P. Time-dependent cracking and brittle creep in crustal rocks: a review. J Struct Geol 2013;52:17–43. https://doi.org/10. 1016/j.jsg.2013.03.007. [34] Meredith PG, Atkinson BK. Stress corrosion and acoustic emission during tensile crack propagation in Whin Sill dolerite and other basic rocks. Geophys J Int 1983;75(1):1–21. [35] Lajtai EZ, Bielus LP. Stress corrosion cracking of Lac du Bonnet granite in tension and compression. Rock Mech Rock Eng 1986;19(2):71–87. [36] Park N. Discrete element modeling of rock fracture behavior: fracture toughness and time-dependent fracture growth [PhD thesis] The University of Texas at Austin; 2006. [37] Potyondy DO. Simulating stress corrosion with a bonded-particle model for rock. Int J Rock Mech Min Sci 2007;44(5):677–91. [38] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41(8):1329–64. [39] Liu G, Sun WC, Lowinger SM, Zhang ZH, Huang M, Peng J. Coupled flow network and discrete element modeling of injection-induced crack propagation and coalescence in brittle rock. Acta Geotech 2019;14(3):843–68. [40] Lan HX, Martin CD, Hu B. Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading. J Geophys Res: Solid Earth 2010;115(B1). [41] Potyondy DO. A grain-based model for rock: approaching the true microstructure. Proceedings of rock mechanics in the Nordic Countries. 2010. [42] Liu G, Cai M, Huang M. Mechanical properties of brittle rock governed by microgeometric heterogeneity. Comput Geotech 2018;104:358–72. [43] Rong G, Peng J, Cai M, Yao M, Zhou C, Sha S. Experimental investigation of thermal cycling effect on physical and mechanical properties of bedrocks in geothermal fields. Appl Therm Eng 2018;141:174–85. [44] Bahrani N, Kaiser PK. Estimation of confined peak strength of crack-damaged rocks. Rock Mech Rock Eng 2016;50(2):309–26. [45] Chen W, Konietzky H, Tan X, Frühwirt T. Pre-failure damage analysis for brittle rocks under triaxial compression. Comput Geotech 2016;74:45–55. [46] Wang X, Cai M. Modeling of brittle rock failure considering inter- and intra-grain contact failures. Comput Geotech 2018;101:224–44. [47] Peng J, Wong LNY, Liu G, Teh CI. Influence of initial micro-crack damage on strength and micro-cracking behavior of an intrusive crystalline rock. Bull Eng Geol Environ 2018;78(4):2957–71. [48] Wong LNY, Zhang Y. An extended grain-based model for characterizing crystalline materials: an example of marble. Adv Theory Simulat 2018;1(8):1800039. [49] Freiman SW. Effects of chemical environments on slow crack growth in glasses and ceramics. J Geophys Res Solid Earth 1984;89(B6):4072–6. [50] Wiederhorn S, Bolz L. Stress corrosion and static fatigue of glass. J Am Ceram Soc 1970;53(10):543–8. [51] Hillig WB, Charles RJ. Surfaces, stress-dependent surface reactions, and strength.

High Strength Mater 1965:682–705. [52] Martin CD. The strength of massive Lac du Bonnet granite around underground openings PhD thesis Winnipeg, Canada: University of Manitoba; 1993. [53] Itasca Consulting Group I. Manual PFC2D (Particle Flow Code), Version 4.0 Users’ Guide. Minneapolis, Minnesota, USA: Itasca, Minneapolis; 2008. [54] Schmidtke RH, Lajtai E. The long-term strength of Lac du Bonnet granite. Int J Rock Mech Min Sci Geomech Abstracts 1985:461–5. [55] Eberhardt E, Stimpson B, Stead D. Effects of grain size on the initiation and propagation thresholds of stress-induced brittle fractures. Rock Mech Rock Eng 1999;32(2):81–99. [56] Mavko G, Mukerji T, Dvorkin J. The rock physics handbook: tools for seismic analysis of porous media. Cambridge University Press; 2009. [57] Read RS. Interpreting excavation-induced displacements around a tunnel in highly stressed granite. PhD Thesis. Winnipeg: University of Manitoba, 1994. [58] Rong G, Liu G, Hou D, Zhou CB. Effect of particle shape on mechanical behaviors of rocks: a numerical study using clumped particle model. Sci World J 2013;2013. [59] Liu G, Rong G, Peng J, Zhou CB. Numerical simulation on undrained triaxial behavior of saturated soil by a fluid coupled-DEM model. Eng Geol 2015;193:256–66. [60] Yoon J. Application of experimental design and optimization to PFC model calibration in uniaxial compression simulation. Int J Rock Mech Min Sci 2007;44(6):871–89. [61] Fakhimi A, Villegas T. Application of dimensional analysis in calibration of a discrete element model for rock deformation and fracture. Rock Mech Rock Eng 2007;40(2):193–211. [62] Bahrani N, Kaiser PK, Valley B. Distinct element method simulation of an analogue for a highly interlocked, non-persistently jointed rockmass. Int J Rock Mech Min Sci 2014;71:117–30. [63] Wong LNY, Peng J, Teh CI. Numerical investigation of mineralogical composition effect on strength and micro-cracking behavior of crystalline rocks. J Nat Gas Sci Eng 2018;53:191–203. [64] Wang X, Cai M. A comprehensive parametric study of grain-based models for rock failure process simulation. Int J Rock Mech Min Sci 2019;115:60–76. [65] Shen J, Wan L, Zuo J. Non-linear Shear Strength Model for Coal Rocks. Rock Mech Rock Eng 2019;52(10):4123–32. [66] Cai M. Practical estimates of tensile strength and Hoek-Brown strength parameter mi of brittle rocks. Rock Mech Rock Eng 2010;43(2):167–84. [67] Cai M. A simple method to estimate tensile strength and Hoek-Brown strength parameter mi of brittle rocks. Proceedings of the 3rd CANUS rock mechanics symposium, Toronto. 2009. [68] Chen Wei, Konietzky Heinz. Simulation of heterogeneity, creep, damage and lifetime for loaded brittle rocks. Tectonophysics 2014;633:164–75. https://doi.org/10. 1016/j.tecto.2014.06.033. [69] Lajtai EZ, Schmidtke RH. Delayed failure in rock loaded in uniaxial compression. Rock Mech Rock Eng 1986;19(1):11–25. [70] Cai M, Kaiser PK, Tasaka Y, Maejima T, Morioka H, Minami M. Generalized crack initiation and crack damage stress thresholds of brittle rock masses near underground excavations. Int J Rock Mech Min Sci 2004;41(5):833–47.

14