Experimental and Numerical Studies on Interaction Mechanism between Joints and a Hole in Rock-like Materials under Uniaxial Compression

Abstract

To investigate the influence of joints on the stability of underground opening, uniaxial compression tests and FE analyses based on a microplane damage model for rocks has been conducted for rock-like models with a circular hole and a set of non-persistent joints. It was found that the peak strength and Young’s modulus decrease with the increase in joint continuity factor k, and variation of them with joint inclination angle β are W or V-shaped curves with the minima and maxima at β = 30° and 90°, respectively. The failure modes of the specimens and the collapse modes of the hole can be related to crack coalescence between the hole and the joints or matrix. Numerical simulation can reproduce the main features of macroscopic mechanical behavior and explain the anisotropic damage mechanism. The strong interaction between the hole and the nearest joint was revealed. During the elastic stage, stress concentration around the hole will be altered by the presence of the joints, and the effect may be strengthened with the increase in k. At the peak strength, the current stress concentration areas will be transferred from the hole surface to the interior due to stress loosening in damage localization bands/zones, and a higher hoop stress concentration factor may lead to lower strength.

1. Introduction

The existence of natural joints in rock mass will influence the damage process of surrounding rock mass due to excavation. Therefore, studying the interaction mechanism of cracks and holes is of great significance for the safe, efficient, and sustainable exploitation of underground energy resources (coal, oil, gas, geothermal, etc.). Excavation of an opening may cause disturbance of the original states of the surrounding rock mass, lead to stress redistribution, irreversible changes of the intrinsic rock structure, and formation of an Excavation Damage Zone (EDZ) [1]. Field experiments were conducted for the underground test tunnel in massive rocks or sparsely jointed rock masses, e.g., in brittle unfractured granite at Canada Underground Research Laboratory [2], in a Mesozoic shale formation of very low permeability at Switzerland [3], and in migmatite in Baishan hydropower station in China [4] have confirmed that the hydraulic conductivity in the EDZ was greatly increased, while the elastic wave velocity in the EDZ was decreased compared with that of the undisturbed rock.

It has been known that rock masses usually contain discontinuities at different scales, such as faults, bedding planes and joints. Many accidents occurring during construction and operation stages of tunnels can be attributed to the movement of faults or coalescence of non-persistent joints around the tunnel [5,6]. To systematically investigate the influence of discontinuities on the stability of underground openings, large-scale geomechanical model tests have been conducted by many researchers. For example, He et al. [7] studied the stability of a rectangular tunnel in horizontally stratified rocks at a depth of 1000 m in the Qishan coal mine by using the infrared thermography technique and found intensive cracks normal to the rock layers. Huang et al. [8] conducted overload tests for an underground tunnel with a weak interlayer and found that it may increase the size of the failure zones and cause asymmetry.

Compared to large-scale geo-mechanical model tests, small-scale physical model tests are more convenient for investigating the failure mechanism of a rock mass around an opening, especially for those with many discontinuous joints. For example, Lajtai and Lajtai [9] carried out triaxial compression tests for plaster blocks containing a circular hole. They found that the collapse of the hole can be divided into three stages, i.e., initiation of primary tensile and normal shear cracks from the hole, formation of inclined shear cracks within the crush zone, and coalescence of them with the secondary tensile cracks far from the hole. Wong et al. [10] conducted excavation for a circular opening in plaster models with non-persistent joints under biaxial compression. They found that tensile mode and shear or combined tensile–shear mode creep failure may occur for lower and higher confining pressure, respectively. Sagong et al. [11] carried out biaxial compression tests for a cement model with an opening and a set of almost continuous close joints to investigate the influence of joint dip angle (30°, 45°, and 60°) on the rock segment fracture and joint sliding behaviors. Yang et al. [12] conducted a uniaxial compression test for a rock-like material model with an opening and a set of non-persistent open joints (the joint continuity factor is fixed at 0.5), and the influence of joint orientation on the mechanical properties and crack coalescence from the hole, and the joints were investigated.

To learn more about the underlying failure mechanism for a rock mass around an opening, numerical studies have been conducted by many researchers, and some novel insights have been obtained. For example, Kawamoto et al. [13] and Swaboda et al. [14] conducted finite element (FE) modeling for underground tunnels in jointed rock masses based on their anisotropic damage model and found that damage analysis gives more satisfactory results with the measured data than those by the conventional analysis. Huang [8] carried out FE analysis for their physical model tests by incorporating an isotropic plastic damage model into ABAQUS and found that the size of the damage zone around the tunnel increased with the increase in the thickness of the weak interlayer. Yeung [15] conducted two-dimensional DDA modeling for a tunnel in a rock mass with two sets of continuous joints and found that rock block size may be more appropriate than joint spacing to measure tunnel stability in a blocky rock mass. Sagong et al. [11] and Yang et al. [12] carried out PFC2D modeling for their physical model tests and found that the damage zone around the hole decreased with the increase in the joint angle and tensile and shear cracks may occur mainly at the peripheral of the preexisting joints and the sidewalls of the circular hole, respectively.

The above-mentioned experimental and numerical studies have improved our understanding of the failure mechanism of jointed rock masses with an opening. However, the influence of joint persistence has received less attention in the literature, even though its effect on rock masses’ strength has long been recognized and investigated [16,17,18,19,20]. Furthermore, consideration of the anisotropic damage effect of joints or cracks under tension and compression has turned out to be rather difficult and not straightforward in the classical macroscopic tensorial approach [21].

To efficiently characterize the anisotropic inelastic response of quasi-brittle materials (whose main damage mechanism is crack evolution, opening/closing, and slip), Bažant [22] proposed a vectorial approach in constitutive modeling called the “microplane model”. The main feature of this vectorial approach is that the constitutive law is formulated as a relation between the stress and strain vectors on a microplane (a plane of arbitrary orientation of the material), while the tensorial stress and strain relation can be obtained by superimposing the contributions from all microplanes in a suitable manner. Different microplane models have been developed for concrete [23,24,25], clays [26], soils [27], and rocks [28]. To consider the resistance of cracks or joints under compression, Chen and Bažant [29] proposed a microplane damage model (MPD model) for rocks by introducing the two-phase concept (the rock matrix and the rock joint phases) on microplanes. Numerical simulations for triaxial compression tests of sandstone [30] and uniaxial tests on jointed plaster mortar [13] have demonstrated that the proposed model not only can characterize the main features of rocks (anisotropy, strain softening, dilatancy, and brittle-ductile transition) but also can explain the underlying anisotropic damage mechanisms. The MPD model for rocks has been implemented into a commercial FEM software ABAQUS through a VUMAT subroutine [31].

In this study, the combined influence of joint persistence and orientation on the response of a rock mass containing an opening will be investigated experimentally and numerically. Uniaxial compression tests will be conducted for rock-like models with a circular hole and a set of non-persistent joints, and the dependence of strength, deformability, and cracking process on the two geometrical parameters will be studied first. Then, the underlying failure mechanism, including stress concentration, damage evolution, stress redistribution, and the interaction between the hole and the joints, will be investigated through FE modeling based on Chen and Bažant’s MPD model for rocks.

2. Experimental Methodology

2.1. Setup of Jointed Specimens with a Hole

The dimensions of the jointed specimens containing a 50 mm diameter cylindrical hole are 300 mm × 300 mm × 50 mm in width, length, and thickness, respectively. Since the length of the specimen is six times the diameter of the hole, the stress concentration caused by the hole can be ignored for the surrounding rock beyond the edge. A set of non-persistent open joints penetrated through the thickness direction are regularly arranged in an area at a distance of 20 mm from the specimen boundary, as shown in Figure 1. The joint spacing s and the joint center distance c were both fixed at 50 mm. The two geometrical parameters of the joint set, i.e., the joint continuity factor k (defined as the ratio of the joint length L_j to the joint center distance c, i.e., k = L_j/c) and the joint inclination angle β (defined as the angle between the joint plane and the horizontal plane), were investigated. Three values of joint continuity factor k, namely, 0.2, 0.5, and 0.8, were investigated and named as groups HB, HC, and HD, respectively. Five values of joint inclination angle β, i.e., 0°, 30°, 45°, 60°, and 90°, were studied for each value of k. For comparison, the intact specimens without or with the hole were also investigated and named A and HA, respectively. Therefore, a total of 17 series of specimens were tested and listed in

Table 1. Joint geometrical parameters of the tested specimens.

Number	Serie	With an Opening	k	β	Lj (mm)
1	A	No	-	-	-
2	HA	Yes	0.0	-	0
3	HB0	Yes	0.2	0°	10
4	HB30	Yes	0.2	30°	10
5	HB45	Yes	0.2	45°	10
6	HB60	Yes	0.2	60°	10
7	HB90	Yes	0.2	90°	10
8	HC0	Yes	0.5	0°	25
9	HC30	Yes	0.5	30°	25
10	HC45	Yes	0.5	45°	25
11	HC60	Yes	0.5	60°	25
12	HC90	Yes	0.5	90°	25
13	HD0	Yes	0.8	0°	40
14	HD30	Yes	0.8	30°	40
15	HD45	Yes	0.8	45°	40
16	HD60	Yes	0.8	60°	40
17	HD90	Yes	0.8	90°	40

.

2.2. Specimen Preparation and Test Procedure

The model material used in this study is a mixture of gypsum, Portland cement, and water at a ratio of gypsum: cement: water = 0.99:0.01:0.6 in weight. In addition, a retarder with a dose of 0.05% in weight was added to slow down the setting of the model material. The physical and mechanical properties of the model material are listed in

Table 2. Physical and mechanical properties of the model material.

Density ρ (g/cm−3)	Young’s Modulus E (GPa)	Poisson’s Ratio v	UCS σc (MPa)	Brazilian Tensile Strength σt (MPa)	Cohesion c (MPa)	Friction Angle φ (°)
1.22	3.17	0.23	10.86	2.25	2.2	34

, which were obtained through the uniaxial compression test, standard triaxial compression test, and Brazilian test.

Figure 2 shows the mold and fabrication process of a specimen. The mode consists of a cover plate with precut slots and a base plate made of PMMA (polymethyl methacrylate), and a frame and a cylinder made of stainless steel. To make the pre-existing joints, a group of 0.3 mm-thick nickel steel sheets was inserted into the mixture through the precut slots and removed after the start of the liquid mixture solidification. The specimens were cured for 21 days at room temperature. In order to observe repeatable results, the mixing, fabrication, and curing of the material were carefully controlled, and at least two samples with the same joint configuration were prepared for the tests. The mechanical parameters of each series were taken as the average of the two specimens.

The uniaxial compression tests were performed in an INSTRON 8506 servo-controlled hydraulic loading system with high stiffness. To obtain complete stress–strain curves, displacement control was applied with a constant loading rate of 0.0025 mm/s. The load and displacement were recorded automatically by the loading system, and fracture propagation on the surface of the specimens were monitored by a digital video and recorded by a high-resolution digital camera.

3. Laboratory Test Results

3.1. Strength and Deformability

Figure 3 shows the laboratory test results for the axial stress–strain curves of the two intact specimens with and without a hole and the jointed specimens with a hole. Chen et al. [21] have classified the stress–strain curves of jointed specimens under uniaxial compression into four types according to their different patterns after the first peak, i.e., Type 1—single-peak curve, Type 2—general strain softening with oscillations, Type 3—yield platform-strain softening, and Type 4—yield platform-strain hardening-strain softening. In this study, except for Type 4 deformation behavior, single-peak curve (Type 1) and multi-peak deformation behaviors of Type 2 and 3 were also observed for the jointed specimens with the hole under uniaxial compression. Type 1 deformation behavior was observed for specimens A, HA, HB45, HD45, HB60, HB90, and HC90, Type 2 deformation behavior for specimens HD0, HC30, and HD30, and Type 3 deformation behavior for specimens HB0, HC0, HB30, HC45, HC60, and HD90.

The peak strength σ_p of specimens A and HA are 11.94 MPa and 10.14 MPa, respectively, and Young’s moduli E (the tangent modulus in the linear elastic stage of an axial stress–strain curve) of the two specimens are 3.84 GPa and 3.02 GPa, respectively. Figure 4 and Figure 5 show the variation of σ_p and E with the two geometrical parameters, i.e., the joint inclination angle β and the joint continuity factor k. It can be seen that: (1) the peak strength and Young’s modulus of the specimen may decrease significantly in the presence of the hole and the non-persistent open joints; (2) for each value of k, the curves of σ_p vs. β and E vs. β are W-shaped or V-shaped with the minima and maxima occurring at β = 30° and 90°, respectively; (3) for a constant value of β, the peak strength σ_p and Young’s modulus E will decrease with increase in the joint continuity factor k, and the fastest and slowest velocity occurs at β = 30°and β = 90°, respectively.

3.2. Cracking Process

The anisotropic nonlinear behavior of the jointed specimens with the hole is related to their different cracking process. Figure 6 shows the observed crack initiation types, which are classified according to their mechanism (tensile or shear, represented by T and S, respectively) and position (represented by their lower index), namely, from the matrix (T_m), joint (T_j1, T_j2, and S_j), and periphery of the hole (T_h and S_h). Figure 7 shows the observed crack coalescence types between the joints and between the hole and the joint or matrix (initial by H), respectively.

Crack initiation type T_m is produced by an axial tensile crack that emanates from the matrix; crack initiation type T_j1 is produced by a tensile crack or a wing crack that emanates from the tip or middle of the joint, propagates perpendicularly first and then toward the loading direction; crack initiation type T_j2 is a quasi-coplanar tensile crack that emanates from the joint tip, which only occurs in specimens with vertical joints; crack initiation type S_j is a quasi-coplanar shear crack that emanates from the joint tip; crack initiation T_h is a tensile crack that emanates from the roof or floor of the hole and propagates toward the loading direction; crack initiation type S_h is a shear crack or crush zone having finite width that emanates from the left or right sides of the hole and propagates perpendicular to the loading direction.

Two types of tensile crack coalescence between the joints were observed, namely, TT₁ and TT₂, which were formed by the linkage of two type T_j1 tensile cracks and two type T_j2 tensile cracks, respectively. Coalescence Type SS is formed by the linkage of two type S_j shear cracks along the same joint plane. Two types of mixed tensile–shear crack coalescence between the joints were observed, namely, TS₁ and TS₂, which were formed by _the linkage of one Type T_j1 tensile crack and another Type S_j shear crack that started from a different joint plane and the same joint plane, respectively.

Three types of tensile crack coalescence between the hole and the joint or matrix were observed, namely, HTT_m, HTT_j1, and HTT_j2, which were formed by the linkage of a Type T_h tensile crack from the hole and a tensile crack of Type T_m from the matrix, or T_j1 and T_j2 from the joint, respectively. Two types of mixed tensile–shear crack coalescence between the hole and the joint or matrix were observed, namely, HST_m and HST_j, which is formed by the linkage of a Type S_h shear crack from the hole and a tensile crack of type T_m or T_j1, respectively. Coalescence type HSJB is produced by the occurrence of the two coalescence types, HST_j and HTT_j1, from the same joint, which may cause the formation of a block.

Table 3. Crack coalescence types and collapse modes of the hole observed in the jointed specimens with a hole (Symbols ★ and ▲ represent crack coalescence types between the joints without or with participation of shear crack, respectively; ● and □ represent crack coalescence types between the hole and the joint or matrix without or with participation of shear crack, respectively).

Types	TT1	TT2	SS	TS1	TS2	HTTm	HTTj1	HTTj2	HSTm	HSTj	HSTB	Collapse Mode
HB0	★								□	□		Ih
HC0	★					●				□		Ih
HD0	★						●			□		IIh
HB30	★					●				□		Ih
HC30	★			▲		●				□		IIh
HD30	★			▲			●			□		IIh
HB45	★					●			□			Ih
HC45	★			▲			●			□	□	IIIh
HD45			▲		▲		●			□	□	IIIh
HB60	★						●			□		IIh
HC60	★		▲	▲			●			□	□	IIIh
HD60	★		▲	▲	▲		●			□		Ih
HB90		★							□	□		IIh
HC90		★						●		□		Ih
HD90		★						●		□		Ih

summarizes crack coalescence types observed in the jointed specimens with a hole. Tensile crack coalescence between joints of type TT₁ and TT₂ was found in specimens with horizontal joints or incline joints and those with vertical joints, respectively. Shear cracks coalescence between joints of type SS was found in specimens HD45, HC60, and HD60. Mixed shear–tensile cracks coalescence of type TS₁ and TS₂ were found in the specimens with inclined joints at medium and larger continuity factors (k = 0.5 and 0.8), i.e., HC30, HD30, HC45, HD45, HC60 and HD60. Tensile cracks coalescence around the hole, i.e., Types HTT_m, HTT_j1 and HTT_j2 were found in all jointed specimens except for specimens HB0 and HB90. Mixed shear–tensile cracks coalescence around the hole of Types HST_m and HST_j were found in all jointed specimens, while coalescence type HSB with formation of blocks was found only in specimens HC45, HD45, and HC60.

3.3. Collapse Modes of the Hole

Crack coalescence around the hole may lead to a collapse of the hole. Lajtai and Lajtai [9] have found that the linkage of the crush zone and secondary tensile cracks in the rock matrix (Coalescence type HST_m) may finally lead to a collapse of the hole in the intact rock samples. The observed collapse modes of the hole can be classified into three types as follows (see Figure 8):

(1)

Mode I_h—without removable blocks, where the linked tensile cracks or shear–tensile cracks around the hole may extend to the boundary of the specimen, and therefore, no removable blocks can be formed;

(2)

Mode II_h—with removable blocks formed by upper joints, where the linked tensile cracks or shear–tensile cracks around the hole may extend to the upper joints, leading to the formation of one or several removable blocks;

(3)

Mode III_h—with removable blocks formed by the nearest inclined joints, where coalescence type HSTB may occur from the nearest inclined joints to the hole, leading to the formation of a pair of removable blocks.

The collapse modes of all jointed specimens with a hole are listed in Table 3. Collapse mode III_h occurs in specimens with coalescence type HSTB, i.e., HC45, HC60, and HD45. Collapse mode II occurs either in specimens with low joint inclination angle (β = 0° and 30°) at the medium or the largest joint continuity factors (k = 0.5 or 0.8), i.e., HD0, HC30, and HD30, or in specimens with the high-joint inclination angles (β = 60° and 90°) at the smallest joint continuity factor (k = 0.2), i.e., HB60 and HB90. Collapse mode I occurs in the other specimens.

Figure 9, Figure 10 and Figure 11 show the cracking process of the specimens with different collapse modes, i.e., specimens HC90, HD0, and HC60. For specimen HC90 of collapse mode I_h (see Figure 9), at the pre-peak stage, Type T_j2 tensile cracks initiated from joints at the right side of the hole first (see Figure 9a), and then Type T_h tensile cracks and Type S_h shear cracks initiated from the hole, as well as Type T_j2 wing crack initiated from joints near the left side(see Figure 9b). At the post-peak stage, the linkage of shear cracks at the two sides of the hole with the wing cracks from the joints (coalescence Type HST_j) finally leads to the collapse of the hole (see Figure 9c).

For specimen HD0 of collapse mode II_h (see Figure 10), at the pre-peak stage, Type T_j1 tensile cracks initiated from joints at the two sides of the hole first (see Figure 10a), and then Type S_h shear cracks initiated from the two sides of the hole (see Figure 10b). At the post-peak stage, Type T_h tensile cracks initiated from the roof and floor of the hole and type T_j1 tensile cracks initiated from the upper joints (see Figure 10c), and the linkage of those cracks (coalescence Types HST_j and HTT_j1) led to the formation of the upper removable blocks and, finally, to the collapse of the hole (see Figure 10d).

For specimen HC60 of collapse mode III_h (see Figure 11), at the pre-peak stage, Type T_j1 wing cracks initiated from the nearest inclined joints, and type T_m tensile cracks initiated in the matrix first (see Figure 11a). At the post-peak stage, Type T_h tensile cracks initiated from the upper right and lower left of the hole, and several type T_j1 tensile cracks initiated from other inclined joints around the hole (see Figure 11b). Type S_h shear cracks initiated from the two sides of the hole, and the linkage of these cracks led to the formation of a pair of removable blocks and, finally, to the collapse of the hole (see Figure 11c).

3.4. Failure Modes

Propagation and coalescence of cracks together with the pre-existing joints may cause the formation of continuous fractures, which may finally lead to the failure of the specimen. The observed failure modes in the jointed specimen with a hole can be classified into four types as follows (see Figure 12):

(1)

Mode I—axial cleavage, in which linked vertical tensile fractures may break the specimens into several parts or pillars;

(2)

Mode II—crushing, in which linked vertical tensile fractures together with preexisting low inclination angle joints may break the specimens into many small blocks;

(3)

Mode III—stepped, in which mixed tensile–shear cracks and the preexisting inclined joints in neighbor columns form one or several stepped failure planes;

(4)

Mode IV-sliding, in which quasi-coplanar shear cracks and the preexisting joints form one or several failure planes and lead to sliding along those joint planes.

It should be noted that the four types of failure modes also have been found for jointed specimens without a hole in many pieces of the literature, for example, Prudencio and Jan [32] and Chen et.al. [33].

Figure 13 presents the failure modes of all jointed specimens with a hole and their fractures. Typical specimens of different failure modes are marked with the red border, i.e., specimens HC90, HD0, HC30, and HD45 with failure modes I to IV, respectively.

Failure mode I occurs in specimens at the smallest joint continuity factor (k = 0.2) or with vertical joints (β = 90°), namely, specimens HB0, HB30, HB45, HB60, HB90, HC90, and HD90. Failure mode II occurs in specimens with horizontal joints (β = 0°) at medium and the largest joint factors (k = 0.5 and 0.8), namely, specimens HC0 and HD0. Failure mode III occurs in specimens with inclined joints (β = 30°, 45°, and 60°) at medium joint factor (k = 0.5) or specimens with β = 30° at the largest joint factor (k = 0.8), namely, specimens HC30, HC45, HC60, and HD30. Failure mode IV occurs in specimens with β = 45° and 60° at the largest joint factor (k = 0.8), namely, specimens HD45 and HD60. In general, specimens with failure modes III or IV that are characterized by involving shear crack propagation may have lower peak strength than those with failure modes I and II.

4. FE Analysis Based on the MPD Model

4.1. A brief Introduction to the MPD Model for Rocks

Figure 14 shows the flowchart of the MPD model for rocks in a VUMAT user subroutine of ABAQUS. For a representative volume element (RVE) of a rock or a rock mass, a kinematic constraint is assumed for correlations between macroscopic and microplane stress/strain. Namely, the microplane strain vector of the rock mass ${\mathbf{\epsilon }}^{mic}$ is the directional projection of the macroscopic strain tensor $\mathbf{\epsilon }$, while the macroscopic stress tensor $\mathbf{\sigma }$ can be obtained by the directional integral of all microplane stress vectors of the rock mass ${\mathbf{\sigma }}^{mic}$, according to the principle of variational work, which can be given by the following:

$${\epsilon }_N={\epsilon }_{ij}{N}_{ij},{\epsilon }_M={\epsilon }_{ij}{M}_{ij},{\epsilon }_L={\epsilon }_{ij}{L}_{ij};$$

(1a,b,c)

$${\sigma }_{ij}=\frac{3}{4\pi }{\displaystyle \int \left({\sigma }_N{N}_{ij}+{\sigma }_M{M}_{ij}+{\sigma }_L{L}_{ij}\right)\mathrm{d}S}$$

(2)

$$N_{ij}=n_i{n}_j,M_{ij}=\left(m_i{n}_j+m_j{n}_i\right)/2,L_{ij}=\left(l_i{n}_j+l_j{n}_i\right)/2.$$

(3)

where ${\sigma }_{ij}$ and ${\epsilon }_{ij}$ are tensorial components of the macroscopic stress and strain tensor, respectively, where the subscripts i, j = 1, 2, 3 refer to the global Cartesian coordinate $x_i$; $n_i$, $l_i$, and $m_i$ are the unit normal vector of the microplane and the two orthogonal unit vectors on the microplane; ${\sigma }_N,{\sigma }_M,{\sigma }_L$, and ${\epsilon }_N,{\epsilon }_M,{\epsilon }_L$ are the normal and the two shear components of microplane stress and strain vectors; ${\displaystyle \int \mathrm{d}S}={\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }\sin \varphi \mathrm{d}\varphi }\mathrm{d}\theta }$ is the integral over the surface of a unit sphere, where $\varphi$ and $\theta$ are the spherical angles. Here, Einstein summation convention is adopted, i.e., repetition of subscripts implies summation over 1, 2, and 3.

A microplane of the rock mass can be seen as a combination of the two phases, namely, the rock matrix and the rock joint. The ratio of the joint phase on the microplane ${\omega }_N=\omega \left(\mathit{n}\right)$ is selected as the damage variable. On each microplane, the two phases are mechanically coupled in parallel; namely, the microplane strain vector of each phase is identical to that of the rock mass ${\mathbf{\epsilon }}^{mic}$, while the microplane stress vector of the rock mass ${\mathbf{\sigma }}^{mic}$ is the summation of that of the two phases ${\mathbf{\sigma }}^{\left(K\right)}(K=R,J)$ weighted by their ratios:

$${\mathbf{\epsilon }}^{mic}={\mathbf{\epsilon }}^{\left(R\right)}={\mathbf{\epsilon }}^{\left(J\right)}$$

(4)

$${\mathbf{\sigma }}^{mic}=\left(1-{\omega }_N\right){\mathbf{\sigma }}^{\left(R\right)}+{\omega }_N{\mathbf{\sigma }}^{\left(J\right)}$$

(5)

On the microplane, the elastic behavior of the two phases can be described by their linear elastic laws. The inelastic behavior of the two phases on the microplane under tension, compression, and shear loading can be described by their corresponding strain-dependent stress boundaries. The details of the linear elastic laws and strain-dependent stress boundaries of the two phases can be found in Chen and Bažant [29]. For each phase, the elastic trail stress will be returned to the stress boundaries if it exceeds the boundary.

To reflect the initiation and propagation of tensile or shear cracks due to the accumulation of volumetric expansion, normal extension, or shear deformation on microplanes, the microplane damage evolution law was given by the power law of the historic strains as follows:

$${\omega }_N=1-e^{-\left[{\left({\epsilon }_V^{-h}/a_1\right)}^{q_{{}_1}}+{\left({\epsilon }_D^{-h}/a_2\right)}^{q_{{}_2}}+{\left({\epsilon }_T^{-h}/a_3\right)}^{q_{{}_3}}\right]}$$

(6)

where ${\epsilon }_V^{-h}, {\epsilon }_D^{-h}, {\epsilon }_T^{-h}$ are the maximum volumetric expansion strain, deviatoric tensile strain, and shear strain reached so far on the microplane, respectively; a₁, a₂, a₃, q₁, q_2, and q₃ material constants control damage evolution related to ${\epsilon }_V^{-h}, {\epsilon }_D^{-h}, {\epsilon }_T^{-h}$.

To measure the overall damage of the RVE, the directional average of ${\omega }_N$ for all microplanes, namely, damage ω₀, was also calculated and can be given by the following:

$${\omega }_0=\frac{1}{4\pi }{\displaystyle \int {\omega }_N}\mathrm{d}S\approx 2{\displaystyle \sum _{m=1}^{N_m}{w}_m{\omega }_N{}^{\left(m\right)}}$$

(7)

where $N_m$ is the number of the numerical integration points on a unit hemisphere surface, $w_m$ is the weight of the m-th microplane, and ${\omega }_N{}^{\left(m\right)}$ is the damage variable on the m-th microplane.

The value of the microplane damage ${\omega }_N$ and the average damage ω₀ are always between 0 and 1. For a microplane, ${\omega }_N=0$ and ${\omega }_N=1$ represent no crack/joint and a continuous crack/joint on the microplane, respectively. For the RVE, ${\omega }_0=0$ and ${\omega }_0=1$ represent no damage (the intact RVE without any crack/joint) and entire damage (the RVE full of joint material), respectively.

4.2. FE Model and Calibration of the Material Parameters

The size of the numerical model is the same as that in the physical experiment. Figure 15 shows the FE model of specimen HC45 (κ = 0.5 and β = 45°). The top and bottom loading plates are rigid bodies, loaded at a constant displacement rate and fixed, respectively. The pre-existing joints are modeled by a rectangle opening with a thickness of 0.3 mm. The element type used here is a six-node linear triangular prism (C3D6), and the total number of elements is about 200,000.

The twenty material parameters of the MPD model used in numerical modeling are calibrated by specimens A and HA (the intact specimens with or without a hole) and are listed in

Table 4. Material parameters of the MPD model used in numerical modeling.

Parameter	Value	Meaning or Function
E(R)	3.2 GPa	Young’s modulus of the rock matrix
E(J)	3.2 GPa	Young’s modulus of the rock joint or crack
v(R)	0.23	Poisson’s ratio of the rock matrix
v(J)	0.23	Poisson’s ratio of the rock joint or crack
T(R)	3.2 MPa	Microplane tensile strength of the rock matrix
T(J)	0 MPa	Microplane tensile strength of the rock joint
α0	5.0	Ratio of the minimum compressive strength to tensile strength of the rock matrix
βc	0.5	Maximum value of the joint mobilization factor of the rock joint
ε0V	0.001	Volumetric strain threshold for strain hardening of the rock matrix
ε0N	0.0005	Normal strain threshold for strain hardening of the rock matrix
c1	0.2	Control dependence of compressive stress boundary of the rock matrix phase on εN
c2	0.005	Control dependence of compressive stress boundary of the rock matrix phase on εV
c3	0.001	Control dependence of shear stress boundary on normal stress of the two phases
c4	0.05	Control dependence of shear stress boundary on the confinement of the two phases
a1	0.007	Control dependence of damage evolution on volumetric expansion
a2	0.006	Control dependence of damage evolution on deviatoric tensile strain accumulation
a3	0.05	Control dependence of damage evolution on shear strain accumulation
q1	1.5	Control speed of damage evolution due to volumetric expansion
q2	1.5	Control speed of damage evolution due to deviatoric tensile strain accumulation
q3	1.0	Control speed of damage evolution due to shear strain accumulation

. The calibration method for these parameters can be given as follows: (1) Young’s modulus and Poisson’s ratio of the rock matrix, $E^{\left(R\right)}$ and ${\nu }^{\left(R\right)}$, can be obtained from Table 2. The same values were taken for those of the rock joint or crack $E^{\left(J\right)}$ and ${\nu }^{\left(J\right)}$ for simplicity; (2) parameters ε⁰_V, ε⁰_N, c₁, c₂, c₃, c₄, a₁, a₂, a₃ q₁, q_2, and q₃ mainly control the shape of the nonlinear stress–strain curve and can be set as the similar values for a certain type of rocks; (3) parameters T^(R), T^(J), α_0, and β_c can be obtained from optimal fitting to the stress–strain curves of the specimens.

4.3. Numerical Results of Macroscopic Mechanical Response

Figure 16 shows the axial stress–strain curves of specimens A and HA (the intact specimens with or without a hole, respectively) as well as all the jointed specimens with a hole obtained by numerical simulation. Figure 17 and Figure 18 show curves of the peak strength and Young’s modulus vs. the two geometrical parameters by numerical simulation.

In Figure 16a, Figure 17a and Figure 18a, results by laboratory tests were also depicted for comparison. In general, the numerical results can reproduce the main features observed in laboratory tests, including dependence of the deformation behavior, peak strength, and Young’s modulus on the joint inclination angle β and the joint continuity factor k.

4.4. Failure Mechanism of the Intact Specimen with the Hole

The nonlinear mechanical response of specimen HA under uniaxial compression can be related to stress concentration and damage evolution in the matrix around the hole. The problem of the disturbance of the elastic stress field caused by the presence of the hole under far-field uniaxial compression was solved by Kirsch in 1898. Accordingly, the radial and hoop stress concentration factor (σ_r/σ₁) and (σ_θ/σ₁) can be given as follows (see the textbook of Jaeger et al. [34]):

$$\begin{array}{l}\frac{{\sigma }_{\mathrm{r}}}{{\sigma }_1}=\frac{1}{2}\left[\left(1-\frac{a^2}{r^2}\right)+\left(1-4\frac{a^2}{r^2}+\frac{3a^4}{r^4}\right)\cos 2\theta \right]\\ \frac{{\sigma }_{\mathsf{\theta }}}{{\sigma }_1}=\frac{1}{2}\left[\left(1+\frac{a^2}{r^2}\right)-\left(1+3\frac{a^4}{r^4}\right)\cos 2\theta \right]\end{array}$$

(8a,b)

where r is the distance to the hole center; a is the radius of the circular hole; θ is the angle of rotation from the axial stress σ₁; σ_r and σ_θ are the radial and hoop stresses in the specimen.

The maximum tensile and compressive stress may occur along the hole periphery, with ${\sigma }_{\mathsf{\theta }}/{\sigma }_1=-1$ at θ = 0° or 180° (the roof or floor) and ${\sigma }_{\mathsf{\theta }}/{\sigma }_1=3$ at θ = 90° or 270° (the right or left side), respectively.

Figure 19a,b shows hoop stress σ_θ in specimen HA at the early elastic deformation stage, when σ₁ = 1 MPa and the peak strength, when σ₁ = σ_p, respectively. Figure 19c,d shows damage ω₀ at σ₁ = σ_p and fractures observed in the experiment. Figure 20 and Figure 21 show the variation of radial and hoop stress concentration factors (σ_r/σ₁) and (σ_θ/σ₁) and the average damage of all microplanes ω₀ vs. distance ratio r/a along a horizontal or vertical plane in specimen HA at σ₁ = 1 MPa and σ₁ = σ_p, respectively.

At the early elastic deformation stage, when σ₁ = 1 MPa, tensile or compressive stress concentrate at the roof/floor and the left/right sides of the hole periphery (see Figure 19a), respectively. The results of (σ_r/σ₁) and (σ_θ/σ₁) vs. r/a by numerical simulation are very close to those by theoretical solution for the elastic solid since damage ω₀ remains almost zero when σ₁ = 1 MPa (see Figure 20a and Figure 21a).

At the peak strength, when σ₁ = σ_p, damage evolution due to tensile or compressive stress concentration is localized in two narrow vertical bands above or below the hole or in two small horizontal zones at the two sides of the hole (see Figure 19c), respectively. As shown in Figure 20b and Figure 21b, the length of the vertical damage bands and the horizontal damage zones are about 1.5a and 0.5a, respectively. Due to stress loosening in damage localization zones, the current tensile or compressive stress concentration areas will move to the front of the damage zones (see Figure 19b). In general, fractures in specimen HA coincide very well with the damage localization areas.

4.5. Failure Mechanism of the Jointed Specimens with the Hole

The nonlinear mechanical response of the jointed specimens with the hole under uniaxial compression may be much different from that of the intact specimen with the hole since the presence of the joints may influence stress concentration, damage evolution, and inelastic deformation development.

4.5.1. Influence of Joints on Stress Concentration at Elastic Stage

For a single elliptic crack in an infinite elastic solid subject to a far-field uniaxial compression stress σ₁, the stress concentration around the crack has been solved by Muskhelishvili in 1963 and can be found in the textbook of Jaeger et al. [34]. Lajtai [35] discussed the influence of the axial ratio b/a (the ratio of the smaller to the larger semiaxis) and the inclination angle β on the stress concentrations along an elliptic crack periphery, and the extent and shape of tensile zone beyond the periphery of a flat crack (b/a = 0).

In this study, the length of the joints in specimen groups HB, HC, and HD is 10, 25, and 40 mm, respectively, while the joint aperture is fixed at 0.3 mm for all specimen groups, and the axial ratio b/a of the equivalent elliptic crack is 0.03, 0.012, and 0.075, respectively. The minimum and maximum principal stress concentration factors σ_min/σ₁ and σ_max/σ₁ around the elliptic crack have been calculated for b/a = 0.03 and β = 0°, 45°, and 90° and are plotted in Figure 22.

For β = 0° and 45°, the elliptic crack produces large tension zones around the whole crack periphery with σ_min/σ₁ = −1 and near the two ends with σ_min/σ₁= −4.7, respectively; it induces small zones of very high compressive stress concentration around the crack tips, with σ_max/σ₁ = 67.7 and 22.3, respectively. For β = 90°, the elliptic crack produces very small tensile zones at the crack tips with σ_min/σ₁ = −1 and nearly no compressive stress concentration around the crack periphery with σ_max/σ₁ = 1~1.06. With the decrease in the axial ratio b/a, the magnitude of the maximum tension for the inclined elliptic crack and the order of magnitude of the maximum compression for the inclined and horizontal elliptic crack will increase. The stress concentration around the elliptic crack can explain the crack initiation around the joints, including tensile crack initiation from the middle of the horizontal joints and wing crack initiation from the inclined joints, and shear crack initiation from the tips of the inclined joints (see Figure 6).

For the specimens with the hole and multiple joints, stress concentration around the hole and the joints will be affected by their interaction. At the early elastic stage when σ₁ = 1 MPa, Figure 23 and Figure 24 show the distribution of hoop stress σ_θ and hoop stress concentration factor (σ_θ/σ₁) along the hole periphery in the specimens with k = 0.2 and 0.8 at β = 0°, 45°, and 90°, respectively.

When the joint continuity factor is the smallest, i.e., k = 0.2, either the interaction among the small joints or the interaction between the joints and the hole is not obvious. The joints in these specimens only produce stress concentration in the small areas around themselves; therefore, the joint orientation may have little influence on the stress distribution in the other area of the specimen. For all values of β, the values of (σ_θ/σ₁)_min are between (−0.94~ −1.18) and occur at θ = 0° (the roof) or 180° (the floor), while the values of (σ_θ/σ₁)_max are between (3.23~3.51) and occur at θ = 90° (the right side) or 270° (the left side). When the joint continuity factor is the largest, i.e., k = 0.8, the interaction among the large joints and the interaction between the joints and the hole may become significant, except for the vertical joints. For specimen HD90 (k = 0.8 and β = 90°), tensile stress concentration increases slightly in the presence of the nearest vertical joints near the top and the bottom, while compressive stress concentration changes a little, with (σ_θ/σ₁)_min = −1.18 and (σ_θ/σ₁)_max = 3.64 occurring at the roof/floor and the two sides of the hole, respectively. For specimen HD0 (k = 0.8 and β = 0°), a small decrease in tensile stress concentration and a significant increase in compressive stress concentration due to adjacent horizontal joints were found, with (σ_θ/σ₁) _min = −0.78 at the roof/floor and (σ_θ/σ₁)_max= 5.29 at the two sides, respectively. For specimen HD45 (k = 0.8 and β = 45°), nearly no changes in tensile stress concentration and a salient increase in compressive stress concentration were found, with (σ_θ/σ₁)_min = −1.04 at θ = 346.4° or 166.4° nearby the roof/floor and (σ_θ/σ₁)_max= 4.42 at θ = 58.7°or 234.2° close to the nearest inclined joints, respectively. In general, the nearest joints, which are located on the intersecting joint plane to the hole, may have a significant influence on the stress concentration around the hole.

4.5.2. Influence of Joints on Damage Evolution and Stress Loosening

At peak strength, when σ₁ = σ_P, Figure 25 and Figure 26 show the distribution of hoop stress σ_θ and damage ω₀ in the specimens, with k = 0.2 and 0.8 at β = 0°, 45°, and 90°, respectively.

Similar to the intact specimen containing the hole (specimen HA), damage developed from the initial tensile or compressive stress concentration areas and propagated into localization zones, which could lead to stress loosening in these zones and transfer of the current stress concentration areas from the hole neighborhood to the interior.

Even for the specimens with the smallest joint continuity factor (k = 0.2), joint orientation has a salient influence on damage localization and stress concentration. For specimens with horizontal or vertical joints, i.e., HB0, HD0, HB90, and HD90, narrow vertical tensile damage localization bands and wide compressive horizontal damage localization zones developed from the roof/floor and the two sides of the hole, respectively. For specimens with inclined joints, i.e., HB45 and HD45, the current tensile stress concentration areas still located near the roof and floor since no damage developed in these areas; oblique compressive damage localization zones developed from the two sides (for specimen HB45) or the upper left/lower right sides (for specimen HD45) of the hole linked with coplanar damage zones in the rock bridges along the intersecting joint plane to form a failure plane. With the increase in joint continuity factor, widths of compressive damage localization zones and the length of tensile damage bands around the hole may decrease due to the decrease in the rock bridge length.

Figure 27 shows the evolution of hoop stress and damage of elements 1 (σ_θ1 and ω₀₁) and 2 (σ_θ2 and ω₀₂) in the specimens with k = 0.2 and 0.8 at β = 0°, 45°, and 90°, respectively. For each specimen, tensile and compressive strain softening happens in elements 1 (at the top) and 2 (at the right side), respectively. Due to the high compressive stress concentration, damage in element 2 grows earlier and faster than damage in element 1 (except for specimen HD90, due to the influence of the adjacent vertical joints); the peak stress in the two elements occurs much earlier than the peak strength of the specimen (except for specimen HD45).

For specimens with horizontal or vertical joints, i.e., HB0, HD0, HB90, and HD90, due to stress concentration in the two elements at the early elastic stage, damages ω₀₁ and ω₀₂ grow rapidly and reach the maximum value 1 (represent entire damage on all microplanes) before the peak strength leads to strain softening and stress loosening of the two elements (see Figure 27a,b,e,f). For specimens with inclined joints, i.e., HB45 and HD45, damage ω₀₂ grows slowly during the whole loading process (damage starts from element 2 firstly and then propagates toward the nearest inclined joint, see Figure 26), leading to partial stress loosening of the two elements at peak strength. For specimen HD45, due to the deflection of stress concentration areas at the early elastic stage, the average damage ω₀ of the two elements grows slowly and is not large (in specimen HB45) or very small (in specimen HD45) during the whole loading process, leading to partial or no stress loosening of these elements at peak strength. Figure 28 shows the variation of the hoop stress concentration factor of elements 1 (σ_θ1/σ₁) and 2 (σ_θ2/σ₁) with joint inclination angle β at peak strength σ₁ = σ_P for all jointed specimens with the hole. For each joint continuity factor, curves of σ_θ1 /σ₁ vs. β and σ_θ2/σ₁ vs. β are anti-V or U, or M shaped with the maxima at β = 30° or 45°, or 60° and the minima at β = 0° or 90°. The hoop stress concentration factor of the two elements σ_θ1 /σ₁ and σ_θ2/σ₁ increases with the increase in the joint continuity factor k. It seems that the hoop stress concentration factor at peak strength can be related to the strength of the specimen; namely, a higher hoop stress concentration factor at σ₁ = σ_P corresponds to a lower strength.

Figure 29 shows damage ω₀ and displacement vector at the final step in the four typical specimens with different failure modes, i.e., specimens HC90, HD0, HC30, and HD45. In general, the damage localization zones by numerical simulation coincide with the fractures observed in the experiment (see Figure 13) very well. The distribution of the displacement vector in specimen HC30 reveals the departure movement along the stepped failure plane, and that in specimen HD 45 reveals a sliding of the upper block along the diagonal failure plane.

Figure 30 shows the displacement vector around the hole at the final step in the three typical specimens with different collapse modes, i.e., specimens HC90, HD0, and HC60. The opposite horizontal displacement component at the roof in specimen HC90 or floor in specimen HD0 implies normal departure movement along the tensile damage zones. Compressive shear movement occurs in damage zones at the two sides of the hole in specimens HC90 and HD0 or in the areas between the hole surface and the nearest inclined joints in specimen HC60. The blocks moving into the hole marked by the rectangle areas agree with those observed in the laboratory test.

5. Conclusions

In this study, the combined influence of joint persistence and orientation on the mechanical behavior of a jointed rock mass with a hole under uniaxial compression has been investigated by physical modeling and FE analysis based on a microplane damage model for rocks.

The peak strength σ_p and Young’s modulus E may decrease with the increase in joint continuity factor k, and variation of σ_p and E with joint inclination angle β are the W-shaped or V-shaped curves with the minima and maxima occurring at β = 30° and 90°, respectively. Compared with specimens of failure mode I-axial cleavage or II-crushing, specimens of the failure mode III-stepped or IV-sliding may have a lower peak strength due to involving of shear crack coalescence and sliding failure, which occurs in specimens with inclined joints (β = 30°, 45°, and 60°) at the medium or the largest joint continuity factor (k = 0.5 or 0.8). For most of the specimens, both tensile cracks and crush zones may initiate at the roof/floor and the left/right sides of the hole, respectively, and their linkage with tensile or shear cracks starting from the joints or matrix may lead to the collapse of the hole by three different modes, namely, Mode I_h—without removable blocks, Mode II_h—with removable blocks formed by the upper joints, and Mode III_h—with removable blocks formed by the nearest inclined joints.

The numerical simulation can not only reproduce the main features of macroscopic mechanical behavior observed in physical modeling but also explain the failure process, such as stress concentration, damage evolution, and stress loosening, as well as the interaction mechanisms between the hole and the joints. The strong interaction between the hole and the nearest joints was revealed. At the elastic stage, stress concentration around the hole will be altered by the presence of the joints, and the effect may be strengthened with the increase in the joint continuity factor k. Vertical joints (β = 90°) will have the least influence on stress concentration around the hole, horizontal joints (β = 0°) will induce the largest compressive stress concentration at the two sides of the hole, and inclined joints (β = 30°, 45 and 60°) may cause salient deflection of compressive stress concentration areas around the hole. At the peak strength, stress loosening may happen much earlier before peak strength due to the damage evolution, and the current stress concentration areas will be transferred from the hole surface to the interior. Hoop stress concentration around the hole can be related to the strength of the specimen, i.e., a higher hoop stress concentration factor may lead to lower strength.