Numerical direct shear tests to model the shear behaviour of rock joints

Computers and Geotechnics 51 (2013) 101–115

Contents lists available at SciVerse ScienceDirect

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

Numerical direct shear tests to model the shear behaviour of rock joints M. Bahaaddini a,b,⇑, G. Sharrock a, B.K. Hebblewhite a a b

School of Mining Engineering, The University of New South Wales, Sydney, Australia Shahid Bahonar University of Kerman, Kerman, Iran

a r t i c l e

i n f o

Article history: Received 19 December 2012 Received in revised form 18 February 2013 Accepted 19 February 2013 Available online 26 March 2013 Keywords: Direct shear test Rock joint Particle flow code Bonded particle model Smooth joint model

a b s t r a c t In this paper, the shear behaviour of rock joints are numerically simulated using the discrete element code PFC2D. In PFC, the intact rock is represented by an assembly of separate particles bonded together where the damage process is represented by the breakage of these bonds. Traditionally, joints have been modelled in PFC by removing the bonds between particles. This approach however is not able to reproduce the sliding behaviour of joints and also results in an unrealistic increase of shear strength and dilation angle due to the inherent micro-scale roughness of the joint surface. Modelling of joints in PFC was improved by the emergence of the smooth joint model. In this model, slip surfaces are applied to contacts between particles lying on the opposite sides of a joint plane. Results from the current study show that this method suffers from particle interlocking which takes place at shear displacements greater than the minimum diameter of the particles. To overcome this problem, a new shear box genesis approach is proposed. The ability of the new method in reproducing the shear behaviour of rock joints is investigated by undertaking direct shear tests on saw-tooth triangular joints with base angles of 15°, 25° and 35° and the standard joint roughness coefficient profiles. A good agreement is found between the results of the numerical models and the Patton, Ladanyi and Archambault and Barton and Choubey models. The proposed model also has the ability to track the damage evolution during the shearing process in the form of tensile and shear fracturing of rock asperities. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The mechanical behaviour of jointed rock masses is strongly dependent on the mechanical and geometric properties of discontinuities. At surface and shallow depth, the rock mass behaviour is predominately controlled by shear displacement along discontinuities rather than failure of intact rock [1]. Surface roughness is known to have a significant effect on the shear behaviour of rock discontinuities. Much has been done to develop empirical and analytical shear strength criteria to account the effect of roughness. Patton [2] and Ladanyi and Archambault [3] were among the first to develop a shear strength criterion for rough rock joints. Patton [2] by undertaking experiments on saw-tooth triangular asperities proposed a bilinear shear strength criterion. In this model, the failure envelope consists of two linear segments intersecting at a normal stress rT, called the transition stress. For normal stresses rn less than rT, the shear strength s is governed by sliding along the joint asperities with an inclination angle of i and friction angle

⇑ Corresponding author at: School of Mining Engineering, The University of New South Wales, Sydney, Australia. Tel.: +61 2 93856118; fax: +61 2 9313 7269. E-mail addresses: [email protected] (M. Bahaaddini), g.sharrock@ unsw.edu.au (G. Sharrock), [email protected] (B.K. Hebblewhite). 0266-352X/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2013.02.003

/l. Above the transition stress, the shear strength is governed by shearing of asperities with cohesion cj and residual friction angle /r.

s ¼ rn tanð/l þ iÞ for rn < rT s ¼ cj þ rn tan /r for rn P rT where

ð1Þ ð2Þ

j rT ¼ ðtan /lctan . /r Þ

Although this analytical model captures two key mechanisms of asperity sliding and shearing, it has several drawbacks such as the lack of quantitative description of joint cohesion cj which is dependent on intact rock properties and asperity geometry. Furthermore, for real rock joints shearing and sliding take place simultaneously [4,5]. Based on the energy principles, Ladanyi and Archambault [3,6] developed a failure model which recognises the simultaneous sliding and shearing mechanisms.

s¼

rn ð1  as Þðv_ þ tan /l Þ þ as SR 1  ð1  as Þv_ tan /l

ð3Þ

where as is the shear area ratio, v_ is dilation rate and SR is the intact rock strength. Ladanyi and Archambault [3,6] suggested that the intact rock strength SR can be obtained from the Fairhurst intact strength criterion [7]:

102

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

Shear Strength (τ)

Patton’s model for sliding

results of tilt, push or pull tests. They also proposed that the peak dilation angle can be estimated as:

Fairhurst’s model for intact material

dn ¼ Ladanyi and Archambault’s model

M¼

)

Fig. 1. Ladanyi and Archambault model (modified from [8]).

Table 1 Micro-properties of balls and bonds. Particle micro-mechanical properties

Parallel bonds micro-mechanical properties

Ball density (kg/m3)

2205

Ec (GPa) Coefficient of friction kn/ks

2.8 0.6 1.45

SR ¼ r c

pffiffiffiffiffiffiffiffiffiffiffiffi   nþ11 rn 0:5 1þn n rc

Ec (GPa) Normal strength (MPa) Shear strength (MPa) s  n =k k

2.8 20 ± 4.5 20 ± 4.5 1.45

ð4Þ

where rc is the uniaxial compressive strength of intact rock and n is the ratio of tensile to uniaxial compressive strength of the intact rock. This model is shown in Fig. 1. Ladanyi and Archambault [6] based on extensive experimental study on saw-tooth triangular joints, proposed the following empirical power laws for estimation of as and v_ :

  rn k1 as ¼ 1  1 

ð5Þ

rT

v_ ¼

  rn k2 1 tan i

ð6Þ

rT

The parameter rT is the transition stress at which the strength of joint is equal to the strength of intact rock and k1 and k2 are empirical constants with values of 1.5 and 4.0, respectively. Although this method is based on sound theoretical concepts, there are several shortcomings, such as the difficulty to determine rT and i for irregular joint geometries and incorrect kinematic consideration of the shearing of irregular joints [4,9]. Barton and Choubey [10], based on extensive experimental studies on real rock joints, proposed the following empirical shear strength model:



s ¼ rn tan /r þ JRC log

  JCS

rn

ð8Þ

where M is a damage coefficient, takes values of 1 or 2 for shearing under low or high normal stress, respectively [11]. This parameter can also be approximated from the following relation [10]:

Residual strength of planar joints

Normal Stress (

  1 JCS JRC log M rn

ð7Þ

where /r is residual friction angle and JRC and JCS are joint roughness coefficient and joint compressive strength, respectively. For unweathered rock joints, the residual friction angle /r is equal to the basic friction angle /b and JCS is equal to the uniaxial compressive of intact rock rc. JRC can be determined either from visual comparison with 10 standard JRC profiles or by back calculation of the

JRC þ 0:70 12 logðJCS=rn Þ

ð9Þ

Several empirical and analytical methods have been proposed in the last few decades to estimate the shear properties of rock joints in direct shear tests. However, a robust method that can consider the effect of joint roughness on the shear behaviour of rock joints is not still available. One part of this problem relates to the lack of suitable method for quantitative representation of the joint roughness. Another part is related to unknown mechanism of asperities degradation during the shearing process. The rapid progress of computer technology has provided the opportunity to numerically simulate the shear behaviour of rock joints. Several researchers have tried to simulate the shear behaviour of rock joints using the constitutive models [12–14]. However, these approaches are unable to trace the asperity degradation and development of cracks inside the intact material. Furthermore, special difficulties exist in modelling the geometries of real rock joints. Karami and Stead [15] by using a hybrid FEM/DEM code, ELFEN software, tried to investigate the process of joint surface degradation in the direct shear test. They carried out direct shear test on three standard JRC profiles of 0–2, 10–12 and 18–20 and studied the joint surface damage degradation in terms of asperity wear and propagation of tensile cracks in intact material along the joint plane. However, the ability of their method in reproducing the shear strength and dilation angle was not studied quantitatively. Particle flow code (PFC) is a promising numerical approach to study the mechanical behaviour of discontinuities explicitly. PFC is a discrete element software in which the intact material is represented by a composite of spherical (in 3D) or circular (in 2D) particles bonded together, called bonded particle model (BPM). The damage process in the BPM is represented by breakage of bonds and propagation and coalescence of these micro-cracks control the macroscopic behaviour of the BPM. Therefore, it provides a scientific tool to track the asperities degradation during the shearing process and investigate the effect of joint roughness on the shear behaviour of rock joints. The BPM also has the ability to reproduce many features of intact rock behaviour [16–21]. The main aim of this paper is to investigate the shear behaviour of rock joints in a direct shear test using PFC2D [22]. The most commonly used method for modelling joints in PFC is the bond removal approach. In this method particles lying on the joint surface are left unbonded. Shortcomings of this method are presented in Section 3. A new approach to model the joints in PFC is the smooth joint (SJ) model, in which slip surfaces are applied at contacts between particles that lie on the opposite sides of the joint plane. A detailed analysis of the SJ model in reproducing the shear

Table 2 Model output for calibrated intact rock (macro-properties). UCS (MPa)

E (GPa)

m

Experimental

Average

27.4

4.2

0.2

Numerical

Average SD CV (%)

27.45 0.59 2.15

4.25 0.01 0.23

0.198 0.007 0.38

SD: Standard Deviation, CV (%) = (Standard Deviation/Average)  100.

103

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

Normal Force Wall 2 xvel = 0.1m/s

40 mm

Wall 1 Loading wall xvel = 0.1m/s

Box gap 1.2 mm

Wall 3 Fixed wall 100 mm Fig. 2. Simulation of a direct shear test under constant normal load.

Fig. 3. Generation of a planar joint by bond removal method (joint particles are shown in green colour). (For interpretation of colors in this figure legend, the reader is referred to the web version of this article.)

behaviour of rock joints is presented in Section 4. Investigations presented in this paper show that this approach suffers from particle interlocking at large shear displacement. To overcome this problem a new method is presented and the ability of proposed method in reproducing the shear behaviour of rock joints is investigated by undertaking direct shear tests on saw-tooth triangular joints and standard JRC profiles.

2. Numerical model set up 2.1. Bonded particle model calibration Two basic bond models are provided in PFC; namely the contact bond (CB) and parallel bond (PB) models. Both CB and PB models approximate the physical behaviour of cement like substance lying between particles and joining them together. The PB is active over a finite rectangular (2D) or circular (3D) cross-section lying on the contact plane and centred at the contact point and can transmit both force and moment between particles while the CB is active only at the contact point and can only transmit force acting at the contact point. In the PB model, bond breakage results in immediate reduction of contact stiffness while in the CB model, the contact stiffness is active as long as the particles stay in contact

[22,23]. In previous numerical modelling of direct shear tests in PFC, the CB model has been predominately used due to its reduced number of micro-parameters compared to the PB model [24–26]. However, it is well understood that the PB model produces a more realistic representation of rock-like materials [23,27]. Therefore, the PB model is used in the current study. The proper selection of the micro-properties for the particles and bonds is an important step in bonded particle modelling. In order to generate a numerical model to simulate the real rock, the micro-mechanical properties of the models are calibrated against physical experiments on Hawkesbury sandstone. The Hawkesbury sandstone dominates the Sydney region in Australia and hosts many significant civil and mining excavations. The strength of the Hawkesbury sandstone is known to vary considerably depending on deposition type, the degree of cementation and weathering. An extensive dataset of intact rock properties of Hawkesbury sandstone exists at the University of New South Wales [28]. Based on the procedure developed by Potyondy and Cundall [29], BPM samples having a width of 42 mm and height of 84 mm are generated. The size of the particles are controlled by predefined minimum and maximum diameters (Dmin = 0.28 mm and Dmax = 0.42 mm) which satisfies a uniform particle-size distribution. Inverse calibration is used to calibrate the micro-properties of the bonds and the particles. This is involved an iterative process

104

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

are compared with physical experiments in Table 2. The numerical tests results are in good agreement with the physical experiments and the low values of coefficients of variation (CV) of the numerical experiments indicate the minor difference between the results of numerical tests. 2.2. Numerical modelling of direct shear test under constant normal load In direct shear tests, shear box specimens with width 100 mm and height 40 mm are generated. Samples consist of around 38,000 particles and their micro-properties are presented in Table 1. A schematic of the direct shear test under constant normal load (CNL) is presented in Fig. 2. The normal load is applied vertically to the upper block and this load is kept constant during the test by using a servo-control mechanism [22]. The lower block is kept stationary during the test and a horizontal velocity of 0.1 m/ s is applied to the upper block. Calculations in PFC are based on an explicit time stepping algorithm. The time step (Dt) in each calculation has small value of Dt  4.139  108 s. This means that the upper block moves horizontally at the rate of 4.139  1010 m/ time step. Sensitivity studies conducted as a part of this study have shown that this displacement rate is slow enough to ensure that the sample remains in quasi-static equilibrium. The shear displacement is measured by tracing the horizontal displacement of the upper block and the shear stress is calculated by dividing the recorded the reaction force on wall 3 by the joint length. Normal displacement is measured by the vertical displacement of the upper wall.

Fig. 4. Effect of particle coefficient of friction on (a) shear stress and (b) dilation of planar joint.

to select the micro-parameters to reproduce the desired macroproperties measured in the laboratory experiments. It is common to calibrate micro-parameters against uniaxial compressive strength UCS, elastic modulus E and Poisson’s ratio v. The first step involves the calibration of E which is controlled by the particle contact modulus Ec,particle normal/shear stiffness ratio kn/ks, parallel bond modulus Ec and parallel bond normal/shear stiffness ratio s . The second step is matching v which is influenced by k /k n =k k n s n =k s and these micro-parameters are calibrated in an iterative and k process with the first step. Finally, the UCS is calibrated by altering the normal and shear strengths of parallel bonds [22]. The microparameters of the calibrated BPM model are listed in Table 1. Ten uniaxial compression tests on are carried out on BPM models with different randomly generated particle assemblies and mechanical properties listed in Table 1. The results of numerical experiments

3. Bond removal method One of the first attempts to model the shear behaviour of a rock joint using PFC2D was carried out by Cundall [30]. In this model, for generation of the joint, the particle contacts within a specific distance (i.e. Dmax) on either side of the joint track were left unbounded. Subsequently, this method has been used by a number of researchers for simulation of the direct shear test in PFC [24– 26,31–34]. A generated planar joint using this method is shown in Fig. 3. Unbonded particles are shown in green colour. Previous studies have shown that the coefficient of friction of joint particles has considerable effect on the shear behaviour of joints and therefore need to be set close to zero to reproduce a realistic shear strength [24,26]. Cundall [30] suggested that the coefficient of friction of joint particles should be calibrated against the friction angle of the joint plane. In order to investigate the effect of the coefficient of friction of joint particles on the shear strength and dilation, a sensitivity analysis is undertaken. Planar joints are generated and the coefficient

Fig. 5. Stress concentration at micro-scale asperities generated the by bond removal method (black line: contact force, red line: tensile crack, blue line: shear crack). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

4.5 4

τ = 0.58σn + 0.93 R2 = 0.996

Triangular

Shear Stress (MPa)

3.5 Planar

3 2.5

τ = 0.43σn + 0.89 R2 = 0.995

τ = 0.96σn + 0.33 R2 = 0.997

2 1.5 1 0.5

τ = 0.82σn + 0.28 R2 = 0.992

0 0

1

2

3

4

5

Normal Stress (MPa) Fig. 6. Peak shear strength envelopes of planar and saw-tooth triangular joint with the base angle of 15° using the bond removal method.

of friction of joint particles l is varied from 0.005 to 0.15. Direct shear tests under constant normal stress of 1 MPa are carried out and the results are presented in Fig. 4. Results show unrealistic shear behaviour for planar joints. Reduction of the coefficient of friction results in a decrease in the peak shear strength and the normal displacement but generation of the joint by the bond removal approach leads to a micro-scale roughness at the joint surface (Fig. 3). Although the average size of particles is very small compare to the length of shear box (0.35 mm to 100 mm), due to the circular shape of the particles, their non-uniform distribution and unequal size, the joint surface is not smooth at the initial stage of the shear test. Therefore, the local shear resistance occurs for riding up of the interlocking particles. Due to stress concentration at these interlocking particles, some bonds break and when the particle jumps over another particle, the stored energy converts to kinetic energy and energy dissipation takes place. As a result, the shear stress drops after reaching the peak shear stress, unlike real planar joints where the shear stress is approximately constant at post peak. The reduction of the coefficient of friction only results in a decrease of the local shear resistance for riding up of interlocking particles. There is also significant dilation in the initial stages of shearing which is due to the riding up of interlocking particles over each other. Fig. 5 shows the contact force distribution during the shearing process. Due to the micro-scale roughness of the joint plane, contact forces are not distributed uniformly over the joint plane and there are unrealistic stress concentrations at interlocking particles. As the coefficient of friction increases, stress concentration at these micro-scale asperities increases and results in increased breakage of adjacent joint particles bonds. Cundall [30] noted that the unbonded band (green particles in Fig. 3) should have a width of several particle diameters; otherwise large resisting forces may occur between the opposing particles on

105

either side of the joint track. Previous studies have shown that the increase of band width results in a reduction in the peak shear strength and normal displacement, but does not lead to realistic shear behaviour of planar joints [35]. Furthermore, free movement of several particles across the joint track reduce the efficiency of this method, especially for modelling rough joints. In order to investigate the ability of the bond removal method to reproduce the shear strength of rock joints, direct shear tests on a saw-tooth triangular joint with a base angle of h = 15° and a planar joint are carried out. The coefficient of friction of joint particles is assumed 0.005. The peak shear strengths of these joints under different normal stresses are presented in Fig. 6. Bilinear peak strength envelopes are found for both saw-tooth triangular and planar joints. It is well understood that the shear strength of planar joints obeys the Coulomb criterion (s = rntan /l). Therefore, its bilinear peak shear strength envelope is related to the existing micro-scale roughness of the joint surface. There are also unrealistic cohesions of 0.28 MPa and 0.33 MPa at low normal stresses for planar and saw-tooth triangular joints, respectively. Furthermore, there is a significant difference in the slope of the peak shear strength envelopes at high normal stresses. The results clearly reveal that the bond removal method is unable to simulate the shear behaviour of rock joints. Generation of joints using this approach results in the inherent micro-scale roughness and recognition of the shear strength is caused by the interlocking particles or the joint asperities becomes complicated, especially for rough joints.

4. Smooth joint model The smooth joint model is a newly developed method for simulation of rock discontinuities in PFC [36]. Limited efforts have been carried out to apply this approach for simulation of direct shear test [37,38]. In this model, rectangular slip surfaces (i.e. smooth joint) are assigned to all contacts between particles lying on the opposite sides of the joint plane. Particle pairs intersected by a smooth joint may overlap and pass through each other rather than forced to move around one another [39,40]. A typical example of the smooth joint model is shown in Fig. 7. The joint plane consists of two coincident planar surfaces (indicated as surface 1 and surface 2). The orientation of these surfaces ^ j , where n ^ j points into is defined by the joint unit normal vector n surface 2. The unit normal vector is defined by the joint dip angle hp:

^ j ¼ ðsin hp ; cos hp Þ n

ð10Þ

The smooth joint is assigned at a contact between particles their centres lying on opposite sides of the joint plane. In order to deter^ j and mine in which surface each particle lies, the dot product of n ^ c is used. The n ^ c of each ball is directed contact unit normal vector n ^j  n ^ c P 0, the from its centre to the centre of connected ball. If n ball lies on surface 1. For example in Fig. 7, ball 1 lies on surface

y

x θp

Surface 2 Surface 1

Joint Plane

Fig. 7. Application of smooth joint (modified from [22]).

106

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

| | |

|

Unload/reload

|

(a)

|

(b)

Fig. 8. Smooth joint force–displacement law (a) Normal force versus normal displacement and (b) Shear force versus shear displacement [22].

2. The smooth joint is applied in the direction parallel to the joint plane and remains active while there is a nonzero overlap between particles [22]. By creation of the smooth joint, the bond between particles (if exist) is removed and a set of elastic spring is applied uniformly over a rectangular cross section. The area of smooth joint cross section A is:

A ¼ 2Rt

ð11Þ

where t is the thickness (t = 1.0) and R is the radius of the smooth joint cross section.

R ¼ k minðRA ; RB Þ

ð12Þ

where  k is radius multiplier, which is usually taken 1.0, and RA, RB are particles radii. The properties of smooth joints are defined by  , shear stiffness k  and coefficient of fricthe SJ normal stiffness k nj sj tion lj. Let force F and displacement U be expressed in the local coordinate system of the joint plane. ^tj is the tangential unit vector of ^ j and ^tj points in the posijoint plane where the cross product of n ^ Therefore, ^ j  ^tj ¼ þk). tive z direction (i.e. n

^ j þ Fs F ¼ Fnn ^ j þ Us U ¼ Un n

particle model and a constant normal stress of 1 MPa is applied to the upper block. The results of a representative direct shear test are shown in Fig. 9. It is found that this model reproduces the shear behaviour of planar joints when the shear displacement is less than the minimum particle diameter (Dmin = 0.28 mm). However, when the shear displacement exceeds Dmin, there is an unrealistic increase in the shear stress and normal displacement. The contact force distribution at the initial stage of the test and after 0.3 mm shear displacement is presented in Fig. 10. At the initial stage of testing, stresses are distributed uniformly along the joint plane while after 0.3 mm shear displacement, stresses are concentrated at lock up points. These stress concentrations lead to an increase in the shear stress and associated initiation and propagation of cracks at the lock up points. In order to find the reason of these

ð13Þ ð14Þ

where Fn and Un are the normal force and displacement and Fs and Us are the shear force and relative shear displacement vectors, respectively. Positive values of Fn and Un denotes compression and overlap, respectively. The strength model of the smooth joint contact is that of the Coulomb sliding model. In each time step, the relative displacement increment between two particles intersected by the smooth joint is decomposed to normal and tangential components to the joint surface. The normal and shear forces are updated via (Fig. 8):

 ADU e F n :¼ F n þ k nj n  ADUe F0 :¼ Fs  k s

sj

s

ð15Þ ð16Þ

where D denotes an increment. An updated shear force value F0s is computed at the end of each time step. The greatest possible value of shear force F s is dependent on the smooth joint coefficient of fric  tion lj F s ¼ F n lj . If jF0s j 6 Fs , then jFs j ¼ jF0s j ; otherwise sliding takes place and jFs j ¼ F s (Fig. 8b). 4.1. Simulation of direct shear test by the smooth joint model In order to study the ability of the SJ model to reproduce the shear behaviour of joints, a direct shear test on a planar joint is undertaken. A planar joint plane is inserted into the bonded

Fig. 9. Results of direct shear test on planar joint using the smooth joint model.

107

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

Particle interlocking

Surface 2 Surface 1

(a)

(b)

Fig. 10. Particle interlocking in the smooth joint model (a) initial stage and (b) after 0.3 mm shear displacement.

Fig. 11. Procedure of shear box genesis (a) vessel generation, (b) assembly generation, (c) isotropic stress installation, (d) elimination of floaters, (e) parallel bond installation and (f) application of smooth joints.

stress concentrations, the process of smooth joint application is investigated at the micro-scale (Fig. 10). At the initial stage, all ball-ball contacts are checked and smooth joints are assigned at contacts where the balls lie on opposite sides of the joint plane. Centre of ball 5803 is above the joint plane (lie on surface 2) and the smooth joints are assigned to its contacts with balls 6910 and 7191 that lie on surface 1. By shear movement of the upper block to the right side, ball id 5803 passes through ball id 7219 and a new contact is generated with ball 3068. The newly generated contacts are checked each cycling step. Due to the applied normal stress to the upper block, ball 5803 moves downward during the shearing process and its centre lies below the joint plane (on surface 1). The newly generated contact of this ball (ball

5803) and ball 3068 is checked and both balls lie on the same side of the joint plane and hence a smooth joint is not assigned at this contact. This results in the interlocking of ball 5803 and the stress concentration around this particle. This particle interlocking leads to the unrealistic increase of shear stress. This result shows that the SJ model cannot reproduce the shear behaviour of joints at shear displacements larger than Dmin.

5. Shear box genesis approach As discussed in Section 4.1, the inability of the smooth joint model to reproduce the shear behaviour of joints is related to the

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

detection of the top and bottom particles of the shear box blocks and the proper application of the smooth joint. To overcome this problem, a shear box genesis method is proposed. In this approach the procedures of sample preparation and smooth joint application are modified. Upper and lower blocks of the shear box are generated separately and the smooth joint is only applied at contacts between the particles of these blocks. The following six steps are employed in the shear box genesis approach (Fig. 11). 1. Material vessels generation: Each upper and lower block consists of four frictionless walls (Fig. 11a). In order to generate a rough joint in the numerical model, the rock joint profile is digitised and its nodes are used to generate the joint surface walls. In PFC, only one side of the standard wall is active and a wall is defined by the order of its nodes [22]. Therefore, two walls are required to produce a joint profile. 2. Initial assembly generation: Assemblies of randomly placed particles are generated and the upper and lower vessels filled separately (Fig. 11b). The particle diameters are controlled by the predefined Dmin = 0.14 mm and Dmax = 0.21 mm and satisfy uniform particle-size distribution. The number of particles for each block is determined by their volumes and the overall porosity. The normal stiffness of walls is made 10% higher than the average normal stiffness of particles to ensure the ball-wall overlap remains small. Then, the system is allowed to rearrange and reach the static equilibrium under zero friction of particles [22,29]. 3. Applying a specified isotropic stress: In order to reduce the magnitude of locked-in stresses, the radii of all particles are reduced uniformly to reach a specified isotropic stress (rt0 ¼ 1% of the uniaxial compressive strength) [29]. Five measurement circles are placed in each block (Fig. 11c). Isotropic stress r0 ¼ r11 þ2 r22 inside these measurement circles are calculated for each step. If the normalised difference between the specified isotropic stress and the measured  t isotropic  stress is less than r r the isotropic stress tolerance 0rt 0 6 0:5 , iteration stops and 0 floating particles are checked, otherwise the iteration continues [22,29]. 4. Elimination of floating particles: The assembly of randomly distributed and non-uniform size particles may have a large number of particles with less than three contacts (red1 particles in Fig. 11d). Before installing bond between particles, it is desirable to reduce the number of these floating particles and produce a dense pack of highly interlocked grains [22,29]. Therefore, in this step floaters of upper and lower blocks are eliminated and assemblies are generated in which all particles away from the boundaries have at least three contacts. 5. Parallel bond installation: Parallel bonds are installed at contacts between particles that have a separation less than 106 times of average radius of two particles (black lines in Fig. 11e). 6. Applying smooth joints: By separate generation of the upper and lower blocks, particles of these blocks can be recognised easily in each stage of test which solves the problem of detection of upper and lower particles of shear box blocks. After generation of upper and lower blocks, the interface walls are removed. A specified constant normal stress is applied to the upper block and allowing the assembly to relax. This is done by stepping the model until the static equilibrium is achieved. During the relaxation process (and also shearing process), new contacts are generated between the particles of the upper and lower blocks. Smooth joint is applied at newly generated contact if the contact lies between the particles defining the upper and

1

For interpretation of color in Fig. 11, the reader is referred to the web version of this article.

lower blocks (green lines in Fig. 11f). The orientations of smooth joints are parallel to the predefined joint profile which is set to the stationary lower block. This approach resolves the problem of particle interlocking and can be easily employed for direct shear tests of irregular joint profiles. 6. Validation study In order to investigate the ability of the shear box genesis approach in reproducing the shear behaviour of rock joints, a validation study is carried out on saw-tooth triangular joints and standard JRC profiles where the results of the numerical models are compared against well-known shear strength criteria. In order to generate an objective model, smooth joint parameters are calibrated before undertaking the validation study. 6.1. Calibration of smooth joint parameters To calibrate the smooth joint (SJ) parameters, experimental normal deformability and direct shear tests are carried out on sawn planar joints of Hawkesbury sandstone. The SJ normal stiffness is calibrated against experimental normal deformability test and the SJ shear stiffness and SJ coefficient of friction are calibrated against direct shear tests. Normal deformability experiments involve loading of singlejointed rectangular blocks of Hawkesbury sandstone with length 100 mm, width 100 mm and height 40 mm and recording the normal force and normal displacement (based on the procedure developed by Bandis et al. [41]). A typical result of an experimental test on intact rock and rock with planar joint is shown in Fig. 12. As the applied normal stress increases, joint normal stiffness increases. The numerical direct shear tests are carried out under the constant normal stress of 0.5–5 MPa. Therefore, the secant normal stiffness

12 Intact rock 10

Normal Stress (MPa)

108

Rock with planar joint 8 6

5 MPa

4

Kn

2 0.5 MPa

0 0

0.1

0.2

0.3

0.4

Normal Displacement (mm) Fig. 12. Experimental normal deformability test (Kn: normal stiffness of the system).

Table 3 Calibrated smooth joint parameters and result of experimental and numerical tests. Calibrated SJ parameters Normal stiffness (MPa/mm) Shear stiffness (MPa/mm) Friction coefficient, lj (tan(/j))

200 50 0.78

Results of experiments Parameter

Physical experiments

Numerical experiments

System normal stiffness (MPa/mm) System shear stiffness (MPa/mm) Friction angle (°)

28.771 6.42 37.62

28.766 6.417 37.81

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

109

Normal stress = 1 MPa

Normal stress = 2 MPa

Fig. 13. Numerical direct shear tests on planar sawn joints under different normal stresses.

Normal stress = 3 MPa

Normal stress = 5 MPa

Fig. 15. Effect of normal stress on asperity degradation of saw-tooth triangular joint with the base angle of 15 degrees after 3 mm shear displacement (red: tensile crack, blue: shear crack). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. Results of numerical direct shear tests on saw-tooth triangular joint with the base angle of 15 degrees, (a) shear stress versus shear displacement and (b) normal displacement versus shear displacement.

of the system in that range is calculated (Kn = 28.77 ± 2.42 MPa/ mm), as shown in Fig. 12. Based on the procedure explained in Section 5, numerical samples with width 100 mm and height 40 mm (the same dimensions as the physical specimens) are generated with planar interface. In the initial step, side walls and interface walls are removed and the assembly is allowed to reach the static equilibrium. The specimen is loaded by a pair of opposing frictionless walls. The top wall acts as a loading platen and the bottom wall is stationary. These walls

are the assembly vessel which is used for sample generation. The normal stiffness of walls is set equal to the average normal stiffness of contacting particles [22]. The vertical load is applied at a constant displacement rate of 0.01 m/s. The loading continues to half of the uniaxial compressive strength of the intact material. The normal stress and normal displacement are recorded during the test and the secant normal stiffness of system in the range of 0.5–5 MPa is calculated. The inverse calibration method is used to calibrate the SJ normal stiffness, in which SJ normal stiffness is varied until the normal stiffness of experimental tests is achieved. Result of this calibration is presented in Table 3. To calibrate the SJ shear stiffness and SJ coefficient of friction, experimental direct shear tests on planar joints of Hawkesbury sandstone are undertaken. First, numerical direct shear test under constant normal stress of 1 MPa is carried out and the SJ shear stiffness is varied to reproduce the system shear stiffness measured in laboratory experiments. Then, numerical direct shear tests under different normal stresses are undertaken and SJ coefficient of friction is calibrated. Calibrated SJ parameters and results of numerical and experimental tests are presented in Table 3. Results of direct shear test under different normal stresses are depicted in Fig. 13. This results show that the proposed model can reproduce the shear behaviour of planar joints. 6.2. Direct shear test on saw tooth triangular joints Three saw-tooth triangular joint profiles with the base angle h = 15°, 25° and 35° are generated and direct shear test under

110

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

Patton Ladanyi & Archambault Numerical model

6 5

Peak Shear Strength (MPa)

Peak Shear Strength (MPa)

7

4 3 2 1 0 1

2

3

4

Patton Ladanyi & Archambault

10

Numerical model

8 6 4 2 0 0

5

1

2

3

4

Normal Stress (MPa)

Normal Stress (MPa)

(a)

(b)

Peak Shear Strength (MPa)

0

12

12

5

Patton Ladanyi & Archambault

10

Numerical model

8 6 4 2 0 0

1

2

3

4

5

Normal Stress (MPa)

(c) Fig. 16. Investigation of the suggested model’s ability to reproduce the peak shear strength of saw-tooth triangular joints, (a) h = 15°, (b) h = 25° and (c) h = 35°.

different normal stresses are carried out. Results of the direct shear tests on the saw-tooth triangular joint with the base angle of 15° are shown in Figs. 14 and 15. As the normal stress increases, the peak shear strength increases. When the normal stress is less than 3 MPa, shearing takes place by sliding along the joint. At a normal stress of 1 MPa, when the shear stress reaches the peak value, plastic behaviour starts and sliding takes place along the joint (Fig. 14a). Only limited number of cracks is observed during shearing (Fig. 15). As the normal stress increases, the stress concentrations at asperities increase and more tensile cracks are observed. At a normal stress of 3 MPa, asperity wearing is observed in the model. For normal stresses higher than 3 MPa, asperities are sheared-off and the shear stress after reaching to the peak stress drops to residual value. As the normal stress increases, the peak shear displacement dpeak (shear displacement at peak shear strength) and shear stiffness increase. The increase of normal stress also results in a decrease in normal displacement. Results of direct shear tests on saw-tooth triangular joints under different normal stresses are presented in Fig. 16. These results are compared with Patton [2] and Ladanyi and Archambault [3,6] models which were developed on triangular joint profiles. Due to lack of data for the joint cohesion cj, comparison of numerical models against Paton model is carried out for only the sliding mechanism (Eq. (1)). The Ladanyi and Archambault model is derived using the mechanical properties of Hawkesbury sandstone. The transition pressure, rT, is obtained 26.3 MPa using Mogi [42] approach. Goodman [43] suggested that in the absence of sufficient data, this transition pressure can be selected as the uniaxial compressive of intact rock (rc = 27.4 MPa). The selected transition pressure is close to the Goodman proposed value. At the asperity base angle of 15°, when the normal stress is less than 3 MPa, shearing takes place in sliding mode and the peak shear strength of the numerical models are in good agreement

with Patton model. When the normal stress exceeds 3 MPa, shearing of asperities occurs and the peak shear strength of numerical

θ = 15°

θ = 25°

θ = 35°

Fig. 17. Asperity degradation of saw-tooth triangular joints under the constant normal stress of 2 MPa (shear displacement = 3 mm).

111

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

profiles corresponding to different JRC values. They suggested that the JRC value can be estimated by visual comparison to these standard profiles. In order to investigate the ability of the suggested model in reproducing the shear strength and dilation angle of real rock joint profiles, numerical direct shear tests on standard JRC profiles are carried out and the results are compared against Barton and Choubey [10] model. The standard JRC profiles were originally obtained by using a profile comb which can measure profile points every 0.5 to 1 mm. Therefore, these standard JRC profiles are digitised with a horizontal spacing of 0.5 mm. Ten series of shear box samples are generated and direct shear tests under different normal stresses are carried out. Results of the direct shear test on JRC 16–18 are presented in Figs. 18 and 19. As the normal stress increases, the shear stiffness increases and the maximum shear strength is achieved at a higher shear displacement. At low normal stresses the second order asperities have an influence and limited number of fracturing occurs on steep-sided roughness features (Fig. 19). As the normal stress increases, shearing of first order and second order asperities occurs and first order asperities take over the shearing behaviour. Therefore, the peak shear strength increases and dilation rate decreases (Fig. 18). Stress concentration at some asperities results in penetration of tensile cracks inside the intact rock material (Fig. 19, normal stress = 5 MPa). Results of the comparison between the peak shear strength of numerical models and Barton and Choubey model for 10 standard JRC profiles are presented in Fig. 20. The JRC values of joint profiles are assumed to be the experimental back-calculated value of Barton and Choubey for each standard profile and the JCS is the

Normal stress = 0.5 MPa

Fig. 18. Results of numerical direct shear tests on JRC 16–18 (a) shear stress versus shear displacement and (b) normal displacement versus shear displacement.

models falls below the Patton sliding model. At the asperity base angle of 25°, transition from sliding to shearing mode takes place at normal stress of 2 MPa while at the asperity base angle of 35°, the transition stress is reduced to 1 MPa. Good agreement is found between the results of numerical models and Paton and Ladanyi and Archambault models in reproducing the peak shear strength. Effect of the asperity angle on the joint surface degradation under the constant normal stress of 2 MPa is shown in Fig. 17. As the asperity angle increases, stress concentration at triangular asperities increases and results in greater damage to the asperities. At triangular base angle h = 15°, shearing takes place by sliding along the asperities and only limited number of cracks occurs in model. At h = 25°, the first asperity is worn out but the shearing is controlled by sliding along other asperities. At h = 35°, the high stress concentration at the asperities results in the asperity shearing-off. This shows that at low asperity base angles, sliding is the dominant shear mechanism while for steeper asperities crushing and asperity shearing is the controlling mechanism. It is found that asperity degradation occurs by initiation of tensile cracks. Propagation of these tensile cracks results in the asperities shearing.

Normal stress = 1 MPa

Normal stress = 3 MPa

Normal stress = 5 MPa

6.3. Standard JRC profiles The most widely used criterion for the estimation of the shear strength of rock discontinuities is Barton and Choubey [10] model. They carried out extensive series of experiments on various rock joint profiles and published a set of 10 standard joint roughness

Fig. 19. Effect of normal stress on asperity degradation of JRC 16–18 after 3 mm shear displacement.

112

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

Fig. 20. Investigation of the suggested model’s ability to reproduce the peak shear strength of standard JRC profiles.

uniaxial compressive strength of intact rock (rc = 27.4 MPa). Results of the numerical models shows that the peak shear strength is slightly overestimated for JRC 6–8 and 8–10 while the peak shear

strength for other JRC profiles are in good agreement with Barton and Choubey model. The results show that this method can reproduce the peak shear strength.

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

113

Fig. 21. Investigation of the suggested model’s ability to reproduce the peak dilation angle of standard JRC profiles.

Normal displacement is recorded during the shearing process and the peak dilation angle is calculated at the peak shear displacement dpeak. The peak dilation angles of numerical models are compared against Barton and Choubey model, as shown in Fig. 21. For JRC profiles of 0–2 and 2–4, the peak dilation angle is very small in both the numerical and Barton and Choubey models. For other joint profiles, a good agreement is found between the numerical models and Barton and Choubey model. The peak dilation angle is underestimated at JRC 4–6 and JRC 6–8 and overestimated at JRC 12–14 and JRC 18–20 under low normal stress, while for other JRC profiles satisfactorily estimated. These results show that the shear box genesis approach has the ability to reproduce the shear strength and the peak dilation angle of real rock joint profiles. Joint surface degradation during the shearing process for JRC 18–20 under the normal stress of 4 MPa is shown in Fig. 22. In

the pre-yield stage (before shear displacement of 0.5 mm) only minor cracks occur in the model and shear displacement takes place by internal deformation of the intact rock and joint surface. When the shear stress exceeds 80% of the peak shear strength, signs of yielding in the shear stress curve are observed by a departure from linear elastic behaviour. The surface damage is observed to commence when the yield stress is passed. From the yield stress, reduction of effective area results in stress concentration at asperities. The stress at some asperities exceeds the bonds normal strength of particles and tensile cracks occur at these asperities. From the peak stress, tensile cracks develop at these asperities which lead to shear strength reduction. From the peak shear stress to residual shear stress the number of cracks increases rapidly. As the shear displacement continues, further crushing of the joint asperities takes place.

114

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115

Fig. 22. Asperity degradation during shearing process for standard JRC 18–20 under normal stress of 4 MPa (a) shear stress, normal displacement and cumulative crack number during the shearing process and (b) corresponding model state (green: smooth joint, red: tensile crack and blue: shear crack). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

7. Conclusions In this paper, the shear behaviour of rock joints in direct shear tests are numerically simulated using PFC2D. Previously, joints were modelled in PFC using two approaches, namely the bond removal method and the smooth joint model. It is found that simulation of joints using the bond removal method results in an unrealistic increase of shear strength and normal displacement due to inherent micro-scale roughness of the joint surface. The smooth joint model also suffers from particle interlocking which occurs at shear displacement greater than the minimum diameter of the particles. In this study, the shear box genesis approach is presented to resolve this problem. The ability of this method in reproducing the shear behaviour of saw-tooth triangular joints and standard JRC profiles is demonstrated. Results of the numerical models are compared against well-known empirical models and good agreement is found between them. The asperity degradation during the shearing process is studied numerically. It is found that in the pre-yield stage only minor cracking takes place and when the shear stress reaches the yield shear stress, reduction of the effective area results in stress concentration and occurrence of tensile cracks at asperities. Propagation

of these tensile cracks and crushing of the asperities leads to reduction of shear stress to residual stress and rapid increase in the number of cracks.

Acknowledgements The authors thank Dr. Xavier Garcia, Dr. David Potyondy and Dr. Matthew Pierce from Itasca Consulting group and Dr. Dion Weatherley from the University of Queensland for their invaluable technical help and comments.

References [1] Hoek E. Shear strength of discontinuities. Practical rock engineering; 2007. p. 14. [2] Patton FD. Multiple modes of shear failure in rock. In: 1st ISRM congress. Lisbon, Portugal; 1966. p. 509–15. [3] Ladanyi B, Archambault G. Simulation of shear behavior of a jointed rock mass. In: The 11th US rock mechanics symposium (USRMS). Berkeley, CA; 1969. p. 105–25. [4] Kodikara JK. Shear behaviour of rock-concrete joints and side resistance of piles in weak rock. PhD thesis, Melbourne, Australia: Monash University; 1989.

M. Bahaaddini et al. / Computers and Geotechnics 51 (2013) 101–115 [5] Seidel JP, Haberfield CM. The application of energy principles to the determination of the sliding resistance of rock joints. Rock Mech Rock Eng 1995;28(4):211–26. [6] Ladanyi B, Archambault G. Direct and indirect determination of shear strength of rock mass. In: Preprint number 80-25 AIME annual meeting. Las Vegas, Nevada; 1980. [7] Fairhurst C. On the validity of the ‘Brazilian’ test for brittle materials. Int J Rock Mech Min Sci Geomech Abstr 1964;1(4):535–46. [8] Hoek E, Bray J. Rock slope engineering. London: The Institution of Mining and Metallurgy; 1977. [9] Seidel JP. The analysis and design of pile shafts in weak rock. PhD thesis, Melbourne, Australia: Monash University; 1993. [10] Barton N, Choubey V. The shear strength of rock joints in theory and practice. Rock Mech 1977;10(1–2):1–54. [11] Olsson R, Barton N. An improved model for hydromechanical coupling during shearing of rock joints. Int J Rock Mech Mining Sci 2001;38(3):317–29. [12] Indraratna B, Haque A. Experimental and numerical modeling of shear behaviour of rock joints. In: GeoEng 2000, An international conference on geotechnical and geological engineering. Pennsylvania, USA; 2000. [13] Vosniakos K. Physical and numerical modelling of shear behaviour of sawtoothed filled rock joint. PhD thesis, Manchester: University of Manchester; 2007. [14] Oh JM. Three dimensional numerical modeling of excavation in rock with dilatant joints. PhD thesis, Urbana, Illinois: University of Illinois; 2005. [15] Karami A, Stead D. Asperity degradation and damage in the direct shear test: a hybrid FEM/DEM approach. Rock Mech Rock Eng 2008;41(2):229–66. [16] Hunt SP, Meyers AG, Louchnikov V. Modelling the Kaiser effect and deformation rate analysis in sandstone using the discrete element method. Comput Geotech 2003;30(7):611–21. [17] 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. [18] Tran TH, Vénier R, Cambou B. Discrete modelling of rock-ageing in rockfill dams. Comput Geotech 2009;36(1–2):264–75. [19] Hazzard JF, Young RP, Maxwell SC. Micromechanical modeling of cracking and failure in brittle rocks. J Geophys Res 2000;105(B7):16683–97. [20] Zhang ZX, Hu XY, Scott KD. A discrete numerical approach for modeling face stability in slurry shield tunnelling in soft soils. Comput Geotech 2011;38(1):94–104. [21] Akram MS, Sharrock GB. Physical and numerical investigation of a cemented granular assembly of steel spheres. Int J Numer Anal Methods Geomech 2010;34(18):1896–934. [22] Itasca Consulting Group Inc. PFC2D manual, version 4.0 Minneapolis, Minnesota; 2008. [23] Cho N, Martin CD, Sego DC. A clumped particle model for rock. Int J Rock Mech Min Sci 2007;44(7):997–1010. [24] Park JW, Song JJ. Numerical simulation of a direct shear test on a rock joint using a bonded-particle model. Int J Rock Mech Min Sci 2009;46(8):1315–28. [25] Asadi M, Rasouli V, Barla G. A bonded particle model simulation of shear strength and asperity degradation for rough rock fractures. Rock Mech Rock Eng 2012;45(5):649–75.

115

[26] Asadi MS. Experimental and PFC2D numerical study of progressive shear behaviour of single rough rock fractures. PhD thesis, Perth, Australia: Curtin University; 2011. [27] Cho N. Discrete element modeling of rock pre-peak fracturing and dilation. PhD thesis, Edmonton, Alberta: University of Alberta; 2008. [28] Sharrock GB, Akram MS, Mitra R. Application of synthetic rock mass modeling to estimate the strength of jointed sandstone. In: 43rd US Rock mechanics symposium and 4th US – Canada rock mechanics symposium. Asheville, North Carolina; 2009. [29] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41(8):1329–64. [30] Cundall PA. Numerical experiments on rough joints in shear using a bonded particle model. In: Aspects of tectonic faulting. Berlin; 2000. p. 1–9. [31] Kusumi H, Matsuoka T, Ashida Y, Tatsumi S. Simulation analysis of shear behavior of rock joint by distinct element method. In: Eurock 2005 – impact of human activity on geological environment. London; 2005. p. 281–6. [32] Jing L, Stephansson O. Fundamentals of discrete element methods for rock engineering: theory and applications. Amsterdam: Elsevier; 2007. [33] Lai ZHF, Douglas KJ, Mostyn G. The strength of rock defects – numerical analysis of scale effects. In: 11th congress of the international society for rock mechanics. Lisbon, Portugal; 2007. p. 481–4. [34] Rasouli V, Harrison JP. Assessment of rock fracture surface roughness using Riemannian statistics of linear profiles. Int J Rock Mech Min Sci 2010;47(6):940–8. [35] Bahaaddini M, Sharrock G, Hebbelwhite B, Mitra R. Direct shear tests to model the shear behaviour of rock joints by PFC2D. In: 46th US rock mechanics/ geomechanics symposium. Chicago, USA; 2012. [36] Pierce M, Cundall P, Potyondy D, Mas Ivars D. A synthetic rock mass model for jointed rock. In: Rock mechanics: meeting society’s challenges and demands. 1st Canada–US rock mechanics symposium. London: Vancouver; 2007. [37] Lambert C, Buzzi O, Giacomini A. Influence of calcium leaching on the mechanical behavior of a rock-mortar interface: a DEM analysis. Comput Geotech 2010;37(3):258–66. [38] Lambert C, Coll C. A DEM approach to rock joint strength estimate. In: Rock slope stability 2009. Santiago, Chile; 2009. [39] Mas Ivars D, Pierce ME, Darcel C, Reyes-Montes J, Potyondy DO, Young RP, et al. The synthetic rock mass approach for jointed rock mass modelling. Int J Rock Mech Min Sci 2011;48(2):219–44. [40] Bahaaddini M, Sharrock G, Hebblewhite BK. Numerical investigation of the effect of joint geometrical parameters on the mechanical properties of a nonpersistent jointed rock mass under uniaxial compression. Comput Geotech 2013;49:206–25. [41] Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci Geomech Abstr 1983;20(6):249–68. [42] Mogi K. Pressure dependence of rock strength and transition from brittle fracture to ductile flow. Bull Earthq Res Inst 1966;44:215–32. [43] Goodman RE. Methods of geological engineering in discontinuous rocks. West Publishing Co.; 1976.