- Email: [email protected]
Projected area-based strength estimation for jointed rock masses in triaxial compression
Computers and Geotechnics 104 (2018) 216–225
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
Projected area-based strength estimation for jointed rock masses in triaxial compression Zhenyu Hana, Dion Weatherleyb, Ruslan Puscasua,c,d,
T
⁎
a
WH Bryan Mining and Geology Research Centre, Sustainable Minerals Institute, The University of Queensland, St Lucia, QLD 4072, Australia Julius Kruttschnitt Mineral Research Centre, Sustainable Minerals Institute, The University of Queensland, Indooroopilly, QLD 4068, Australia c Mining3, Pinjarra Hills, Qld 4069, Australia d CSIRO Mineral Resources, Queensland Centre for Advanced Technologies, Brisbane, Australia b
A R T I C LE I N FO
A B S T R A C T
Keywords: Compressive strength Triaxial compression Discrete Element Method Fracture mechanics Rock mechanics
Reliable evaluation of the strength of rock masses is required for failure analysis in rock engineering. This paper describes the failure of rock mass specimens under triaxial compression via numerical testing based on the Discrete Element Method (DEM). It investigates the fracturing phenomena in rock masses with pre-existing joint sets and outlines a simple yet practical method for estimating compressive strength under complex joint geometric and loading conditions. The simulations capture different failure regimes, including intact, sliding and orthogonal failure, as joint geometric parameters are varied. Sensitivity studies demonstrate that the joint orientation, joint radius and continuity factor are essential geometrical terms affecting the material strength. The results reveal a linear relationship between the projected area of joint transections – which is introduced as a proxy to represent the three geometrical parameters – and the vertical strength. Based on this finding, a hypothesis for failure of jointed brittle materials is proposed, prescribing the influence of two factors: the complexity of the joint configuration and the spacing between the sample surface and the joint network. Such an approach, if validated, provides practitioners with a simple method for rapid estimation of compressive strength of rock-like materials via measurements of the joint geometry.
1. Introduction The main aim of fracture mechanics is to investigate material resistance to fractures. Fractures within rock masses dominate failure processes and the mechanical properties of rock. The lower strength of rock masses as compared with intact rock is caused by the weaker component – joints, the most common discontinuities – transecting the rock into pieces. Understanding the fracture mechanics of rock masses is essential to the design and performance prediction of constructions built in and on rock masses. It has significant effects on mining [1], mineral processing [2,3], civil engineering [4–6] and many other disciplines of science and engineering [7,8]. Although the evidence for complex behaviours of rock masses under loading is overwhelming, the dynamics of rock masses under triaxial compression is still obscure. While it has long been known that fracture networks produce complexity, the propagation of fractures and their influence on the rock strength remain poorly understood. The complex arrangement of joints makes it extremely difficult to determine the mechanical properties of jointed materials [9]. The
⁎
mechanical response of rock masses under different applied loads has different sensitivities to various geometrical parameters such as joint persistence, orientation, continuity, etc. [10]. Due to the large number of potential interactions among these parameters, it can prove difficult to quantitatively determine their effect on rock failure and strength. Applications of current failure criteria in rock strength estimation are not convenient. Usually, practitioners need to identify the discontinuity types using a particular system of rock mass classification and determine the parameter values required by the criterion. During this process, manual errors often arise in rock mass classification, often caused by subjective definitions or a wide value-range for each joint type. Moreover, specific failure criteria can only be used for particular forms of rock conditions from which they were deduced; sometimes being misused beyond their original intended application range [11,12]. Estimating material strength based on the discontinuity geometry enjoys many benefits. Firstly, this approach is objective because the strength assessment is based on measurable properties, which eliminates the manual errors associated with subjective rock mass
Corresponding author at: CSIRO Mineral Resources, Queensland Centre for Advanced Technologies, Brisbane, Australia. E-mail address: [email protected] (R. Puscasu).
https://doi.org/10.1016/j.compgeo.2018.08.020 Received 5 January 2018; Received in revised form 12 August 2018; Accepted 28 August 2018 0266-352X/ © 2018 Elsevier Ltd. All rights reserved.
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
2.2. Particle interactions and loading mechanisms
classification. Secondly, the geometry-based strength estimation could be applied to a range of jointed materials; eliminating the need to deduce specific failure criteria for certain materials. Much effort has been made to reveal the mechanical properties of jointed rocks via laboratory experimentation and numerical simulation. Effects of persistent joints on the material properties have been widely studied in the literature [13–17]. However, understanding of the behaviour of non-persistent jointed materials is limited. Although, some scholars demonstrate the influence of the joint geometry on material features [18–21,10], most specimens used in these studies are thin slices which lead to results similar to those obtained in a 2D environment. Although laboratory experiments have been broadly used to understand the complex mechanical features of rock masses, they are difficult to assemble and conduct. Preparing specimens with artificial joints is challenging. Furthermore, the direct observation and measurement of fracturing phenomena within most materials is laborious due to the rapidity of crack propagation and physical opaqueness. Also, stress field perturbations caused by sample grips or boundaries cannot be ignored. Equipment use, experiment configurations and sample preparation are also relatively costly. In all, current practical impediments limit data acquisition for analysis of the failure behaviour of jointed rock masses. Numerical modelling, on the other hand, overcomes the shortcomings of laboratory experiments and provides an alternative avenue to reveal the insights into the factors governing the strength and fracturing response of jointed rock masses. Numerical methods permit direct observation and measurement of fracturing phenomena with systematically varying parameters. In addition, simulations provide adequate flexibility in dealing with complex material features involving inhomogeneity, anisotropy and boundary conditions [22]. Moreover, the cost of conducting numerical studies is negligible compared with corresponding laboratory tests. Herein, a bonded particle model is employed to investigate deformation processes in jointed material under triaxial compression. The DEM enables simulation of unconstrained nucleation and propagation of cracks and can model complicated material features such as heterogeneities and voids without prohibitive increases in computational complexity [23].
The interactions between two adjacent DEM particles involve six degrees of freedoms – tension, compression, shearing, torsion and bending – which can be expressed by four Hookean elastic interactions; under the isotropic assumption for shearing and bending deformation [26]. The failure of such elastic beam interactions is governed by a generalized Mohr-Coulomb failure criterion:
σs ⩾ c + σN tanθ ,
(1)
where σs and σN are the shear and normal stress within the elastic beam; c and θ are respectively the cohesion and the friction angle of the interaction. The bond adjoining particles i and j breaks forming a macroscopic fracture surface. The normal stress (σN ; positive under tension) and shear stress (σs ) are computed in each time step as follows
σN =
FN |MB | Rij , + A I
(2)
σs =
|FS | |MT | Rij , + A J
(3)
where FN and FS mean normal force and shear force between the particle i and j, MB and MT represent bending moment and torsion moment, A, I and J denote cross-sectional area, moment of inertia and polar moment of inertia of the interaction, and Rij is the half length between the two centers of the adjacent particles i and j. Six rigid platens are placed and bonded to the six surfaces of cubic samples to apply boundary forces. Compression or tension within the DEM specimen is generated by platen movement towards or away from the specimens respectively. For simulations of triaxial compression tests, six platens initially accelerate uniformly to a prescribed stress. The four lateral platens thereafter maintain the compression stress at a prescribed value via a servo control mechanism. The two platens atop and below the specimen continue the compression at a constant rate until the specimen fails macroscopically. 2.3. Model description 2.3.1. Intact rock A rock mass usually comprises intact rock pieces transected by discontinuity planes. Via a geometrical space-filling sphere packing algorithm, around 110,000 non-overlapping spheres are packed into a 30 mm × 30 mm × 30 mm cubic region with beveled edges, as shown in Fig. 1. The radii of spheres range between rmin = 0.2 mm and rmax = 0.6 mm . Microscopic model parameters for this specimen are listed in Table 1. The uniaxial compressive strength of the intact specimen is σU = 127.71 MPa , measured via simulations with zero lateral stress. For the simulation of triaxial compression, the samples are compressed vertically and laterally by six movable platens. During the simulations, the six platens are initially moved with a stress increasing linearly for 20,000 time steps. After that the four lateral platens are fixed at a stress of 22.68 MPa and the top and below platens vertically compress specimens with a constant velocity.
2. Numerical method and model arrangements 2.1. Discrete element method The DEM is a numerical technique used to simulate the behaviour of composite materials formed mainly of granular components, wherein the specimen macroscopic behaviour is depicted as an assembly of microscopic motions of individual particles [24,25]. The force interactions on and between discrete elements depend upon the particular physical scenario to be simulated. Simulations herein were performed using ESyS-Particle (https://launchpad.net/esys-particle), a generalpurpose and open-source DEM software package. ESyS-Particle represents the solid material as an assembly of discrete spherical particles linked to adjacent particles with bonds that are spring-like connections. The macroscopic properties of the DEM specimen are governed by the micromechanical parameters regulating the strength and elasticity of bonds [26]. When the specimen deforms under external loading, individual bonded interactions break when the accumulated stress within the interaction exceeds a limitation governed by a prescribed failure criterion. These interactions, subsequently, are replaced with repulsive frictional interactions, simulating crack onset and subsequent frictional sliding. This approach has proven to be suitable for modelling deformation processes in brittle discontinuous materials subjected to a variety of external loading conditions.
2.3.2. Joint configuration A conceptual model is used to represent a non-ubiquitous rock mass with a single joint set distributed fracture network (DFN) [27,28]. All joints are penny-shaped and non-persistent. The geometrical parameters of the joint system follow the definition by Prudencio and Van Sint Jan [21]. β is the joint orientation relative to the σ1-direction; γ is the joint step angle; α is the joint tip to tip angle; r is the radius of the penny-shaped joint; Lr is the length of the rock bridge; d is the spacing between joint layers. Additionally, Lb is the boundary spacing between the specimen surface and the joint network and nj is the number of joint layers. The continuity factor k represents the ratio of the joint 217
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 1. Three dimensional geometrical configuration for intact rock with principal stresses. The particles near the edges of the top and bottom boundaries are normally bonded with two or three platens, which leads to very high stress concentrations in those regions. Therefore, the original sides of the cubic sample are removed to avoid the early nucleation of fractures along specimen edged corners.
random joints with Lb = 0 mm (Fig. 3(a) and (b)) and random joints with Lb = 4 mm (Fig. 3(c) and (d)). Geometrical parameters of joint distributions are chosen based on Baecher’s theory [29]. The values of joint parameters including location, radius, and orientation are obtained using a Monte-Carlo method. The parameter values of each joint are assigned independently and are blind to those of other joints. Specifically, the distribution of random joints is generated using the following considerations: the positions of joint centres form a Poisson distribution; the radius of joints follow an exponential distribution; the joint orientations follow a normal distribution with variance 0 or 5; the joint center, radius and orientation distributions are uncorrelated.
Table 1 The microscopic parameters of DEM numerical model. Item
Value
Unit
Constant parameters Microscopic Young’s modulus (E) Poisson’s ratio (ν ) Density ( ρ )
1.00 × 105 0.25
MPa −
3.00 × 10−3 0.6 0.4 0.20 0.60
g/mm3 − −
Static frictional coefficient ( μs ) Dynamic frictional coefficient ( μd ) Minimum particle radius (rmin ) Maximum particle radius (rmax ) Timestep increment (Δt ) Variable parameters Cohesion (c ) Friction angle (θ )
2.31 × 10−9
mm mm s
[10, 220] [5.7, 84.2]
MPa °
2.4. Test of scale invariance A distinct advantage of DEM rock specimens is the inherent scale independence of their mechanical properties. Whilst rock mechanics practitioners often discuss a representative element volume for rock masses, such a concept only applies to ubiquitously jointed rock masses. However, in cases of non-ubiquitous jointing such a representative volume cannot be defined. Since joints span many orders of magnitude in spatial scale and have self-similar geometry, no “characteristic scale” exists that could be identified as defining the minimum dimensions of a representative volume. At any given scale of observation, there will always be one or a few large joints subtending that volume. Such
persistence/ (joint persistence + rock bridge), which is introduced to investigate the combined effect of joint persistence and rock bridge on rock mass properties defined as follow
k=
2r . 2r + Lr
(4)
Four kinds of DFN were constructed including parallel joints with Lb = 0 mm (Fig. 2(a)), parallel joints with Lb = 4 mm (Fig. 2(b)),
Fig. 2. Joint geometrical parameters for defining a DFN comprised of parallel joints. 218
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 3. Schematic diagrams for DFN: (a) whole view of a sample with parallel joints and Lb = 0 mm ; (b) right view of a sample with parallel joints and Lb = 0 mm ; (c) whole view of a sample with parallel joints and Lb = 4 mm ; (d) right view of a sample with parallel joints and Lb = 4 mm ; (e) whole view of a sample with random joints and Lb = 0 mm ; (f) right view of a sample with random joints and Lb = 0 mm ; (g) whole view of a sample with random joints and Lb = 4 mm ; (h) right view of a sample with random joints and Lb = 4 mm .
differing particle-scale heterogeneity of the specimens. These test results attest that the scale of the DEM specimens employed in this study chosen to suit laboratory testing does not materially alter the outcomes or implications expounded in the discussion. Likewise, the term “joint” as apposed to “crack” is used to specifically highlight the relevance of the results to both laboratory-scale specimens and rock masses.
structures dominate the mechanical response of the rock volume; their geometrical size and orientation relative to the principal stress directions being of first-order importance. For a geotechnical engineer wishing to estimate the failure strength of a bench or tunnel drift, the mechanical response of a so-called representative element volume is of secondary importance in the presence of joints whose dimensions approach that of the volume under consideration. Since DEM specimens can be scaled to suit specific size and scope, we first conduct three sets of simulations of differing spatial scale to test the scale invariance of mechanical response. Each set of simulations comprises four specimens with similar particle numbers and the same joint pattern, albeit of differing scales ranging from 30 mm to 3 m. The loading conditions (principal stresses and strain rates) are the same for all specimens considered. The detailed configurations for the nine simulations are shown in Table 2. Plots of specimen jointing patterns (in trans-section) and stress-strain response curves are also provided in Figs. 4–6. It is evident that the strain-stress curves for specimens with same joint patterns are similar. The numerical differences in mechanical properties (Young’s modulus and peak stress) are attributable to the differing arrangements of DEM spheres within each specimen; the
3. Simulation results This section presents the numerical results of triaxial compression tests, whose aim is to elucidate the influence of joint geometrical parameters on the failure mechanisms and strength of rock masses. Firstly, a sensitivity study is conducted wherein 215 simulations of triaxial compression tests are performed on a specimen with parallel joints or a DFN. The parallel joints are constructed for 2 values of the joint radius in the range [3, 5] mm , 19 different values of the joint orientation in the range [0°, 90°], 6 values of the continuity factor in the range [0.50, 0.97], 2 values of the spacing (Lb ) between the specimen surface and the joint area [0 mm, 4 mm], and 2 values of the spacing between adjacent joint layers [3 mm, 7 mm]. A DFN is constructed with 20 pre-existing joints for 20 random values of the joint position (average is 15 mm), 20 random values of joint orientations and 20 random values of the joint radius.
Table 2 Model description for scale invariance tests. Specimen
Joint orientation
Size (mm × mm × mm )
Particle number
Particle radius (mm)
Intact rock
–
30 × 30 × 30 300 × 300 × 300 3000 × 3000 × 3000
109115 111954 110565
[0.2, 0.6] [2, 6] [20, 60]
A rock mass with one joint layer
5°
30 × 30 × 30 300 × 300 × 300 3000 × 3000 × 3000
107892 110965 110142
[0.2, 0.6] [2, 6] [20, 60]
A rock mass with three joint layers
45°
30 × 30 × 30 300 × 300 × 300 3000 × 3000 × 3000
107019 110124 109757
[0.2, 0.6] [2, 6] [20, 60]
3.1. Failure modes The failure modes are identified based on visual observations of fracture locations after 1.0 × 105 simulation time-steps. It is observed that fractures tend to concentrate at the vertices of the cubic specimen and fracturing zones generally occur along the lines of the vertices. This phenomenon results in the appearance of conjugate or orthogonal fractures. Although there are complex and mixed modes, it is possible to identify three dominant fracture modes. Schematic sketches of the failure modes are presented in Fig. 7. The failure modes are:
• Intact failure: Initially, most fractures occur at corners and margins
of the top and bottom surfaces. Following the increase in the principal stress σ1, fractures are uniformly generated within the specimen volume.
219
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 4. Trans-sections of DEM models: (a) intact rock; (b) a rock mass with one joint layer; (c) a rock mass with three joint layers.
Fig. 5. Timestep-strain curves of DEM models for varying model scales: (a) intact rock; (b) a rock mass with one joint layer; (c) a rock mass with three joint layers.
Fig. 6. Stress-strain curves of DEM models for varying model scales: (a) intact rock; (b) a rock mass with one joint layer; (c) a rock mass with three joint layers.
Fig. 7. Schematic diagrams of fracture distribution in the middle part of specimens for varying joint orientations β .
• Sliding failure: Fractures firstly occur on the boundary of the top
• Orthogonal failure: In the initial stage of loading, fractures occur on
and bottom surface. As the vertical stress increases, fractures are generated around the edges of pre-existing joints. Thereafter fractures around the joints coalesce into cracks linking co-planar joints into a failure surface.
the boundary of the top and bottom surface. Subsequent propagation of fractures forms a fracturing zone linking the failure areas of the top and bottom. Note that the propagating direction of the fracturing zone in a direction perpendicular to the joint orientation.
220
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 8. Front (the first line) and right (the second line) views of fracture patterns observed in simulations of triaxial compression tests on discrete element method specimens with one joint layer. The number in parentheses indicate the joint orientation β .
Fig. 9. Front (the first line) and right (the second line) views of fracture patterns observed in simulations of triaxial compression tests on discrete element method specimens with three joint layers. Joint orientation is kept constant at 45° and the value of continuity factor varies from 0.50 to 0.91. The number in parentheses indicates the value of the continuity factor k.
vertical surfaces of specimens. In some cases (such as Fig. 8(g)), wingcracks and conjugate fractures are observed. For 40° ⩽ β ⩽ 60° sliding failure occurs: fractures are initiated at the top and bottom corners of specimens and subsequently coalesce and propagate along the joint layer. For 70° ⩽ β ⩽ 90° orthogonal failure occurs: after crack initiation, cracks coalesce and combine the pre-existing joint transections into a persistent failure surface whose orientation is perpendicular to the joint orientation. Based on a sensitivity study of the continuity factor k, which varies from 0.50 to 0.91, it is shown that failure patterns are very sensitive to k as shown in Fig. 9. At k = 0.50 , stress is distributed uniformly and the failure regime is dominated by intact failure. Due to the increase in k, stress concentrations move from the zone between parallel joints to joint tips. High stresses concentrate at joint tips, which subsequently drive tensile cracks from the joint tips. Thereafter, tensile cracks break rock bridges and grow through joint surfaces due to the shear stress. The parameter study of spacing (d) reveals that joint spacing
These modes can be observed for samples with various combinations of joint geometrical parameters. The following numerical results are collected to study the effect of these various parameters on the failure mode of the specimen. 3.2. Parameter study for failure regimes A sensitivity study of joint orientation β is conducted to qualitatively investigate its effect on failure regimes. As β increases from 0° to 90°, the failure pattern changes. Two joint patterns, one with one joint layer and the other with three joint layers, were used in this sensitivity study. Since the failure regimes are similar for both joint patterns, for simplicity only results with one joint layer are discussed herein. Numerical tests demonstrate that deformation patterns strongly depend on joint orientation as shown in Fig. 8. For 0° ⩽ β ⩽ 30° intact fracture occurs: fractures initially appear around the vertices, nucleate around the direction of pre-existing joints, but seldom occur along the 221
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 10. Front (the first line) and right (the second line) views of fracture patterns observed in simulations of triaxial compression tests on discrete element method specimens with three joint layers. Joint orientation is kept constant at 45° and the value of spacing d varies from 3 mm to 7 mm. The number in parentheses indicates the value of spacing d.
Fig. 11. Strength ratio Rs of specimens with one or three joint layers as a function of joint orientation β for different values of the boundary spacing Lb (mm) and the continuity factor k.
the increase in β , the ratio Rs decreases gradually and reaches its minimum value around β = 60°. Thereafter, the ratio steadily increases until β reaches 90°. To investigate its effect on vertical strength, the value of k is increased from 0.50 to 0.91. The resulting relationship between the ratio Rs and continuity factor k are plotted in Fig. 12. As can be seen in these two diagrams, the values of the ratio Rs display a decreasing trend as the continuity factor k increases. Based on Eq. (4), as the continuity factor k increases, radii of joints increase and the bridge length between joints decreases. Based on the fracture mechanics theory, the stress intensity factors increase due to the increase in the joint length, which leads to the decrease in the rock mass strength. Strength sensitivities to the continuity factor k are varied for different values of joint orientation (β ). The larger the joint orientation is, the more sensitive rock mass strength is to the continuity factor k. Numerical results indicate that not only is strength sensitive to continuity factor but it also strongly depends on joint radius r (Fig. 13). When the value of r increases from 3 mm to 5 mm, Rs appears to be more sensitive to β and k. As r increases Rs decreases, or in other words,
influences the failure mode somewhat but the effect is not as pronounced as it is for β and k (see Fig. 10). For values of d in the range [3, 7] mm , the less the spacing is, the more fractures occur under the same triaxial compression conditions. As spacing increases, stress concentration migrates from the zone between parallel joint layers to the tips of the middle joint layer; fractures tend to be distributed from the area along the parallel joint layers to the whole interior of specimens.
3.3. Vertical strength For generalization, a ratio Rs is introduced to represent σ1/ σU where σ1 and σU are the vertical strength and uniaxial compression strength, respectively. The analysis of numerical results shows that the ratio Rs is sensitive to β , as seen in Fig. 11. By varying β in the range [0°, 90°], Fig. 11 indicates that not only is the strength highly sensitive to joint orientation β but the curves of the ratio Rs to β hold similar trends for different values of Lb and k. At β = 0°, the values of the ratio Rs for different curves are normally maximum in the entire range, evincing that vertical joints influence the sample strength the least. Following 222
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 12. Strength ratio Rs of specimens with three joint layers as a function of continuity factor k for different values of joint orientation β and boundary spacing Lb (mm).
Fig. 14. Strength ratio Rs of specimens with three joint layers as a function of joint orientation β for constant joint radius r = 5 mm and continuity factor k = 0.67 , and for different values of spacing d (mm).
increases from 3 mm to 7 mm . Therefore, the rock mass strength is sensitive to joint spacing when joints are not vertical or horizontal and generally increases as joint spacing widens.
Fig. 13. Strength ratio Rs of specimens with three joint layers as a function of continuity factor k for different values of joint radius r and joint orientation β . Solid and dashed lines represent the tests in which r = 3 mm and 5 mm, respectively.
4. Discussion The results show that the geometric configuration of joints does influence rock mass strength in a rather complex manner. However, the question remains: How to explicitly relate the geometrical configurations to strength? For simplicity, let us consider an alternative approach whereby a new geometrical term, the projected area sp of joint transections (shown in Fig. 15), is introduced as the total projected area of joint transections in the perpendicular plane of the major principal stress, which can be calculated by the following equation
larger joints decrease the strength of the jointed rock. This result may be explained from fracture mechanics principles, since stress intensity factors increase with crack length. The results shown in Fig. 14 indicate that there is little sample strength dependency on joint spacing d when β is equal to 0° and 90°. In contrast, the strength ratio Rs is more strongly dependent on d, when β is 30°, 45° or 60°. Moreover, Fig. 14 shows that a linear increase of d yields different dependencies of Rs as a function of β . At β = 90°, the Rs variation is negligible suggesting that the strength of jointed rock is not sensitive to spacing d for joints perpendicular to σ1. At β = 0°, the trends in strength ratio Rs are mixed with both weakening and hardening as d
nj
sp =
∑ j=0
223
sj sinβj ,
(5)
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
It is observed that for similar joint distributions the relationship between the strength ratio Rs = σ1/ σU and the area ratio Ra = sp/ sa is linear. Note that sa represents the lateral cross-sectional area of the specimen, as depicted in Fig. 15, and its value is 30 × 30 = 900 mm2 . A general linear fit to the data yields the following equation
Rs = Ksa·Ra + R 0 ,
(6)
where Ksa is the slope of the fitting equation which is a negative constant and R 0 is the vertical axial intercept of the fitting equation. The fit coefficients for each joint distribution are included in Fig. 16. It can be observed in Fig. 16 that more complex joint geometries generally lead to an increase in Ksa , from Ksa = −0.16 in the case of one joint layer (E1) to Ksa = −1.04 in the case of random joints (E5). Also, it is seen that the linear relationship is dependent on the boundary spacing Lb . Specifically, the larger Lb is, the smaller the slope of the fitting equations. For Ra ≈ 1.0 , the fitting equations including E2, E3 and E5 intersect, within a rectangular area depicted in the figure. Note that the fitting points of E4 on the left side of Ra = 1.0 also tend to cross this region.
Fig. 15. Illustration of the projected area sp and the cross section sa .
5. Conclusion
where βj is the joint orientation of the j-th joint and nj is the total number of joints. sj is the area of the j-th joint, which equals to πr j2 (r j is the radius of j-th joint) when the joint is circular. The quantitive analysis of the relationship between sp and σ1 reveals a monotonic almost linear dependency between the projected area and peak strength. To shed some light on this finding, a series of numerical true triaxial compression experiments are conducted for various geometric joint configurations. 307 triaxial compression test results (excluding the cases where early failure occurs) are represented in Fig. 16.
This paper is concerned with elucidating the conditions under which DEM-based rock mass specimens with pre-existing joints fail under triaxial compression. Simulations of samples containing both parallel and randomly distributed joints are conducted to capture various failure regimes including intact failure, sliding failure and orthogonal failure. A sensitivity study in which the joint geometries are systematically changed reveals that joint direction, joint radius, and the continuity factor are essential parameters influencing the specimen
Fig. 16. Strength ratio Rs as a function of area ratio Ra . r 2 denotes the goodness of fit. 224
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
strength. A projected area-based approach to estimate the compressive strength of jointed materials is also introduced. The rock mass compressive strength is found to be linearly related to the projected area of joint transections. Moreover, it shown that not only is strength a function of joint network configuration but also a function of joint boundary spacing or joint persistence. This finding, if validated, provides a simple and practical method to estimate rock mass strength requiring only knowledge of the joint geometry and intact rock strength (UCS). It enables material engineering practitioners to rapidly estimate the strength of materials with similar discontinuity networks subjected to triaxial compression, as it is relatively easy to determine the strength slope with only a few trials.
[11] Bieniawski Z. Classification of rock masses for engineering: the rmr system and future trends; 1993. p. 553–73. [12] Bieniawski ZT. Engineering rock mass classifications: a complete manual for engineers and geologists in mining, civil, and petroleum engineering. John Wiley & Sons; 1989. [13] Cai M, Kaiser P, Uno H, Tasaka Y, Minami M. Estimation of rock mass deformation modulus and strength of jointed hard rock masses using the {GSI} system. Int J Rock Mech Min Sci 2004;41:3–19. [14] Singh M, Rao SK, Ramamurthy T. Strength and deformational behaviour of a jointed rock mass. Rock Mech Rock Eng 2002;35:45–64. [15] Singh M, Rao KS. Empirical methods to estimate the strength of jointed rock masses. Eng Geol 2005;77:127–37. [16] Wang T-T, Huang T-H. A constitutive model for the deformation of a rock mass containing sets of ubiquitous joints. Int J Rock Mech Min Sci 2009;46:521–30. [17] Zhao W, Huangb R, Yanc M. Mechanical and fracture behavior of rock mass with parallel concentrated joints with different dip angle and number based on pfc simulation. Geomech Eng 2015;8:757–67. [18] Niu S, Jing H, Hu K, Yang D. Numerical investigation on the sensitivity of jointed rock mass strength to various factors. Min Sci Technol (China) 2010;20:530–4. [19] Chen X, Liao Z, Peng X. Deformability characteristics of jointed rock masses under uniaxial compression. Int J Min Sci Technol 2012;22:213–21. [20] Fan X, Kulatilake P, Chen X. Mechanical behavior of rock-like jointed blocks with multi-non-persistent joints under uniaxial loading: a particle mechanics approach. Eng Geol 2015;190:17–32. [21] Prudencio M, Jan MVS. Strength and failure modes of rock mass models with nonpersistent joints. Int J Rock Mech Min Sci 2007;44:890–902. [22] Jing L. A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. Int J Rock Mech Min Sci 2003;40:283–353. [23] Wang Y, Mora P. The esys_particle: a new 3-d discrete element model with single particle rotation. Advances in geocomputing. Springer; 2009. p. 183–228. [24] Cundall PA, Strack OD. A discrete numerical model for granular assemblies. Geotechnique 1979;29:47–65. [25] Hancock W, Weatherley D, Chitombo G. Modeling the gravity flow of rock using the discrete element method. In: MassMin 2012: 6th international conference and exhibition on mass mining, Canadian Institute of Mining, Metallurgy and Petroleum. [26] Wang Y, Mora P. Macroscopic elastic properties of regular lattices. J Mech Phys Solids 2008;56:3459–74. [27] Cacas MC, Ledoux E, de Marsily G, Tillie B, Barbreau A, Durand E, et al. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model. Water Resour Res 1990;26:479–89. [28] Cacas MC, Ledoux E, de Marsily G, Barbreau A, Calmels P, Gaillard B, et al. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 2. The transport model. Water Resour Res 1990;26:491–500. [29] Baecher G, Lanney N, Einstein H, et al. Statistical description of rock properties and sampling. In: The 18th US symposium on rock mechanics (USRMS), American Rock Mechanics Association.
References [1] Singh RN, Sun G. Applications of fracture mechanics to some mining engineering problems. Min Sci Technol 1990;10:53–60. [2] Atkins A. Modelling metal cutting using modern ductile fracture mechanics: quantitative explanations for some longstanding problems. Int J Mech Sci 2003;45:373–96. [3] Floreen S. The creep fracture of wrought nickel-base alloys by a fracture mechanics approach. Metall Trans A 1975;6:1741–9. [4] Bažant ZP, Cedolin L. Fracture mechanics of reinforced concrete. J Eng Mech Divis 1980;106:1287–306. [5] Han Z, Weatherley D, Puscasu R. Application of discrete element method to model crack propagation. In: ISRM Congress 2015 proceedings. [6] Han Z, Weatherley D, Puscasu R. A relationship between tensile strength and loading stress governing the onset of mode i crack propagation obtained via numerical investigations using a bonded particle model. Int J Numer Anal Methods Geomech 2017:1–13. [7] Rossmanith HP, Daehnke A, Nasmillner REK, Kouzniak N, Ohtsu M, Uenishi K. Fracture mechanics applications to drilling and blasting. Fatigue Fract Eng Mater Struct 1997;20:1617–36. [8] Gao H. Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials. Int J Fract 2006;138:101–37. [9] Kulatilake P, Ucpirti H, Wang S, Radberg G, Stephansson O. Use of the distinct element method to perform stress analysis in rock with non-persistent joints and to study the effect of joint geometry parameters on the strength and deformability of rock masses. Rock Mech Rock Eng 1992;25:253–74. [10] Bahaaddini M, Sharrock G, Hebblewhite B. Numerical investigation of the effect of joint geometrical parameters on the mechanical properties of a non-persistent jointed rock mass under uniaxial compression. Comput Geotech 2013;49:206–25.
225
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
Projected area-based strength estimation for jointed rock masses in triaxial compression Zhenyu Hana, Dion Weatherleyb, Ruslan Puscasua,c,d,
T
⁎
a
WH Bryan Mining and Geology Research Centre, Sustainable Minerals Institute, The University of Queensland, St Lucia, QLD 4072, Australia Julius Kruttschnitt Mineral Research Centre, Sustainable Minerals Institute, The University of Queensland, Indooroopilly, QLD 4068, Australia c Mining3, Pinjarra Hills, Qld 4069, Australia d CSIRO Mineral Resources, Queensland Centre for Advanced Technologies, Brisbane, Australia b
A R T I C LE I N FO
A B S T R A C T
Keywords: Compressive strength Triaxial compression Discrete Element Method Fracture mechanics Rock mechanics
Reliable evaluation of the strength of rock masses is required for failure analysis in rock engineering. This paper describes the failure of rock mass specimens under triaxial compression via numerical testing based on the Discrete Element Method (DEM). It investigates the fracturing phenomena in rock masses with pre-existing joint sets and outlines a simple yet practical method for estimating compressive strength under complex joint geometric and loading conditions. The simulations capture different failure regimes, including intact, sliding and orthogonal failure, as joint geometric parameters are varied. Sensitivity studies demonstrate that the joint orientation, joint radius and continuity factor are essential geometrical terms affecting the material strength. The results reveal a linear relationship between the projected area of joint transections – which is introduced as a proxy to represent the three geometrical parameters – and the vertical strength. Based on this finding, a hypothesis for failure of jointed brittle materials is proposed, prescribing the influence of two factors: the complexity of the joint configuration and the spacing between the sample surface and the joint network. Such an approach, if validated, provides practitioners with a simple method for rapid estimation of compressive strength of rock-like materials via measurements of the joint geometry.
1. Introduction The main aim of fracture mechanics is to investigate material resistance to fractures. Fractures within rock masses dominate failure processes and the mechanical properties of rock. The lower strength of rock masses as compared with intact rock is caused by the weaker component – joints, the most common discontinuities – transecting the rock into pieces. Understanding the fracture mechanics of rock masses is essential to the design and performance prediction of constructions built in and on rock masses. It has significant effects on mining [1], mineral processing [2,3], civil engineering [4–6] and many other disciplines of science and engineering [7,8]. Although the evidence for complex behaviours of rock masses under loading is overwhelming, the dynamics of rock masses under triaxial compression is still obscure. While it has long been known that fracture networks produce complexity, the propagation of fractures and their influence on the rock strength remain poorly understood. The complex arrangement of joints makes it extremely difficult to determine the mechanical properties of jointed materials [9]. The
⁎
mechanical response of rock masses under different applied loads has different sensitivities to various geometrical parameters such as joint persistence, orientation, continuity, etc. [10]. Due to the large number of potential interactions among these parameters, it can prove difficult to quantitatively determine their effect on rock failure and strength. Applications of current failure criteria in rock strength estimation are not convenient. Usually, practitioners need to identify the discontinuity types using a particular system of rock mass classification and determine the parameter values required by the criterion. During this process, manual errors often arise in rock mass classification, often caused by subjective definitions or a wide value-range for each joint type. Moreover, specific failure criteria can only be used for particular forms of rock conditions from which they were deduced; sometimes being misused beyond their original intended application range [11,12]. Estimating material strength based on the discontinuity geometry enjoys many benefits. Firstly, this approach is objective because the strength assessment is based on measurable properties, which eliminates the manual errors associated with subjective rock mass
Corresponding author at: CSIRO Mineral Resources, Queensland Centre for Advanced Technologies, Brisbane, Australia. E-mail address: [email protected] (R. Puscasu).
https://doi.org/10.1016/j.compgeo.2018.08.020 Received 5 January 2018; Received in revised form 12 August 2018; Accepted 28 August 2018 0266-352X/ © 2018 Elsevier Ltd. All rights reserved.
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
2.2. Particle interactions and loading mechanisms
classification. Secondly, the geometry-based strength estimation could be applied to a range of jointed materials; eliminating the need to deduce specific failure criteria for certain materials. Much effort has been made to reveal the mechanical properties of jointed rocks via laboratory experimentation and numerical simulation. Effects of persistent joints on the material properties have been widely studied in the literature [13–17]. However, understanding of the behaviour of non-persistent jointed materials is limited. Although, some scholars demonstrate the influence of the joint geometry on material features [18–21,10], most specimens used in these studies are thin slices which lead to results similar to those obtained in a 2D environment. Although laboratory experiments have been broadly used to understand the complex mechanical features of rock masses, they are difficult to assemble and conduct. Preparing specimens with artificial joints is challenging. Furthermore, the direct observation and measurement of fracturing phenomena within most materials is laborious due to the rapidity of crack propagation and physical opaqueness. Also, stress field perturbations caused by sample grips or boundaries cannot be ignored. Equipment use, experiment configurations and sample preparation are also relatively costly. In all, current practical impediments limit data acquisition for analysis of the failure behaviour of jointed rock masses. Numerical modelling, on the other hand, overcomes the shortcomings of laboratory experiments and provides an alternative avenue to reveal the insights into the factors governing the strength and fracturing response of jointed rock masses. Numerical methods permit direct observation and measurement of fracturing phenomena with systematically varying parameters. In addition, simulations provide adequate flexibility in dealing with complex material features involving inhomogeneity, anisotropy and boundary conditions [22]. Moreover, the cost of conducting numerical studies is negligible compared with corresponding laboratory tests. Herein, a bonded particle model is employed to investigate deformation processes in jointed material under triaxial compression. The DEM enables simulation of unconstrained nucleation and propagation of cracks and can model complicated material features such as heterogeneities and voids without prohibitive increases in computational complexity [23].
The interactions between two adjacent DEM particles involve six degrees of freedoms – tension, compression, shearing, torsion and bending – which can be expressed by four Hookean elastic interactions; under the isotropic assumption for shearing and bending deformation [26]. The failure of such elastic beam interactions is governed by a generalized Mohr-Coulomb failure criterion:
σs ⩾ c + σN tanθ ,
(1)
where σs and σN are the shear and normal stress within the elastic beam; c and θ are respectively the cohesion and the friction angle of the interaction. The bond adjoining particles i and j breaks forming a macroscopic fracture surface. The normal stress (σN ; positive under tension) and shear stress (σs ) are computed in each time step as follows
σN =
FN |MB | Rij , + A I
(2)
σs =
|FS | |MT | Rij , + A J
(3)
where FN and FS mean normal force and shear force between the particle i and j, MB and MT represent bending moment and torsion moment, A, I and J denote cross-sectional area, moment of inertia and polar moment of inertia of the interaction, and Rij is the half length between the two centers of the adjacent particles i and j. Six rigid platens are placed and bonded to the six surfaces of cubic samples to apply boundary forces. Compression or tension within the DEM specimen is generated by platen movement towards or away from the specimens respectively. For simulations of triaxial compression tests, six platens initially accelerate uniformly to a prescribed stress. The four lateral platens thereafter maintain the compression stress at a prescribed value via a servo control mechanism. The two platens atop and below the specimen continue the compression at a constant rate until the specimen fails macroscopically. 2.3. Model description 2.3.1. Intact rock A rock mass usually comprises intact rock pieces transected by discontinuity planes. Via a geometrical space-filling sphere packing algorithm, around 110,000 non-overlapping spheres are packed into a 30 mm × 30 mm × 30 mm cubic region with beveled edges, as shown in Fig. 1. The radii of spheres range between rmin = 0.2 mm and rmax = 0.6 mm . Microscopic model parameters for this specimen are listed in Table 1. The uniaxial compressive strength of the intact specimen is σU = 127.71 MPa , measured via simulations with zero lateral stress. For the simulation of triaxial compression, the samples are compressed vertically and laterally by six movable platens. During the simulations, the six platens are initially moved with a stress increasing linearly for 20,000 time steps. After that the four lateral platens are fixed at a stress of 22.68 MPa and the top and below platens vertically compress specimens with a constant velocity.
2. Numerical method and model arrangements 2.1. Discrete element method The DEM is a numerical technique used to simulate the behaviour of composite materials formed mainly of granular components, wherein the specimen macroscopic behaviour is depicted as an assembly of microscopic motions of individual particles [24,25]. The force interactions on and between discrete elements depend upon the particular physical scenario to be simulated. Simulations herein were performed using ESyS-Particle (https://launchpad.net/esys-particle), a generalpurpose and open-source DEM software package. ESyS-Particle represents the solid material as an assembly of discrete spherical particles linked to adjacent particles with bonds that are spring-like connections. The macroscopic properties of the DEM specimen are governed by the micromechanical parameters regulating the strength and elasticity of bonds [26]. When the specimen deforms under external loading, individual bonded interactions break when the accumulated stress within the interaction exceeds a limitation governed by a prescribed failure criterion. These interactions, subsequently, are replaced with repulsive frictional interactions, simulating crack onset and subsequent frictional sliding. This approach has proven to be suitable for modelling deformation processes in brittle discontinuous materials subjected to a variety of external loading conditions.
2.3.2. Joint configuration A conceptual model is used to represent a non-ubiquitous rock mass with a single joint set distributed fracture network (DFN) [27,28]. All joints are penny-shaped and non-persistent. The geometrical parameters of the joint system follow the definition by Prudencio and Van Sint Jan [21]. β is the joint orientation relative to the σ1-direction; γ is the joint step angle; α is the joint tip to tip angle; r is the radius of the penny-shaped joint; Lr is the length of the rock bridge; d is the spacing between joint layers. Additionally, Lb is the boundary spacing between the specimen surface and the joint network and nj is the number of joint layers. The continuity factor k represents the ratio of the joint 217
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 1. Three dimensional geometrical configuration for intact rock with principal stresses. The particles near the edges of the top and bottom boundaries are normally bonded with two or three platens, which leads to very high stress concentrations in those regions. Therefore, the original sides of the cubic sample are removed to avoid the early nucleation of fractures along specimen edged corners.
random joints with Lb = 0 mm (Fig. 3(a) and (b)) and random joints with Lb = 4 mm (Fig. 3(c) and (d)). Geometrical parameters of joint distributions are chosen based on Baecher’s theory [29]. The values of joint parameters including location, radius, and orientation are obtained using a Monte-Carlo method. The parameter values of each joint are assigned independently and are blind to those of other joints. Specifically, the distribution of random joints is generated using the following considerations: the positions of joint centres form a Poisson distribution; the radius of joints follow an exponential distribution; the joint orientations follow a normal distribution with variance 0 or 5; the joint center, radius and orientation distributions are uncorrelated.
Table 1 The microscopic parameters of DEM numerical model. Item
Value
Unit
Constant parameters Microscopic Young’s modulus (E) Poisson’s ratio (ν ) Density ( ρ )
1.00 × 105 0.25
MPa −
3.00 × 10−3 0.6 0.4 0.20 0.60
g/mm3 − −
Static frictional coefficient ( μs ) Dynamic frictional coefficient ( μd ) Minimum particle radius (rmin ) Maximum particle radius (rmax ) Timestep increment (Δt ) Variable parameters Cohesion (c ) Friction angle (θ )
2.31 × 10−9
mm mm s
[10, 220] [5.7, 84.2]
MPa °
2.4. Test of scale invariance A distinct advantage of DEM rock specimens is the inherent scale independence of their mechanical properties. Whilst rock mechanics practitioners often discuss a representative element volume for rock masses, such a concept only applies to ubiquitously jointed rock masses. However, in cases of non-ubiquitous jointing such a representative volume cannot be defined. Since joints span many orders of magnitude in spatial scale and have self-similar geometry, no “characteristic scale” exists that could be identified as defining the minimum dimensions of a representative volume. At any given scale of observation, there will always be one or a few large joints subtending that volume. Such
persistence/ (joint persistence + rock bridge), which is introduced to investigate the combined effect of joint persistence and rock bridge on rock mass properties defined as follow
k=
2r . 2r + Lr
(4)
Four kinds of DFN were constructed including parallel joints with Lb = 0 mm (Fig. 2(a)), parallel joints with Lb = 4 mm (Fig. 2(b)),
Fig. 2. Joint geometrical parameters for defining a DFN comprised of parallel joints. 218
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 3. Schematic diagrams for DFN: (a) whole view of a sample with parallel joints and Lb = 0 mm ; (b) right view of a sample with parallel joints and Lb = 0 mm ; (c) whole view of a sample with parallel joints and Lb = 4 mm ; (d) right view of a sample with parallel joints and Lb = 4 mm ; (e) whole view of a sample with random joints and Lb = 0 mm ; (f) right view of a sample with random joints and Lb = 0 mm ; (g) whole view of a sample with random joints and Lb = 4 mm ; (h) right view of a sample with random joints and Lb = 4 mm .
differing particle-scale heterogeneity of the specimens. These test results attest that the scale of the DEM specimens employed in this study chosen to suit laboratory testing does not materially alter the outcomes or implications expounded in the discussion. Likewise, the term “joint” as apposed to “crack” is used to specifically highlight the relevance of the results to both laboratory-scale specimens and rock masses.
structures dominate the mechanical response of the rock volume; their geometrical size and orientation relative to the principal stress directions being of first-order importance. For a geotechnical engineer wishing to estimate the failure strength of a bench or tunnel drift, the mechanical response of a so-called representative element volume is of secondary importance in the presence of joints whose dimensions approach that of the volume under consideration. Since DEM specimens can be scaled to suit specific size and scope, we first conduct three sets of simulations of differing spatial scale to test the scale invariance of mechanical response. Each set of simulations comprises four specimens with similar particle numbers and the same joint pattern, albeit of differing scales ranging from 30 mm to 3 m. The loading conditions (principal stresses and strain rates) are the same for all specimens considered. The detailed configurations for the nine simulations are shown in Table 2. Plots of specimen jointing patterns (in trans-section) and stress-strain response curves are also provided in Figs. 4–6. It is evident that the strain-stress curves for specimens with same joint patterns are similar. The numerical differences in mechanical properties (Young’s modulus and peak stress) are attributable to the differing arrangements of DEM spheres within each specimen; the
3. Simulation results This section presents the numerical results of triaxial compression tests, whose aim is to elucidate the influence of joint geometrical parameters on the failure mechanisms and strength of rock masses. Firstly, a sensitivity study is conducted wherein 215 simulations of triaxial compression tests are performed on a specimen with parallel joints or a DFN. The parallel joints are constructed for 2 values of the joint radius in the range [3, 5] mm , 19 different values of the joint orientation in the range [0°, 90°], 6 values of the continuity factor in the range [0.50, 0.97], 2 values of the spacing (Lb ) between the specimen surface and the joint area [0 mm, 4 mm], and 2 values of the spacing between adjacent joint layers [3 mm, 7 mm]. A DFN is constructed with 20 pre-existing joints for 20 random values of the joint position (average is 15 mm), 20 random values of joint orientations and 20 random values of the joint radius.
Table 2 Model description for scale invariance tests. Specimen
Joint orientation
Size (mm × mm × mm )
Particle number
Particle radius (mm)
Intact rock
–
30 × 30 × 30 300 × 300 × 300 3000 × 3000 × 3000
109115 111954 110565
[0.2, 0.6] [2, 6] [20, 60]
A rock mass with one joint layer
5°
30 × 30 × 30 300 × 300 × 300 3000 × 3000 × 3000
107892 110965 110142
[0.2, 0.6] [2, 6] [20, 60]
A rock mass with three joint layers
45°
30 × 30 × 30 300 × 300 × 300 3000 × 3000 × 3000
107019 110124 109757
[0.2, 0.6] [2, 6] [20, 60]
3.1. Failure modes The failure modes are identified based on visual observations of fracture locations after 1.0 × 105 simulation time-steps. It is observed that fractures tend to concentrate at the vertices of the cubic specimen and fracturing zones generally occur along the lines of the vertices. This phenomenon results in the appearance of conjugate or orthogonal fractures. Although there are complex and mixed modes, it is possible to identify three dominant fracture modes. Schematic sketches of the failure modes are presented in Fig. 7. The failure modes are:
• Intact failure: Initially, most fractures occur at corners and margins
of the top and bottom surfaces. Following the increase in the principal stress σ1, fractures are uniformly generated within the specimen volume.
219
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 4. Trans-sections of DEM models: (a) intact rock; (b) a rock mass with one joint layer; (c) a rock mass with three joint layers.
Fig. 5. Timestep-strain curves of DEM models for varying model scales: (a) intact rock; (b) a rock mass with one joint layer; (c) a rock mass with three joint layers.
Fig. 6. Stress-strain curves of DEM models for varying model scales: (a) intact rock; (b) a rock mass with one joint layer; (c) a rock mass with three joint layers.
Fig. 7. Schematic diagrams of fracture distribution in the middle part of specimens for varying joint orientations β .
• Sliding failure: Fractures firstly occur on the boundary of the top
• Orthogonal failure: In the initial stage of loading, fractures occur on
and bottom surface. As the vertical stress increases, fractures are generated around the edges of pre-existing joints. Thereafter fractures around the joints coalesce into cracks linking co-planar joints into a failure surface.
the boundary of the top and bottom surface. Subsequent propagation of fractures forms a fracturing zone linking the failure areas of the top and bottom. Note that the propagating direction of the fracturing zone in a direction perpendicular to the joint orientation.
220
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 8. Front (the first line) and right (the second line) views of fracture patterns observed in simulations of triaxial compression tests on discrete element method specimens with one joint layer. The number in parentheses indicate the joint orientation β .
Fig. 9. Front (the first line) and right (the second line) views of fracture patterns observed in simulations of triaxial compression tests on discrete element method specimens with three joint layers. Joint orientation is kept constant at 45° and the value of continuity factor varies from 0.50 to 0.91. The number in parentheses indicates the value of the continuity factor k.
vertical surfaces of specimens. In some cases (such as Fig. 8(g)), wingcracks and conjugate fractures are observed. For 40° ⩽ β ⩽ 60° sliding failure occurs: fractures are initiated at the top and bottom corners of specimens and subsequently coalesce and propagate along the joint layer. For 70° ⩽ β ⩽ 90° orthogonal failure occurs: after crack initiation, cracks coalesce and combine the pre-existing joint transections into a persistent failure surface whose orientation is perpendicular to the joint orientation. Based on a sensitivity study of the continuity factor k, which varies from 0.50 to 0.91, it is shown that failure patterns are very sensitive to k as shown in Fig. 9. At k = 0.50 , stress is distributed uniformly and the failure regime is dominated by intact failure. Due to the increase in k, stress concentrations move from the zone between parallel joints to joint tips. High stresses concentrate at joint tips, which subsequently drive tensile cracks from the joint tips. Thereafter, tensile cracks break rock bridges and grow through joint surfaces due to the shear stress. The parameter study of spacing (d) reveals that joint spacing
These modes can be observed for samples with various combinations of joint geometrical parameters. The following numerical results are collected to study the effect of these various parameters on the failure mode of the specimen. 3.2. Parameter study for failure regimes A sensitivity study of joint orientation β is conducted to qualitatively investigate its effect on failure regimes. As β increases from 0° to 90°, the failure pattern changes. Two joint patterns, one with one joint layer and the other with three joint layers, were used in this sensitivity study. Since the failure regimes are similar for both joint patterns, for simplicity only results with one joint layer are discussed herein. Numerical tests demonstrate that deformation patterns strongly depend on joint orientation as shown in Fig. 8. For 0° ⩽ β ⩽ 30° intact fracture occurs: fractures initially appear around the vertices, nucleate around the direction of pre-existing joints, but seldom occur along the 221
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 10. Front (the first line) and right (the second line) views of fracture patterns observed in simulations of triaxial compression tests on discrete element method specimens with three joint layers. Joint orientation is kept constant at 45° and the value of spacing d varies from 3 mm to 7 mm. The number in parentheses indicates the value of spacing d.
Fig. 11. Strength ratio Rs of specimens with one or three joint layers as a function of joint orientation β for different values of the boundary spacing Lb (mm) and the continuity factor k.
the increase in β , the ratio Rs decreases gradually and reaches its minimum value around β = 60°. Thereafter, the ratio steadily increases until β reaches 90°. To investigate its effect on vertical strength, the value of k is increased from 0.50 to 0.91. The resulting relationship between the ratio Rs and continuity factor k are plotted in Fig. 12. As can be seen in these two diagrams, the values of the ratio Rs display a decreasing trend as the continuity factor k increases. Based on Eq. (4), as the continuity factor k increases, radii of joints increase and the bridge length between joints decreases. Based on the fracture mechanics theory, the stress intensity factors increase due to the increase in the joint length, which leads to the decrease in the rock mass strength. Strength sensitivities to the continuity factor k are varied for different values of joint orientation (β ). The larger the joint orientation is, the more sensitive rock mass strength is to the continuity factor k. Numerical results indicate that not only is strength sensitive to continuity factor but it also strongly depends on joint radius r (Fig. 13). When the value of r increases from 3 mm to 5 mm, Rs appears to be more sensitive to β and k. As r increases Rs decreases, or in other words,
influences the failure mode somewhat but the effect is not as pronounced as it is for β and k (see Fig. 10). For values of d in the range [3, 7] mm , the less the spacing is, the more fractures occur under the same triaxial compression conditions. As spacing increases, stress concentration migrates from the zone between parallel joint layers to the tips of the middle joint layer; fractures tend to be distributed from the area along the parallel joint layers to the whole interior of specimens.
3.3. Vertical strength For generalization, a ratio Rs is introduced to represent σ1/ σU where σ1 and σU are the vertical strength and uniaxial compression strength, respectively. The analysis of numerical results shows that the ratio Rs is sensitive to β , as seen in Fig. 11. By varying β in the range [0°, 90°], Fig. 11 indicates that not only is the strength highly sensitive to joint orientation β but the curves of the ratio Rs to β hold similar trends for different values of Lb and k. At β = 0°, the values of the ratio Rs for different curves are normally maximum in the entire range, evincing that vertical joints influence the sample strength the least. Following 222
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
Fig. 12. Strength ratio Rs of specimens with three joint layers as a function of continuity factor k for different values of joint orientation β and boundary spacing Lb (mm).
Fig. 14. Strength ratio Rs of specimens with three joint layers as a function of joint orientation β for constant joint radius r = 5 mm and continuity factor k = 0.67 , and for different values of spacing d (mm).
increases from 3 mm to 7 mm . Therefore, the rock mass strength is sensitive to joint spacing when joints are not vertical or horizontal and generally increases as joint spacing widens.
Fig. 13. Strength ratio Rs of specimens with three joint layers as a function of continuity factor k for different values of joint radius r and joint orientation β . Solid and dashed lines represent the tests in which r = 3 mm and 5 mm, respectively.
4. Discussion The results show that the geometric configuration of joints does influence rock mass strength in a rather complex manner. However, the question remains: How to explicitly relate the geometrical configurations to strength? For simplicity, let us consider an alternative approach whereby a new geometrical term, the projected area sp of joint transections (shown in Fig. 15), is introduced as the total projected area of joint transections in the perpendicular plane of the major principal stress, which can be calculated by the following equation
larger joints decrease the strength of the jointed rock. This result may be explained from fracture mechanics principles, since stress intensity factors increase with crack length. The results shown in Fig. 14 indicate that there is little sample strength dependency on joint spacing d when β is equal to 0° and 90°. In contrast, the strength ratio Rs is more strongly dependent on d, when β is 30°, 45° or 60°. Moreover, Fig. 14 shows that a linear increase of d yields different dependencies of Rs as a function of β . At β = 90°, the Rs variation is negligible suggesting that the strength of jointed rock is not sensitive to spacing d for joints perpendicular to σ1. At β = 0°, the trends in strength ratio Rs are mixed with both weakening and hardening as d
nj
sp =
∑ j=0
223
sj sinβj ,
(5)
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
It is observed that for similar joint distributions the relationship between the strength ratio Rs = σ1/ σU and the area ratio Ra = sp/ sa is linear. Note that sa represents the lateral cross-sectional area of the specimen, as depicted in Fig. 15, and its value is 30 × 30 = 900 mm2 . A general linear fit to the data yields the following equation
Rs = Ksa·Ra + R 0 ,
(6)
where Ksa is the slope of the fitting equation which is a negative constant and R 0 is the vertical axial intercept of the fitting equation. The fit coefficients for each joint distribution are included in Fig. 16. It can be observed in Fig. 16 that more complex joint geometries generally lead to an increase in Ksa , from Ksa = −0.16 in the case of one joint layer (E1) to Ksa = −1.04 in the case of random joints (E5). Also, it is seen that the linear relationship is dependent on the boundary spacing Lb . Specifically, the larger Lb is, the smaller the slope of the fitting equations. For Ra ≈ 1.0 , the fitting equations including E2, E3 and E5 intersect, within a rectangular area depicted in the figure. Note that the fitting points of E4 on the left side of Ra = 1.0 also tend to cross this region.
Fig. 15. Illustration of the projected area sp and the cross section sa .
5. Conclusion
where βj is the joint orientation of the j-th joint and nj is the total number of joints. sj is the area of the j-th joint, which equals to πr j2 (r j is the radius of j-th joint) when the joint is circular. The quantitive analysis of the relationship between sp and σ1 reveals a monotonic almost linear dependency between the projected area and peak strength. To shed some light on this finding, a series of numerical true triaxial compression experiments are conducted for various geometric joint configurations. 307 triaxial compression test results (excluding the cases where early failure occurs) are represented in Fig. 16.
This paper is concerned with elucidating the conditions under which DEM-based rock mass specimens with pre-existing joints fail under triaxial compression. Simulations of samples containing both parallel and randomly distributed joints are conducted to capture various failure regimes including intact failure, sliding failure and orthogonal failure. A sensitivity study in which the joint geometries are systematically changed reveals that joint direction, joint radius, and the continuity factor are essential parameters influencing the specimen
Fig. 16. Strength ratio Rs as a function of area ratio Ra . r 2 denotes the goodness of fit. 224
Computers and Geotechnics 104 (2018) 216–225
Z. Han et al.
strength. A projected area-based approach to estimate the compressive strength of jointed materials is also introduced. The rock mass compressive strength is found to be linearly related to the projected area of joint transections. Moreover, it shown that not only is strength a function of joint network configuration but also a function of joint boundary spacing or joint persistence. This finding, if validated, provides a simple and practical method to estimate rock mass strength requiring only knowledge of the joint geometry and intact rock strength (UCS). It enables material engineering practitioners to rapidly estimate the strength of materials with similar discontinuity networks subjected to triaxial compression, as it is relatively easy to determine the strength slope with only a few trials.
[11] Bieniawski Z. Classification of rock masses for engineering: the rmr system and future trends; 1993. p. 553–73. [12] Bieniawski ZT. Engineering rock mass classifications: a complete manual for engineers and geologists in mining, civil, and petroleum engineering. John Wiley & Sons; 1989. [13] Cai M, Kaiser P, Uno H, Tasaka Y, Minami M. Estimation of rock mass deformation modulus and strength of jointed hard rock masses using the {GSI} system. Int J Rock Mech Min Sci 2004;41:3–19. [14] Singh M, Rao SK, Ramamurthy T. Strength and deformational behaviour of a jointed rock mass. Rock Mech Rock Eng 2002;35:45–64. [15] Singh M, Rao KS. Empirical methods to estimate the strength of jointed rock masses. Eng Geol 2005;77:127–37. [16] Wang T-T, Huang T-H. A constitutive model for the deformation of a rock mass containing sets of ubiquitous joints. Int J Rock Mech Min Sci 2009;46:521–30. [17] Zhao W, Huangb R, Yanc M. Mechanical and fracture behavior of rock mass with parallel concentrated joints with different dip angle and number based on pfc simulation. Geomech Eng 2015;8:757–67. [18] Niu S, Jing H, Hu K, Yang D. Numerical investigation on the sensitivity of jointed rock mass strength to various factors. Min Sci Technol (China) 2010;20:530–4. [19] Chen X, Liao Z, Peng X. Deformability characteristics of jointed rock masses under uniaxial compression. Int J Min Sci Technol 2012;22:213–21. [20] Fan X, Kulatilake P, Chen X. Mechanical behavior of rock-like jointed blocks with multi-non-persistent joints under uniaxial loading: a particle mechanics approach. Eng Geol 2015;190:17–32. [21] Prudencio M, Jan MVS. Strength and failure modes of rock mass models with nonpersistent joints. Int J Rock Mech Min Sci 2007;44:890–902. [22] Jing L. A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. Int J Rock Mech Min Sci 2003;40:283–353. [23] Wang Y, Mora P. The esys_particle: a new 3-d discrete element model with single particle rotation. Advances in geocomputing. Springer; 2009. p. 183–228. [24] Cundall PA, Strack OD. A discrete numerical model for granular assemblies. Geotechnique 1979;29:47–65. [25] Hancock W, Weatherley D, Chitombo G. Modeling the gravity flow of rock using the discrete element method. In: MassMin 2012: 6th international conference and exhibition on mass mining, Canadian Institute of Mining, Metallurgy and Petroleum. [26] Wang Y, Mora P. Macroscopic elastic properties of regular lattices. J Mech Phys Solids 2008;56:3459–74. [27] Cacas MC, Ledoux E, de Marsily G, Tillie B, Barbreau A, Durand E, et al. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model. Water Resour Res 1990;26:479–89. [28] Cacas MC, Ledoux E, de Marsily G, Barbreau A, Calmels P, Gaillard B, et al. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 2. The transport model. Water Resour Res 1990;26:491–500. [29] Baecher G, Lanney N, Einstein H, et al. Statistical description of rock properties and sampling. In: The 18th US symposium on rock mechanics (USRMS), American Rock Mechanics Association.
References [1] Singh RN, Sun G. Applications of fracture mechanics to some mining engineering problems. Min Sci Technol 1990;10:53–60. [2] Atkins A. Modelling metal cutting using modern ductile fracture mechanics: quantitative explanations for some longstanding problems. Int J Mech Sci 2003;45:373–96. [3] Floreen S. The creep fracture of wrought nickel-base alloys by a fracture mechanics approach. Metall Trans A 1975;6:1741–9. [4] Bažant ZP, Cedolin L. Fracture mechanics of reinforced concrete. J Eng Mech Divis 1980;106:1287–306. [5] Han Z, Weatherley D, Puscasu R. Application of discrete element method to model crack propagation. In: ISRM Congress 2015 proceedings. [6] Han Z, Weatherley D, Puscasu R. A relationship between tensile strength and loading stress governing the onset of mode i crack propagation obtained via numerical investigations using a bonded particle model. Int J Numer Anal Methods Geomech 2017:1–13. [7] Rossmanith HP, Daehnke A, Nasmillner REK, Kouzniak N, Ohtsu M, Uenishi K. Fracture mechanics applications to drilling and blasting. Fatigue Fract Eng Mater Struct 1997;20:1617–36. [8] Gao H. Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials. Int J Fract 2006;138:101–37. [9] Kulatilake P, Ucpirti H, Wang S, Radberg G, Stephansson O. Use of the distinct element method to perform stress analysis in rock with non-persistent joints and to study the effect of joint geometry parameters on the strength and deformability of rock masses. Rock Mech Rock Eng 1992;25:253–74. [10] Bahaaddini M, Sharrock G, Hebblewhite B. Numerical investigation of the effect of joint geometrical parameters on the mechanical properties of a non-persistent jointed rock mass under uniaxial compression. Comput Geotech 2013;49:206–25.
225