# Study of SCA-Induced Rock Crack Propagation under Different Stress Conditions Using a Modified Cohesive Element Method

Authors: Shen Wang; Hani Mitri; Huamin Li; Dongyin Li; Wen Wang
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

When inducing cracks, soundless cracking agents (SCAs) do not generate vibration, harmful gas, dust, nor flying rock fragment, making them suitable for hard rock roof breaking, rock burst prevention, oil or gas reservoir stimulation, and building demolition. In this study, SCA-induced crack initiation and propagation in different stress conditions were modelled using a modified cohesive element method. A new traction-separation law for describing rock compressive shear strength was proposed. The crack length and direction in bidirectional isobaric and unequal stress fields were analyzed in detail. The crack initiation pressure and the incremental ratio of crack length to SCA expansion pressure were proposed as two indicators to evaluate the difficulty in rock breaking in deep underground. Results indicate that (1) the modified cohesive element method used in this study is feasible to model crack propagation in deep rocks; (2) the maximum expansion pressure of SCAs depends on rock elastic modulus and geostress field and should be measured under a condition similar to deep underground prior to SCA borehole spacing design; when using the SCAs with a maximum expansion pressure of 100 MPa in 600 m underground, the suggested borehole spacing is less than 220 mm; and (3)σ3 dominates the crack initiation pressure while the principal stress ratio σ3/σ1 and notch direction control the direction of crack propagation.


---

## Body

## 1. Introduction

The morphology, direction, and propagation length of artificially induced cracks are the important parameters for rock or concrete fracturing projects [1–3]. The common methods of rock and concrete breaking include rock blasting [4, 5], diamond wire cutting [6], controlled blasting [7], hydraulic fracturing [8], and SCAs [9]. SCAs are hopeful to be used in hard rock roof breaking in coal mines (Figure 1(a)), rock burst prevention in metal mines, oil or gas reservoir stimulation (Figure 1(b)) [10], and building demolition [11, 12]. In addition, SCAs are good prospects for rock breaking in subsea tunnel projects as they do not generate vibration, harmful gas, dust, and flying rock fragments when inducing cracks [11, 13].Figure 1
(a) Schematic of SCA application in hard rock breaking in coal mine. (b) CT scanning images of SCA-induced fractures in shale [10].
(a)
(b)SCAs primarily consist of calcium oxide (CaO) and some other additives such as magnesium (MgO), aluminum oxide (Al2O3), and silicon dioxide (SiO2). CaO is the main constituent of SCA. When adding water into SCAs, the hydration reaction between water and CaO generates calcium hydroxide (Ca(OH)2) and expands in volume. The main role of the additives is to control the hydration process. In general, SCAs are filled into boreholes to expand and push the borehole wall [14, 15]. The expansion induces radial compressive stress and tangential tensile stress in the surrounding rock. It is commonly believed that tangential tensile stress is the main factor of crack initiation and propagation. Thus, different multihole arrangement patterns create different rock crack morphologies to meet various engineering requirements.SCAs originated in the 1970s, and since then, many researchers have studied its rock breaking mechanisms. SCA-induced stress distributions under different borehole sizes, rock strengths, and temperatures were studied by Laefer et al. [16] and Mateus and Araújo [17] independently. Dar et al. [18] proposed the optimal spacing of the SCA borehole for breaking marble. De Silva et al. [19] summarized SCA utilities in actual rock breaking projects and put forward its feasible usage in oil formation stimulation. Hanif and Al-Maghrabi [20] developed a SCA application for the granite and marble processing industry. Arshadnejad et al. [21] investigated the crack propagation pathway induced by multihole filled with SCAs. Zhongzhe et al. [22] and Jana [23] proposed a method to determine optimal SCA borehole spacing. Tang et al. [11] analyzed SCA-induced crack initiation and propagation in heterogeneous materials using the RFPA software and proposed that SCA expansion pressure depends on stiffness and confining pressure of the surrounding rock. Tang et al. [24] investigated its feasibility to improve coal permeability. Temperature effect on SCA expansion pressure was also studied numerous times [15, 25–27]. For deep rock projects, the initiation, direction, and length of SCA-induced cracks are influenced by stress field. However, very few studies focused on the effects of geostress field on SCA-induced crack initiation and propagation.Crack propagation modelling issues in rock have been addressed by many studies. Many numerical methods have been performed to model fracture propagation in rock, such as the solid element deletion method [28], the extended finite element method (XFEM) [29], the virtual crack closure technique (VCCT) [30], the cohesive element method [31], the discrete element method (DEM) [32], the smoothed particle hydrodynamics method (SPH) [33], and the discontinuous deformation analysis (DDA) method [34]. The solid element deletion method, as its name suggests, simulates crack propagation by deleting solid elements, so it cannot accurately model contact and slip behaviors between crack surfaces. XFEM is usually used to model brittle crack propagation, but it is difficult to simulate coalescence and branching between a large number of cracks in rock mass. When using VCCT, the location of initial crack is usually required, and besides, it is only able to simulate crack propagation along a prespecified cracking path. Therefore, VCCT is not suitable for modelling rock fracturing. The cohesive element method avoids stress singularity in crack tip by adopting a nonlinear damage zone around the crack tip, resulting in good numerical convergence. Compared to the other methods mentioned above, the cohesive element method can also simulate complex contacts and interactions between cracks after rock fracturing.This study aims to investigate the characteristics of SCA-induced crack propagation in rock under different stress conditions. As shown in Figure2, the notch direction is a significant parameter for controlling the location of crack initiation and direction of crack propagation, so the notch direction was also considered in this study. This paper developed a modified cohesive element method that uses rock compressive shear strength to model SCA-induced crack initiation and propagation. This study is of significance to deepen the understandings of the characteristics and mechanism of SCA-induced crack propagation in deep underground.Figure 2
Schematic of notches around the SCA borehole.


## 2. Methodology for Modelling Crack Propagation in Rock

### 2.1. Cohesive Element Method

Based on the cohesive element method, crack initiation and propagation can be simulated using the finite element method. As shown in Figure3(a), the 2D zero-thickness cohesive element has 4 nodes and 2 integration points. Its geometric thickness is zero, but its constitutive thickness, used to calculate its deformation stiffness, is not equal to zero. The traction forces on the top and bottom faces induce the relative separation displacement between them, and the relative separation displacement components in local 1 direction and 2 direction (Figure 3(a)) represent crack opening and slipping, respectively.Figure 3
(a) Schematic of zero-thickness cohesive element. (b) Process of generating zero-thickness cohesive elements at every solid element interface.
(a)
(b)Figure3(b) illustrates the process of generating zero-thickness cohesive elements in a solid element mesh. The algorithm used was written in the Python programming language and implemented as an ABAQUS plugin. By executing the plugin, cohesive elements are generated at every solid element interface; therefore, the final model consists of solid elements (the yellow elements) and cohesive elements (the blue elements). A solid element and its neighboring cohesive element share two common nodes, and the two adjacent cohesive elements share one node. Crack initiation and propagation are represented by the damage and failure of cohesive elements. Rock elastic-plastic deformation is dominated by solid elements.The mesh should be fine and irregular to avoid mesh structure effect on crack propagation paths and to improve crack propagation randomness. In this study, the irregular mesh was generated using the advancing front technique. The final mesh model is shown in Section4.3.


### 2.2. Initial Stiffness of Cohesive Element

As shown in Figure4(a), when two solid elements are in series, their total stiffness Ko can be calculated by(1)Ko=11/Ks+1/Ks=Ks2,where Ks is the initial stiffness of the solid element.Figure 4
Schematic of (a) two solid elements in series and (b) two solid elements and one cohesive element in series.
(a)
(b)When embedding one cohesive element between two solid elements, the new total stiffnessKo′ can be calculated by(2)Ko′=11/Ks+1/Kc+1/Ks=Ks2+Ks/Kc,where Kc is the initial stiffness of the cohesive element.According to Equations (1) and (2), the total stiffness of the finite element model is affected by embedding cohesive elements. To avoid this effect, Ks must be much greater than Kc in Equation (2).The relationship between element stiffnessK and elastic modulus E is(3)K=Eh,where h is the constitutive element thickness. For a solid element, h is the same as its geometric thickness. For a zero-thickness cohesive element, h is not equal to its geometric thickness (that is, zero). Equation (3) indicates that the constitutive thickness of a cohesive element must be much smaller than the geometric thickness of a solid element to not influence the global mesh stiffness. In this study, the elastic modulus of cohesive element and solid element are the same, and the initial constitutive thickness of cohesive element hoc is one thousandths of the constitutive thickness of solid element hs.


## 2.1. Cohesive Element Method

Based on the cohesive element method, crack initiation and propagation can be simulated using the finite element method. As shown in Figure3(a), the 2D zero-thickness cohesive element has 4 nodes and 2 integration points. Its geometric thickness is zero, but its constitutive thickness, used to calculate its deformation stiffness, is not equal to zero. The traction forces on the top and bottom faces induce the relative separation displacement between them, and the relative separation displacement components in local 1 direction and 2 direction (Figure 3(a)) represent crack opening and slipping, respectively.Figure 3
(a) Schematic of zero-thickness cohesive element. (b) Process of generating zero-thickness cohesive elements at every solid element interface.
(a)
(b)Figure3(b) illustrates the process of generating zero-thickness cohesive elements in a solid element mesh. The algorithm used was written in the Python programming language and implemented as an ABAQUS plugin. By executing the plugin, cohesive elements are generated at every solid element interface; therefore, the final model consists of solid elements (the yellow elements) and cohesive elements (the blue elements). A solid element and its neighboring cohesive element share two common nodes, and the two adjacent cohesive elements share one node. Crack initiation and propagation are represented by the damage and failure of cohesive elements. Rock elastic-plastic deformation is dominated by solid elements.The mesh should be fine and irregular to avoid mesh structure effect on crack propagation paths and to improve crack propagation randomness. In this study, the irregular mesh was generated using the advancing front technique. The final mesh model is shown in Section4.3.


## 2.2. Initial Stiffness of Cohesive Element

As shown in Figure4(a), when two solid elements are in series, their total stiffness Ko can be calculated by(1)Ko=11/Ks+1/Ks=Ks2,where Ks is the initial stiffness of the solid element.Figure 4
Schematic of (a) two solid elements in series and (b) two solid elements and one cohesive element in series.
(a)
(b)When embedding one cohesive element between two solid elements, the new total stiffnessKo′ can be calculated by(2)Ko′=11/Ks+1/Kc+1/Ks=Ks2+Ks/Kc,where Kc is the initial stiffness of the cohesive element.According to Equations (1) and (2), the total stiffness of the finite element model is affected by embedding cohesive elements. To avoid this effect, Ks must be much greater than Kc in Equation (2).The relationship between element stiffnessK and elastic modulus E is(3)K=Eh,where h is the constitutive element thickness. For a solid element, h is the same as its geometric thickness. For a zero-thickness cohesive element, h is not equal to its geometric thickness (that is, zero). Equation (3) indicates that the constitutive thickness of a cohesive element must be much smaller than the geometric thickness of a solid element to not influence the global mesh stiffness. In this study, the elastic modulus of cohesive element and solid element are the same, and the initial constitutive thickness of cohesive element hoc is one thousandths of the constitutive thickness of solid element hs.


## 3. Damage Initiation and Evolution of Cohesive Element

The damage initiation and damage evolution of cohesive elements are described by the traction-separation law [31]. Crack opening and slipping are caused by tractions on the top and bottom faces of cohesive elements including tensile and shear stresses. In this study, tensile stress and tensile displacement are represented by positive values, whereas compressive stress and compressive displacement are flagged as negative values. As shown in Figure 5, once the traction exceeds the strength, the cohesive element begins to damage, and its stiffness decreases as separation displacement increases. When the cohesive element is completely damaged, its stiffness reduces to zero, indicating forward crack propagation. The traction-separation law can be divided into three stages: elastic stage, damage initiation stage, and damage evolution stage.Figure 5
Schematic of the cohesive zone around the crack tip and the traction-separation law for the cohesive element.


### 3.1. Linear Elastic Traction-Separation Equation

It is assumed that traction-separation behavior in the elastic stage obeys the linear elastic law which is the relationship between nominal stresses and nominal strains. Nominal stresses are force components divided by original area at each integration point of a cohesive element, whereas nominal strains are separations divided by initial constitutive thickness, given by(4)εnεs=1hocδnδs,where εn and εs are the nominal tensile and shear strain components, respectively, and δn and δs are the corresponding separation displacement components, respectively.The nominal traction stress vectort consists of two components in 2D problems: tn and ts, which represent normal traction (along the local 2 direction in Figure 3(a)) and shear traction (along the local 1 direction in Figure 3(a)), respectively. The elastic relationship between nominal stresses and nominal strains is represented as(5)t=tnts=En00Esεnεs,where En and Es are the tensile and shear elastic moduli, respectively.


### 3.2. Damage Initiation

In deep underground surrounding a SCA borehole, stress field is affected by the SCA’s expansion pressure and crustal stress field. Rock shear strength usually increases with confining pressure; therefore, compressive shear strength characteristics are crucial for crack propagation modelling.


#### 3.2.1. Compressive Shear Damage Criterion (tn < 0)

Whentn < 0, the damage initiation of cohesive elements is assumed to obey the Mohr–Coulomb criterion. The damage criterion function is given by(6)ts=S−tntanϕ,tn<0,where S is the pure shear strength (i.e., the cohesion) of rock and φ is the friction angle. In the tn-ts plane, Equation (6) is plotted as two radials, as illustrated in Figure 6.Figure 6
Damage criterion in thetn-ts plane.


#### 3.2.2. Tensile-Shear Damage Criterion (tn ≥ 0)

The cohesive element begins to damage when the nominal stress ratios in the quadratic interaction function below reaches 1. This criterion is represented as(7)tnT2+tsS2=1,tn≥0,where T is the tensile strength of cohesive elements. For rock materials, T is usually smaller than S. Therefore, Equation (7) represents a semiellipse in the tn-ts plane, as shown in Figure 6. Equation (7) describes the mixed mode of tensile and shear stresses that causes damage initiation.


### 3.3. Damage Evolution

It is assumed that a linear softening evolution appears after damage initiation. A scalar damage variable,D, represents the overall damage in one cohesive element and captures the combined effects of all the active damage mechanisms. D is zero initially and then gradually evolves from 0 to 1 upon further loading. The stress components of the traction-separation law are affected by the damage according to the following equation:(8)tn=1−Dt¯n,t¯n≥0,tn=t¯n,t¯n<0,ts=1−Dt¯s,where t¯n and t¯s are the stress components predicted by the elastic law at the current strains, as illustrated in Figure 7. Note that the tensile stiffness is degraded when t¯n≥0 but unchanged when t¯n<0.Figure 7
The concept of damage variableD.The mixed mode happens when cohesive elements undergo normal and shear tractions simultaneously. To describe damage evolution of this mode, an effective displacementδm [35] is defined as(9)δm=δn2+δs2,where the operator  is defined as(10)λ=λ,λ≥0,0,λ<0.Therefore,δm is equal to δs when δn<0 (i.e., tn<0).According to Equations (4) and (5), tractions at the damage initiation point are calculated by(11)tno=Enhocδno,tso=Eshocδso,where tno and tso are the normal and shear stress components at the damage initiation point, respectively, and δno and δso are normal and shear displacement components at the damage initiation point, respectively (Figure 7).The mode mixity ratioβ is defined as(12)β=δsδn.According to Equation (9), the effective displacement at the damage initiation point, δmo, is given by(13)δmo=δno2+δso2.According to Equations (12) and (13), δno and δso can be represented by β and δmo as(14)δno=δmo1+β2,δso=βδmo1+β2,where δmo can be obtained by substituting Equations (11) and (14) into (6) and (7), and solving for δmo gives(15)δmo=hocTS1+β2En2S2+β2Es2T2,tn≥0hocS1+β2βEs+Entanϕ,tn<0.Equation (15) indicates that δmo changes with β, which is indirectly determined from the damage initiation criteria.For cohesive elements under tensile-shear loading, the most used damage evolution law is the power law criterion [35], given by(16)GIGICα+GIIGIICα=1,tn≥0,where GI and GII are tension and shear energy components, respectively, and GIC and GIIC are I-mode and II-mode fracture energies, respectively. When GI and GII satisfy Equation (16), the damage variable D reaches 1. GI and GII are calculated as(17)GI=∫0δnftndδn,GII=∫0δsftsdδs,where δnf and δsf are the tensile and shear displacement components at the complete damage point.G
IC and GIIC are the area of triangle O-tno-δnf and triangle O-tso-δsf, respectively (Figure 7), and are calculated by(18)GIC=Enδnoδnf2hoc,GIIC=Esδsoδsf2hoc.According to Equations (9) and (12), δnf and δsf can be represented by β and δmf:(19)δnf=δmf1+β2,δsf=βδmf1+β2,where δmf is the effective displacement at the complete damage point. When α = 1, δmf is obtained by substituting Equations (18) and (19) into (16), given by(20)δmf=2hocGICGIIC1+β2δmoEnGIIC+β2EsGIC.Under compressive shear loading, according to Equation (12), β tends to infinity. Thus, δmf can be solved by letting β → ∞ in Equation (20), which is given by(21)δmf=limβ⟶∞2hocGICGIIC1+β2δmoEnGIIC+β2EsGIC=limβ⟶∞2hocGICGIICδmoEnGIIC+β2EsGIC+limβ⟶∞2hocGICGIICβ2δmoEnGIIC+β2δmoEsGIC=0+2hocGICGIICδmoEsGIC=2hocGIICδmoEs.For the linear softening evolution,D is expressed as(22)D=δmfδmmax−δmoδmmaxδmf−δmo,where δmmax is the maximum effective displacement value attained during loading. As shown in Figure 8, the loading history refers to the loading path (Path 1) or the reloading path (Path 2).Figure 8
Schematic of the concept ofδmmax.D is calculated using Equation (22), and then the stresses are updated by substituting D into Equation (8). Figure 9 illustrates the traction-separation law when considering compressive shear strength characteristics. The damage initiation and evolution laws were written in the Fortran programming language as a UMAT subroutine.Figure 9
Illustration of the whole traction-separation law considering rock’s compressive shear strength characteristics.


## 3.1. Linear Elastic Traction-Separation Equation

It is assumed that traction-separation behavior in the elastic stage obeys the linear elastic law which is the relationship between nominal stresses and nominal strains. Nominal stresses are force components divided by original area at each integration point of a cohesive element, whereas nominal strains are separations divided by initial constitutive thickness, given by(4)εnεs=1hocδnδs,where εn and εs are the nominal tensile and shear strain components, respectively, and δn and δs are the corresponding separation displacement components, respectively.The nominal traction stress vectort consists of two components in 2D problems: tn and ts, which represent normal traction (along the local 2 direction in Figure 3(a)) and shear traction (along the local 1 direction in Figure 3(a)), respectively. The elastic relationship between nominal stresses and nominal strains is represented as(5)t=tnts=En00Esεnεs,where En and Es are the tensile and shear elastic moduli, respectively.


## 3.2. Damage Initiation

In deep underground surrounding a SCA borehole, stress field is affected by the SCA’s expansion pressure and crustal stress field. Rock shear strength usually increases with confining pressure; therefore, compressive shear strength characteristics are crucial for crack propagation modelling.


### 3.2.1. Compressive Shear Damage Criterion (tn < 0)

Whentn < 0, the damage initiation of cohesive elements is assumed to obey the Mohr–Coulomb criterion. The damage criterion function is given by(6)ts=S−tntanϕ,tn<0,where S is the pure shear strength (i.e., the cohesion) of rock and φ is the friction angle. In the tn-ts plane, Equation (6) is plotted as two radials, as illustrated in Figure 6.Figure 6
Damage criterion in thetn-ts plane.


### 3.2.2. Tensile-Shear Damage Criterion (tn ≥ 0)

The cohesive element begins to damage when the nominal stress ratios in the quadratic interaction function below reaches 1. This criterion is represented as(7)tnT2+tsS2=1,tn≥0,where T is the tensile strength of cohesive elements. For rock materials, T is usually smaller than S. Therefore, Equation (7) represents a semiellipse in the tn-ts plane, as shown in Figure 6. Equation (7) describes the mixed mode of tensile and shear stresses that causes damage initiation.


## 3.2.1. Compressive Shear Damage Criterion (tn < 0)

Whentn < 0, the damage initiation of cohesive elements is assumed to obey the Mohr–Coulomb criterion. The damage criterion function is given by(6)ts=S−tntanϕ,tn<0,where S is the pure shear strength (i.e., the cohesion) of rock and φ is the friction angle. In the tn-ts plane, Equation (6) is plotted as two radials, as illustrated in Figure 6.Figure 6
Damage criterion in thetn-ts plane.


## 3.2.2. Tensile-Shear Damage Criterion (tn ≥ 0)

The cohesive element begins to damage when the nominal stress ratios in the quadratic interaction function below reaches 1. This criterion is represented as(7)tnT2+tsS2=1,tn≥0,where T is the tensile strength of cohesive elements. For rock materials, T is usually smaller than S. Therefore, Equation (7) represents a semiellipse in the tn-ts plane, as shown in Figure 6. Equation (7) describes the mixed mode of tensile and shear stresses that causes damage initiation.


## 3.3. Damage Evolution

It is assumed that a linear softening evolution appears after damage initiation. A scalar damage variable,D, represents the overall damage in one cohesive element and captures the combined effects of all the active damage mechanisms. D is zero initially and then gradually evolves from 0 to 1 upon further loading. The stress components of the traction-separation law are affected by the damage according to the following equation:(8)tn=1−Dt¯n,t¯n≥0,tn=t¯n,t¯n<0,ts=1−Dt¯s,where t¯n and t¯s are the stress components predicted by the elastic law at the current strains, as illustrated in Figure 7. Note that the tensile stiffness is degraded when t¯n≥0 but unchanged when t¯n<0.Figure 7
The concept of damage variableD.The mixed mode happens when cohesive elements undergo normal and shear tractions simultaneously. To describe damage evolution of this mode, an effective displacementδm [35] is defined as(9)δm=δn2+δs2,where the operator  is defined as(10)λ=λ,λ≥0,0,λ<0.Therefore,δm is equal to δs when δn<0 (i.e., tn<0).According to Equations (4) and (5), tractions at the damage initiation point are calculated by(11)tno=Enhocδno,tso=Eshocδso,where tno and tso are the normal and shear stress components at the damage initiation point, respectively, and δno and δso are normal and shear displacement components at the damage initiation point, respectively (Figure 7).The mode mixity ratioβ is defined as(12)β=δsδn.According to Equation (9), the effective displacement at the damage initiation point, δmo, is given by(13)δmo=δno2+δso2.According to Equations (12) and (13), δno and δso can be represented by β and δmo as(14)δno=δmo1+β2,δso=βδmo1+β2,where δmo can be obtained by substituting Equations (11) and (14) into (6) and (7), and solving for δmo gives(15)δmo=hocTS1+β2En2S2+β2Es2T2,tn≥0hocS1+β2βEs+Entanϕ,tn<0.Equation (15) indicates that δmo changes with β, which is indirectly determined from the damage initiation criteria.For cohesive elements under tensile-shear loading, the most used damage evolution law is the power law criterion [35], given by(16)GIGICα+GIIGIICα=1,tn≥0,where GI and GII are tension and shear energy components, respectively, and GIC and GIIC are I-mode and II-mode fracture energies, respectively. When GI and GII satisfy Equation (16), the damage variable D reaches 1. GI and GII are calculated as(17)GI=∫0δnftndδn,GII=∫0δsftsdδs,where δnf and δsf are the tensile and shear displacement components at the complete damage point.G
IC and GIIC are the area of triangle O-tno-δnf and triangle O-tso-δsf, respectively (Figure 7), and are calculated by(18)GIC=Enδnoδnf2hoc,GIIC=Esδsoδsf2hoc.According to Equations (9) and (12), δnf and δsf can be represented by β and δmf:(19)δnf=δmf1+β2,δsf=βδmf1+β2,where δmf is the effective displacement at the complete damage point. When α = 1, δmf is obtained by substituting Equations (18) and (19) into (16), given by(20)δmf=2hocGICGIIC1+β2δmoEnGIIC+β2EsGIC.Under compressive shear loading, according to Equation (12), β tends to infinity. Thus, δmf can be solved by letting β → ∞ in Equation (20), which is given by(21)δmf=limβ⟶∞2hocGICGIIC1+β2δmoEnGIIC+β2EsGIC=limβ⟶∞2hocGICGIICδmoEnGIIC+β2EsGIC+limβ⟶∞2hocGICGIICβ2δmoEnGIIC+β2δmoEsGIC=0+2hocGICGIICδmoEsGIC=2hocGIICδmoEs.For the linear softening evolution,D is expressed as(22)D=δmfδmmax−δmoδmmaxδmf−δmo,where δmmax is the maximum effective displacement value attained during loading. As shown in Figure 8, the loading history refers to the loading path (Path 1) or the reloading path (Path 2).Figure 8
Schematic of the concept ofδmmax.D is calculated using Equation (22), and then the stresses are updated by substituting D into Equation (8). Figure 9 illustrates the traction-separation law when considering compressive shear strength characteristics. The damage initiation and evolution laws were written in the Fortran programming language as a UMAT subroutine.Figure 9
Illustration of the whole traction-separation law considering rock’s compressive shear strength characteristics.


## 4. Numerical Simulation Procedure

### 4.1. Mechanical Properties of Rock

The mechanical properties used in this study were obtained from laboratory tests. The rock samples are medium-grained sandstone obtained from the roof of coal seam No. 22 in the Burtai coal mine in the Ordos Basin, China. The elastic modulus and Poisson’s ratio were measured by uniaxial compression tests, the tensile strength was obtained from Brazilian split tests, the cohesion and friction angle were tested by triaxial compression tests, and the mode-I and mode-II fracture energies were determined by three-point bending tests and four-point shear tests, respectively. Table1 lists the mechanical properties. Note that the elastic modulus Eo, cohesion c, friction angle φo, and tensile strength st of the solid elements are equal to the elastic modulus En, pure shear strength S, friction angle φ, and tensile strength T of the cohesive elements.Table 1
Rock mechanical properties used in this study.
Mechanical parameters
Values
Elastic modulus of the solid element,Eo (GPa)
10.16
Poisson’s ratio of the solid element,v
0.28
Cohesion of the solid element,c (MPa)
18.3
Friction angle of the solid element,φo (°)
35.0
Tensile strength,st (MPa)
3.97
Elastic modulus of the cohesive element,En (GPa)
10.16
Shear modulus of the cohesive element,Es (GPa)
3.91
Pure shear strength of the cohesive element,S (MPa)
18.3
Tensile strength of the cohesive element,T (MPa)
3.97
Equivalent friction angle of the cohesive element,φ (°)
35.0
Mode-I fracture energy,GIC (Pa·m)
103
Mode-II fracture energy,GIIC (Pa·m)
154


### 4.2. SCA Expansion Property

The variation of SCA expansion pressure over time can be measured in a confined state. As shown in Figure10, SCAs are put into a thick-walled steel cylinder tube, and then the SCA expansion property will be reflected on the outer surface deformation of the cylinder tube. According to the literature [9] and [24], for any point (ρ, θ, and z) in the cylinder tube, its stress state is calculated by(23)σρ=r2ρ2−R2ρ2R2−r2Pin−R2ρ2−r2ρ2R2−r2Pex,σθ=r2ρ2+R2ρ2R2−r2Pin−R2ρ2+r2ρ2R2−r2Pex,τρθ=0,σz=Pin,where σρ, σθ, and σz are the stress components along the radical, tangential, and axial directions, respectively; R and r are the inner and outer radii of the steel cylinder tube; and Pin and Pex are the pressures of the inner surface and outside surface.Figure 10
Schematic of how SCA expansion pressure is measured in a thick-walled steel cylinder tube.According to Equation (23), the stress components on the outside surface of the steel cylinder tube are calculated by(24)σρρ=R=−Pex,σθρ=R=2r2R2−r2Pin−R2+r2R2−r2Pex,τρθρ=R=0,σzρ=R=Pin.According to the generalized Hooke law, the tangential strain component is(25)εθ=σθ−νσρEp,where Ep is the elastic modulus of the surrounding materials.By substituting Equation (24) into Equation (25), Pin can be obtained as(26)Pin=R2−r22r2εθEp+R2+r2R2−r2−vPex.Equation (26) shows that Pin is influenced by both Ep and Pex. Borehole expansion pressure changes with Pex, so the SCA expansion property measured under Pex = 0 is not suitable for numerical simulation when ground stress field is taken into account.According to Tang et al. [11] and De Silva et al. [9], the SCA expansion process is assumed to be a static process because SCAs usually take 20 h or longer to reach its maximum expansion pressure. Therefore, the loading rate of expansion pressure has very few effects on the mechanical response of its surrounding materials. In this study, the expansion pressure was applied to the borehole wall using a static loading process from 0 MPa to 150 MPa.


### 4.3. Mesh Model

In fact, SCA boreholes must have an opening that allows it to expand freely along the axial direction. The free expansion behavior along the axial direction is not what mainly causes rock fracturing, so the interaction between SCAs and surrounding rock along the axial direction was neglected in this study. Additionally, in order to reduce computation complexity in ABAQUS, a 2D model was used, which is equivalent to a 3D model by setting an out-of-plane dimension. As shown in Figure11, the mesh is a 500 mm × 500 mm square. The SCA borehole diameter is 30 mm, and the distance between the two notch tips is 56 mm. The mesh, initially made of 193413 triangular elements, was generated using the advancing front technique, and then 273441 zero-thickness cohesive elements were embedded into the mesh.Figure 11
The mesh model.The bottom of the model was fixed in the vertical direction. The maximum principal stressσ1 and the minimum principal stress σ3 were applied in both the vertical and horizontal directions, respectively. The numerical simulation consists of two steps: first, σ1 and σ3 were applied to the corresponding model boundaries, and the geostatic stress was equilibrated. Then, SCA expansion pressure was gradually applied to the borehole wall from 0 MPa to 150 MPa.


### 4.4. Test Schemes

A series of numerical tests were carried out to study the characteristics of SCA-induced crack initiation and propagation under different stress conditions and notch directions. Stress conditions were categorized intoσ1 = σ3 = 0 MPa (free boundary conditions where the rock sample is not fixed and not loaded with any external forces), σ1 = σ3 ≠ 0 MPa (bidirectional isobaric stress field), and σ1−σ3 = 5 MPa (bidirectional unequal stress field), as shown in Table 2, where θ is the angle between notch direction and σ1 direction (Figure 2).Table 2
The stress loading conditions (unit: MPa).
Notch conditions
Loading conditions
(σ1, σ3)
No notches
σ
1 = σ3
(0,0)
θ = 0°, 45° and 60°
σ
1 = σ3
(5,5)
(7.5,7.5)
(10,10)
(12.5,12.5)
(15,15)
(17.5,17.5)
(20,20)
θ = 0°, 45° and 90°
σ
1 − σ3 = 5
(5,0)
(7.5,2.5)
(10,5)
(12.5,7.5)
(15,10)
(17.5,12.5)
(20,15)


## 4.1. Mechanical Properties of Rock

The mechanical properties used in this study were obtained from laboratory tests. The rock samples are medium-grained sandstone obtained from the roof of coal seam No. 22 in the Burtai coal mine in the Ordos Basin, China. The elastic modulus and Poisson’s ratio were measured by uniaxial compression tests, the tensile strength was obtained from Brazilian split tests, the cohesion and friction angle were tested by triaxial compression tests, and the mode-I and mode-II fracture energies were determined by three-point bending tests and four-point shear tests, respectively. Table1 lists the mechanical properties. Note that the elastic modulus Eo, cohesion c, friction angle φo, and tensile strength st of the solid elements are equal to the elastic modulus En, pure shear strength S, friction angle φ, and tensile strength T of the cohesive elements.Table 1
Rock mechanical properties used in this study.
Mechanical parameters
Values
Elastic modulus of the solid element,Eo (GPa)
10.16
Poisson’s ratio of the solid element,v
0.28
Cohesion of the solid element,c (MPa)
18.3
Friction angle of the solid element,φo (°)
35.0
Tensile strength,st (MPa)
3.97
Elastic modulus of the cohesive element,En (GPa)
10.16
Shear modulus of the cohesive element,Es (GPa)
3.91
Pure shear strength of the cohesive element,S (MPa)
18.3
Tensile strength of the cohesive element,T (MPa)
3.97
Equivalent friction angle of the cohesive element,φ (°)
35.0
Mode-I fracture energy,GIC (Pa·m)
103
Mode-II fracture energy,GIIC (Pa·m)
154


## 4.2. SCA Expansion Property

The variation of SCA expansion pressure over time can be measured in a confined state. As shown in Figure10, SCAs are put into a thick-walled steel cylinder tube, and then the SCA expansion property will be reflected on the outer surface deformation of the cylinder tube. According to the literature [9] and [24], for any point (ρ, θ, and z) in the cylinder tube, its stress state is calculated by(23)σρ=r2ρ2−R2ρ2R2−r2Pin−R2ρ2−r2ρ2R2−r2Pex,σθ=r2ρ2+R2ρ2R2−r2Pin−R2ρ2+r2ρ2R2−r2Pex,τρθ=0,σz=Pin,where σρ, σθ, and σz are the stress components along the radical, tangential, and axial directions, respectively; R and r are the inner and outer radii of the steel cylinder tube; and Pin and Pex are the pressures of the inner surface and outside surface.Figure 10
Schematic of how SCA expansion pressure is measured in a thick-walled steel cylinder tube.According to Equation (23), the stress components on the outside surface of the steel cylinder tube are calculated by(24)σρρ=R=−Pex,σθρ=R=2r2R2−r2Pin−R2+r2R2−r2Pex,τρθρ=R=0,σzρ=R=Pin.According to the generalized Hooke law, the tangential strain component is(25)εθ=σθ−νσρEp,where Ep is the elastic modulus of the surrounding materials.By substituting Equation (24) into Equation (25), Pin can be obtained as(26)Pin=R2−r22r2εθEp+R2+r2R2−r2−vPex.Equation (26) shows that Pin is influenced by both Ep and Pex. Borehole expansion pressure changes with Pex, so the SCA expansion property measured under Pex = 0 is not suitable for numerical simulation when ground stress field is taken into account.According to Tang et al. [11] and De Silva et al. [9], the SCA expansion process is assumed to be a static process because SCAs usually take 20 h or longer to reach its maximum expansion pressure. Therefore, the loading rate of expansion pressure has very few effects on the mechanical response of its surrounding materials. In this study, the expansion pressure was applied to the borehole wall using a static loading process from 0 MPa to 150 MPa.


## 4.3. Mesh Model

In fact, SCA boreholes must have an opening that allows it to expand freely along the axial direction. The free expansion behavior along the axial direction is not what mainly causes rock fracturing, so the interaction between SCAs and surrounding rock along the axial direction was neglected in this study. Additionally, in order to reduce computation complexity in ABAQUS, a 2D model was used, which is equivalent to a 3D model by setting an out-of-plane dimension. As shown in Figure11, the mesh is a 500 mm × 500 mm square. The SCA borehole diameter is 30 mm, and the distance between the two notch tips is 56 mm. The mesh, initially made of 193413 triangular elements, was generated using the advancing front technique, and then 273441 zero-thickness cohesive elements were embedded into the mesh.Figure 11
The mesh model.The bottom of the model was fixed in the vertical direction. The maximum principal stressσ1 and the minimum principal stress σ3 were applied in both the vertical and horizontal directions, respectively. The numerical simulation consists of two steps: first, σ1 and σ3 were applied to the corresponding model boundaries, and the geostatic stress was equilibrated. Then, SCA expansion pressure was gradually applied to the borehole wall from 0 MPa to 150 MPa.


## 4.4. Test Schemes

A series of numerical tests were carried out to study the characteristics of SCA-induced crack initiation and propagation under different stress conditions and notch directions. Stress conditions were categorized intoσ1 = σ3 = 0 MPa (free boundary conditions where the rock sample is not fixed and not loaded with any external forces), σ1 = σ3 ≠ 0 MPa (bidirectional isobaric stress field), and σ1−σ3 = 5 MPa (bidirectional unequal stress field), as shown in Table 2, where θ is the angle between notch direction and σ1 direction (Figure 2).Table 2
The stress loading conditions (unit: MPa).
Notch conditions
Loading conditions
(σ1, σ3)
No notches
σ
1 = σ3
(0,0)
θ = 0°, 45° and 60°
σ
1 = σ3
(5,5)
(7.5,7.5)
(10,10)
(12.5,12.5)
(15,15)
(17.5,17.5)
(20,20)
θ = 0°, 45° and 90°
σ
1 − σ3 = 5
(5,0)
(7.5,2.5)
(10,5)
(12.5,7.5)
(15,10)
(17.5,12.5)
(20,15)


## 5. Numerical Simulation Results and Discussion

### 5.1. Crack Characteristics under the Free Boundary Conditions

To predict crack propagation using the cohesive element method, the characteristics of SCA-induced crack propagation under the free boundary conditions were computed first. The numerical results were compared with the laboratory results of Tang et al. [11], as shown in Figure 12. Figures 12(a) and 12(b) show the test results from two different samples by Tang et al. [11]. The cracks inside the samples are similar to those on the surfaces of the intact rock samples, which is confirmed by Tang et al.’s results [24]. This is because in each of their tests [11, 24], at least one end of the SCA borehole was unrestricted, so the SCA expansion pressure along the axial direction was released and unable to affect the intact rock samples. However, if a rock sample contains obvious natural fractures, the SCA-induced cracks inside it may be different from those on its surface. The effect of natural fractures on SCA-induced cracks is not considered in this study and will be addressed in the future. Figure 12 shows that the crack angles are significantly different for the two experiments but similar in the numerical results. The crack angle difference in Figures 12(a) and 12(b) is due to different mineral grain and natural microfracture distribution. The simulations in this study did not focus on the mesoscopic heterogeneity of rock materials. The difference between the results of the simulation and that of the experiment is due to the rock parameters used, but both results show three main cracks in each sample. In ABAQUS, tensile stress is positive. In Figure 12(c), tensile stress concentrates around the crack tips, and the cracks propagate gradually as expansion pressure increases. Generally, the crack propagation characteristics simulated by the cohesive element method is similar to those in the experimental results in the literature [11], which suggests that the cohesive element method based on the modified traction-separation law can model SCA-induced crack initiation and propagation.Figure 12
The SCA-induced crack characteristics of the experiment results and that of the numerical results under the free boundary conditions. (a) and (b) show the experiment results of Tang et al. [10]. (c) shows the numerical results.
(a)
(b)
(c)


### 5.2. Crack Characteristics in Bidirectional Isobaric Stress Fields

The numerical results of crack propagation underσ1 = σ3 = 10 MPa conditions are illustrated in Figure 13. There are three main cracks under the free boundary conditions, but only two main cracks in the bidirectional isobaric stress field where the SCA-induced cracks are initiated at the notch tips and the crack directions are roughly parallel to the notch direction. Crack length gradually increases with the SCA expansion pressure P. According to Figure 13, the extent of tensile stress concentrations around the crack tips also increase as the cracks propagate. At the beginning of crack propagation, the tensile stress concentrations only appear around the crack tips. However, as P increases, the crack openings and borehole diameter expand, causing bulge and tensile stress concentration in the middle of the boundaries.Figure 13
Crack characteristics underσ1 = σ3 = 10 MPa and under different notch conditions. (a) θ = 0°; (b) θ = 45°; (c) θ = 60°. P is the SCA expansion pressure.
(a)
(b)
(c)To apply SCAs in deep underground, the spacing of SCA boreholes, which depends on the crack lengthl, is crucial, so further analysis of the relationship between P and l under different stress conditions should be undertaken. P-l linear fit in different bidirectional isobaric stress fields is plotted in Figure 14, and the corresponding equations are listed in Table 3. According to Figure 14, it is obvious that in each stress field, SCA-induced crack length increases with its expansion pressure. The SCA crack initiation pressure P0 is a proper parameter for evaluating SCAs’ rock breaking ability (the higher the P0, the more difficult it is to break rock). Under the free boundary conditions, P0 is only approximately 4–6 MPa [11]. Under σ1 = σ3 = 5 MPa, P0 is 16.6236 MPa. However, when σ1 = σ3 = 20 MPa, P0 increases to 37.8456 MPa, which is approximately 6.3–9.5 times that under the free boundary conditions.Figure 14
(a) The variations of crack length with increasing SCA expansion pressure under different bidirectional isobaric stress fields. (b) Variations of crack initiation pressure P0 and the slope of P-l linear fits al as confining stress increases.
(a)
(b)Table 3
TheP-l linear fit equations under the conditions of σ1 = σ3.
Confining stress (MPa)
P-l linear fitting equations
AdjustedR2
Expansion pressure of crack initiation,P0 (MPa)
5
l = 11.908P − 197.95429
0.94012
16.6236
7.5
l = 7.09694P − 144.0286
0.95673
20.2945
10
l = 5.33562P − 129.77332
0.96778
24.3221
12.5
l = 4.12966P − 122.09284
0.96812
29.5647
15
l = 3.52583P − 122.94985
0.97831
34.8712
17.5
l = 2.71927P − 99.53406
0.9667
36.6032
20
l = 2.18073P − 82.53101
0.97728
37.8456The slope of theP-l linear fit (unit: mm/MPa) al is the incremental ratio of crack length to SCA expansion pressure. As shown in Figure 14(b), al decreases as confining stress increases. The lower the al, the more difficult it is to break rock. The reduction rate of al gradually slows down as confining stress increases. l is approximately 500 mm under P = 60 MPa and σ1 = σ3 = 5 MPa, but 50 mm under the same P but σ1 = σ3 = 20 MPa. Overall, the results suggest that the bidirectional isobaric stress field does not significantly affect the crack direction, but can significantly affect crack initiation pressure P0 and crack length l.According to Equation (26), the maximum SCA expansion pressure depends on the rock elastic modulus and geostress field, so it should be measured in a condition similar to deep underground before designing borehole spacing. According to Figure 14, at 600 m underground under 100 MPa of maximum expansion pressure (vertical stress of approximately 15 MPa), SCA borehole spacing should be smaller than 220 mm.


### 5.3. Crack Characteristics in Bidirectional Unequal Stress Fields

The crack characteristics under the conditions ofσ1–σ3 = 5 MPa are shown in Figure 15. When the notch direction θ = 0° (Figure 15(a)), the cracks are roughly parallel to the σ1 direction in each stress field. However, when θ = 45°, although the stress difference is just 5 MPa in each 2field, the crack direction is significantly different (Figure 15(b)). When σ1 = 5 MPa, the angle between the crack direction and σ3 direction θc is 82.8°, whereas it is 65.5°, 58.5°, and 53.9° when σ1 = 10 MPa, 15 MPa, and 20 MPa respectively. Figure 16 illustrates the variations of θc with respect to the stress field when θ = 45°. Figure 16(a) shows that the crack direction rotates gradually from the σ1 direction to the notch direction as both σ1 and σ3 increase. In Figure 16(b), it is worth noting that θc decreases linearly as the stress ratio σ3/σ1 increases. When σ3/σ1 approaches 1, θc = 45°, indicating that the crack direction tends toward the notch direction. Conversely, when σ3/σ1 approaches 0, the crack direction tends to the σ1 direction. Their linear relationship is(27)θc=83.61466−38.69622σ3σ1.Figure 15
The crack characteristics under the conditions ofσ1 = σ3 = 10 MPa and of different notch directions. (a) θ = 0°; (b) θ = 45°; (c) θ = 90°.
(a)
(b)
(c)Figure 16
Variation of the SCA-induced crack direction with respect to (a)σ1 and σ3 and (b) σ3/σ1 when θ = 45°.
(a)
(b)Whenθ = 90°, four cracks are initiated, two of which (named cracks N) appear at the notch tips, and the other two (named cracks M) appear at the center of the SCA borehole wall, as shown in Figure 15(c). Whether cracks can continuously propagate depends on the stress field. Cracks N only propagate continuously along the σ1 direction when σ3 = 0; however, as σ1 and σ3 increase, only cracks M can grow parallel to the notch direction. According to the above results under the conditions of θ = 45° and θ = 90°, it can be concluded that the crack propagates roughly along the σ1 direction when σ3 is close to zero; however, as σ3/σ1 increases, both the notch direction and stress field affect the direction of crack propagation.Figure17 shows the relationship of P-l in different bidirectional unequal stress fields and under different notch conditions. Under the condition of θ = 45°, generating 200 mm cracks requires approximately 13 MPa of expansion pressure when (σ1, σ3) = (5 MPa, 0 MPa), compared to 100 MPa when (σ1, σ3) = (20 MPa, 15 MPa). Figures 15 and 17 show that P gradually increases with respect to σ3 when cracks of the same length are generated. Compared with the numerical results under the σ1 = σ3 conditions, it can be concluded that σ3 is the main factor affecting the crack initiation pressure P0, while both the stress ratio σ3/σ1 and notch direction control crack direction.Figure 17
Variations of the crack length versus SCA expansion pressure under different stress and notch direction conditions. (a)θ = 0°; (b) θ = 45°; (c) θ = 90°.
(a)
(b)
(c)The correspondingP-l linear fit equations are listed in Table 4. Figure 18 shows the variations of P0 and al versus σ3 when σ1 − σ3 is kept constant at 5 MPa. According to Figure 18, P0 increases with σ3 under all notch direction conditions; P0 increases with θ under all stress conditions. The SCA expansion pressure required for crack initiation that lowers as the notch direction tends toward the σ1 direction. When θ = 0° and θ = 45°, the trends of P0-σ3 curves are approximately linear, but when θ = 90°, the trend flatlines as σ3 > 12.5 MPa. Under all notch conditions, al is inversely proportional to σ3. The results suggest that crack generation becomes more difficult as σ3 and θ increase.Table 4
TheP-l linear fit equations under the conditions of σ1–σ3 = 5 MPa.
Notch direction,θ (°)
Stress condition (MPa)
P-l linear fitting equations
AdjustedR2
Crack initiation pressure,P0 (MPa)
0
(5,0)
l = 127.54345P − 704.87343
0.80767
5.5265
(7.5,2.5)
l = 20.40528P − 201.65121
0.86493
9.8823
(10,5)
l = 11.65431P − 183.93688
0.94212
15.7827
(12.5,7.5)
l = 7.56828P − 150.47596
0.95267
19.8825
(15,10)
l = 5.61437P − 142.98323
0.96422
25.4674
(17.5,12.5)
l = 4.60247P − 142.83603
0.95935
31.0346
(20,15)
l = 3.70937P − 131.56198
0.97198
35.4675
45
(5,0)
l = 87.80175P − 705.51258
0.86615
8.0353
(7.5,2.5)
l = 17.68538P − 231.72743
0.89638
13.1028
(10,5)
l = 9.49899P − 176.70845
0.96058
18.6029
(12.5,7.5)
l = 6.35449P − 150.48197
0.9683
23.6812
(15,10)
l = 5.23749P − 159.45243
0.98428
30.4444
(17.5,12.5)
l = 4.10701P − 142.96515
0.9775
34.8100
(20,15)
l = 3.35472P − 132.33265
0.96648
39.4467
90
(5,0)
l = 22.25843P − 69.97357
0.31616
3.1437
(7.5,2.5)
l = 7.89522P − 163.5916
0.95204
20.7203
(10,5)
l = 5.82953P − 153.31147
0.95795
26.2991
(12.5,7.5)
l = 4.51084P − 138.39713
0.96233
30.6810
(15,10)
l = 3.53828P − 120.16012
0.96753
33.9600
(17.5,12.5)
l = 2.86423P − 106.77622
0.96722
37.2792
(20,15)
l = 2.23186P − 84.21847
0.98777
37.7346Figure 18
The variations of crack initiation pressureP0 and the slope of P-l linear fit al versus σ3 under different notch direction conditions.


## 5.1. Crack Characteristics under the Free Boundary Conditions

To predict crack propagation using the cohesive element method, the characteristics of SCA-induced crack propagation under the free boundary conditions were computed first. The numerical results were compared with the laboratory results of Tang et al. [11], as shown in Figure 12. Figures 12(a) and 12(b) show the test results from two different samples by Tang et al. [11]. The cracks inside the samples are similar to those on the surfaces of the intact rock samples, which is confirmed by Tang et al.’s results [24]. This is because in each of their tests [11, 24], at least one end of the SCA borehole was unrestricted, so the SCA expansion pressure along the axial direction was released and unable to affect the intact rock samples. However, if a rock sample contains obvious natural fractures, the SCA-induced cracks inside it may be different from those on its surface. The effect of natural fractures on SCA-induced cracks is not considered in this study and will be addressed in the future. Figure 12 shows that the crack angles are significantly different for the two experiments but similar in the numerical results. The crack angle difference in Figures 12(a) and 12(b) is due to different mineral grain and natural microfracture distribution. The simulations in this study did not focus on the mesoscopic heterogeneity of rock materials. The difference between the results of the simulation and that of the experiment is due to the rock parameters used, but both results show three main cracks in each sample. In ABAQUS, tensile stress is positive. In Figure 12(c), tensile stress concentrates around the crack tips, and the cracks propagate gradually as expansion pressure increases. Generally, the crack propagation characteristics simulated by the cohesive element method is similar to those in the experimental results in the literature [11], which suggests that the cohesive element method based on the modified traction-separation law can model SCA-induced crack initiation and propagation.Figure 12
The SCA-induced crack characteristics of the experiment results and that of the numerical results under the free boundary conditions. (a) and (b) show the experiment results of Tang et al. [10]. (c) shows the numerical results.
(a)
(b)
(c)


## 5.2. Crack Characteristics in Bidirectional Isobaric Stress Fields

The numerical results of crack propagation underσ1 = σ3 = 10 MPa conditions are illustrated in Figure 13. There are three main cracks under the free boundary conditions, but only two main cracks in the bidirectional isobaric stress field where the SCA-induced cracks are initiated at the notch tips and the crack directions are roughly parallel to the notch direction. Crack length gradually increases with the SCA expansion pressure P. According to Figure 13, the extent of tensile stress concentrations around the crack tips also increase as the cracks propagate. At the beginning of crack propagation, the tensile stress concentrations only appear around the crack tips. However, as P increases, the crack openings and borehole diameter expand, causing bulge and tensile stress concentration in the middle of the boundaries.Figure 13
Crack characteristics underσ1 = σ3 = 10 MPa and under different notch conditions. (a) θ = 0°; (b) θ = 45°; (c) θ = 60°. P is the SCA expansion pressure.
(a)
(b)
(c)To apply SCAs in deep underground, the spacing of SCA boreholes, which depends on the crack lengthl, is crucial, so further analysis of the relationship between P and l under different stress conditions should be undertaken. P-l linear fit in different bidirectional isobaric stress fields is plotted in Figure 14, and the corresponding equations are listed in Table 3. According to Figure 14, it is obvious that in each stress field, SCA-induced crack length increases with its expansion pressure. The SCA crack initiation pressure P0 is a proper parameter for evaluating SCAs’ rock breaking ability (the higher the P0, the more difficult it is to break rock). Under the free boundary conditions, P0 is only approximately 4–6 MPa [11]. Under σ1 = σ3 = 5 MPa, P0 is 16.6236 MPa. However, when σ1 = σ3 = 20 MPa, P0 increases to 37.8456 MPa, which is approximately 6.3–9.5 times that under the free boundary conditions.Figure 14
(a) The variations of crack length with increasing SCA expansion pressure under different bidirectional isobaric stress fields. (b) Variations of crack initiation pressure P0 and the slope of P-l linear fits al as confining stress increases.
(a)
(b)Table 3
TheP-l linear fit equations under the conditions of σ1 = σ3.
Confining stress (MPa)
P-l linear fitting equations
AdjustedR2
Expansion pressure of crack initiation,P0 (MPa)
5
l = 11.908P − 197.95429
0.94012
16.6236
7.5
l = 7.09694P − 144.0286
0.95673
20.2945
10
l = 5.33562P − 129.77332
0.96778
24.3221
12.5
l = 4.12966P − 122.09284
0.96812
29.5647
15
l = 3.52583P − 122.94985
0.97831
34.8712
17.5
l = 2.71927P − 99.53406
0.9667
36.6032
20
l = 2.18073P − 82.53101
0.97728
37.8456The slope of theP-l linear fit (unit: mm/MPa) al is the incremental ratio of crack length to SCA expansion pressure. As shown in Figure 14(b), al decreases as confining stress increases. The lower the al, the more difficult it is to break rock. The reduction rate of al gradually slows down as confining stress increases. l is approximately 500 mm under P = 60 MPa and σ1 = σ3 = 5 MPa, but 50 mm under the same P but σ1 = σ3 = 20 MPa. Overall, the results suggest that the bidirectional isobaric stress field does not significantly affect the crack direction, but can significantly affect crack initiation pressure P0 and crack length l.According to Equation (26), the maximum SCA expansion pressure depends on the rock elastic modulus and geostress field, so it should be measured in a condition similar to deep underground before designing borehole spacing. According to Figure 14, at 600 m underground under 100 MPa of maximum expansion pressure (vertical stress of approximately 15 MPa), SCA borehole spacing should be smaller than 220 mm.


## 5.3. Crack Characteristics in Bidirectional Unequal Stress Fields

The crack characteristics under the conditions ofσ1–σ3 = 5 MPa are shown in Figure 15. When the notch direction θ = 0° (Figure 15(a)), the cracks are roughly parallel to the σ1 direction in each stress field. However, when θ = 45°, although the stress difference is just 5 MPa in each 2field, the crack direction is significantly different (Figure 15(b)). When σ1 = 5 MPa, the angle between the crack direction and σ3 direction θc is 82.8°, whereas it is 65.5°, 58.5°, and 53.9° when σ1 = 10 MPa, 15 MPa, and 20 MPa respectively. Figure 16 illustrates the variations of θc with respect to the stress field when θ = 45°. Figure 16(a) shows that the crack direction rotates gradually from the σ1 direction to the notch direction as both σ1 and σ3 increase. In Figure 16(b), it is worth noting that θc decreases linearly as the stress ratio σ3/σ1 increases. When σ3/σ1 approaches 1, θc = 45°, indicating that the crack direction tends toward the notch direction. Conversely, when σ3/σ1 approaches 0, the crack direction tends to the σ1 direction. Their linear relationship is(27)θc=83.61466−38.69622σ3σ1.Figure 15
The crack characteristics under the conditions ofσ1 = σ3 = 10 MPa and of different notch directions. (a) θ = 0°; (b) θ = 45°; (c) θ = 90°.
(a)
(b)
(c)Figure 16
Variation of the SCA-induced crack direction with respect to (a)σ1 and σ3 and (b) σ3/σ1 when θ = 45°.
(a)
(b)Whenθ = 90°, four cracks are initiated, two of which (named cracks N) appear at the notch tips, and the other two (named cracks M) appear at the center of the SCA borehole wall, as shown in Figure 15(c). Whether cracks can continuously propagate depends on the stress field. Cracks N only propagate continuously along the σ1 direction when σ3 = 0; however, as σ1 and σ3 increase, only cracks M can grow parallel to the notch direction. According to the above results under the conditions of θ = 45° and θ = 90°, it can be concluded that the crack propagates roughly along the σ1 direction when σ3 is close to zero; however, as σ3/σ1 increases, both the notch direction and stress field affect the direction of crack propagation.Figure17 shows the relationship of P-l in different bidirectional unequal stress fields and under different notch conditions. Under the condition of θ = 45°, generating 200 mm cracks requires approximately 13 MPa of expansion pressure when (σ1, σ3) = (5 MPa, 0 MPa), compared to 100 MPa when (σ1, σ3) = (20 MPa, 15 MPa). Figures 15 and 17 show that P gradually increases with respect to σ3 when cracks of the same length are generated. Compared with the numerical results under the σ1 = σ3 conditions, it can be concluded that σ3 is the main factor affecting the crack initiation pressure P0, while both the stress ratio σ3/σ1 and notch direction control crack direction.Figure 17
Variations of the crack length versus SCA expansion pressure under different stress and notch direction conditions. (a)θ = 0°; (b) θ = 45°; (c) θ = 90°.
(a)
(b)
(c)The correspondingP-l linear fit equations are listed in Table 4. Figure 18 shows the variations of P0 and al versus σ3 when σ1 − σ3 is kept constant at 5 MPa. According to Figure 18, P0 increases with σ3 under all notch direction conditions; P0 increases with θ under all stress conditions. The SCA expansion pressure required for crack initiation that lowers as the notch direction tends toward the σ1 direction. When θ = 0° and θ = 45°, the trends of P0-σ3 curves are approximately linear, but when θ = 90°, the trend flatlines as σ3 > 12.5 MPa. Under all notch conditions, al is inversely proportional to σ3. The results suggest that crack generation becomes more difficult as σ3 and θ increase.Table 4
TheP-l linear fit equations under the conditions of σ1–σ3 = 5 MPa.
Notch direction,θ (°)
Stress condition (MPa)
P-l linear fitting equations
AdjustedR2
Crack initiation pressure,P0 (MPa)
0
(5,0)
l = 127.54345P − 704.87343
0.80767
5.5265
(7.5,2.5)
l = 20.40528P − 201.65121
0.86493
9.8823
(10,5)
l = 11.65431P − 183.93688
0.94212
15.7827
(12.5,7.5)
l = 7.56828P − 150.47596
0.95267
19.8825
(15,10)
l = 5.61437P − 142.98323
0.96422
25.4674
(17.5,12.5)
l = 4.60247P − 142.83603
0.95935
31.0346
(20,15)
l = 3.70937P − 131.56198
0.97198
35.4675
45
(5,0)
l = 87.80175P − 705.51258
0.86615
8.0353
(7.5,2.5)
l = 17.68538P − 231.72743
0.89638
13.1028
(10,5)
l = 9.49899P − 176.70845
0.96058
18.6029
(12.5,7.5)
l = 6.35449P − 150.48197
0.9683
23.6812
(15,10)
l = 5.23749P − 159.45243
0.98428
30.4444
(17.5,12.5)
l = 4.10701P − 142.96515
0.9775
34.8100
(20,15)
l = 3.35472P − 132.33265
0.96648
39.4467
90
(5,0)
l = 22.25843P − 69.97357
0.31616
3.1437
(7.5,2.5)
l = 7.89522P − 163.5916
0.95204
20.7203
(10,5)
l = 5.82953P − 153.31147
0.95795
26.2991
(12.5,7.5)
l = 4.51084P − 138.39713
0.96233
30.6810
(15,10)
l = 3.53828P − 120.16012
0.96753
33.9600
(17.5,12.5)
l = 2.86423P − 106.77622
0.96722
37.2792
(20,15)
l = 2.23186P − 84.21847
0.98777
37.7346Figure 18
The variations of crack initiation pressureP0 and the slope of P-l linear fit al versus σ3 under different notch direction conditions.


## 6. Conclusions

To better understand the mechanism of SCA-induced crack propagation in deep underground, crack initiation and propagation under different stress conditions were modelled using a modified cohesive element method. A new traction-separation law for describing rock compressive shear strength was proposed. Crack length and crack direction in bidirectional isobaric and unequal stress fields were analyzed in detail. According to the numerical results, the following conclusions are drawn:(1)
The crack morphology modelled using the modified cohesive element method is similar to that of the experimental results, indicating that the method proposed in this study can model SCA-induced crack initiation and propagation.(2)
The crack initiation pressureP0 and the incremental ratio of crack length to SCA expansion pressure al can be used as useful indicators to evaluate the difficulty to break rock in deep underground. High P0 and low al make rock breaking more difficult.(3)
The maximum expansion pressure of SCAs depends on the elastic modulus and geostress field of the rock and therefore should be measured under conditions similar to deep underground before designing SCA borehole spacing, which should be less than 220 mm at 600 m deep with a maximum SCA expansion pressure of 100 MPa. To break rock more effectively, borehole spacing should be reduced under higher stress condition.(4)
σ3 is the main factor that affects the crack initiation pressure P0. Both the stress ratio σ3/σ1 and notch direction control the direction of crack propagation. When the notch direction is oblique with the σ1 direction at 45°, the angle between the crack and the σ3 direction θc decreases linearly as σ3/σ1 increases, and their linear relationship is(28)θc=83.61466−38.69622σ3σ1.



---

*Source: 7936043-2018-11-01.xml*