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

School of Energy Science and Engineering, Henan Polytechnic University, Jiaozuo 454000, Henan, China Department of Mining and Materials Engineering, McGill University, Montreal, QC, Canada H3A 0E8 Henan Polytechnic University, Jiaozuo 454000, Henan, China State and Local Joint Engineering Laboratory for Gas Drainage and Ground Control of Deep Mines, Henan Polytechnic University, Jiaozuo 454000, Henan, China


Introduction
e morphology, direction, and propagation length of artificially induced cracks are the important parameters for rock or concrete fracturing projects [1][2][3]. e 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].
SCAs primarily consist of calcium oxide (CaO) and some other additives such as magnesium (MgO), aluminum oxide (Al 2 O 3 ), and silicon dioxide (SiO 2 ). 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. e 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]. e 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. us, 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 lled 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 sti ness and con ning pressure of the surrounding rock. Tang et al. [24] investigated its feasibility to improve coal permeability. Temperature e ect on SCA expansion pressure was also studied numerous times [15,[25][26][27]. For deep rock projects, the initiation, direction, and length of SCA-induced cracks are in uenced by stress eld. However, very few studies focused on the e ects of geostress eld 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 nite 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]. e 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 di cult 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 prespeci ed cracking path. erefore, VCCT is not suitable for modelling rock fracturing. e 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.
is study aims to investigate the characteristics of SCA-induced crack propagation in rock under di erent stress conditions. As shown in Figure 2, the notch direction is a signi cant parameter for controlling the location of crack initiation and direction of crack propagation, so the notch direction was also considered in this study. is paper developed a modi ed cohesive element method that uses rock compressive shear strength to model SCA-induced crack initiation and propagation. is study is of signicance to deepen the understandings of the characteristics and mechanism of SCA-induced crack propagation in deep underground.

Methodology for Modelling Crack
Propagation in Rock Based on the cohesive element method, crack initiation and propagation can be simulated using the nite element method. As shown in Figure 3(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 sti ness, is not equal to zero. e 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(b) illustrates the process of generating zerothickness cohesive elements in a solid element mesh. e algorithm used was written in the Python programming language and implemented as an ABAQUS plugin. By executing the plugin, cohesive elements are generated (1) (1) at every solid element interface; therefore, the nal 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. e mesh should be ne and irregular to avoid mesh structure e ect on crack propagation paths and to improve crack propagation randomness. In this study, the irregular mesh was generated using the advancing front technique. e nal mesh model is shown in Section 4.3.

Initial
where K s is the initial sti ness of the solid element. When embedding one cohesive element between two solid elements, the new total sti ness K o ′ can be calculated by where K c is the initial sti ness of the cohesive element. According to Equations (1) and (2), the total sti ness of the nite element model is a ected by embedding cohesive elements. To avoid this e ect, K s must be much greater than K c in Equation (2). e relationship between element sti ness K and elastic modulus E is where h is the constitutive element thickness. For a solid element, h is the same as its geometric thickness. For a zerothickness 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 in uence the global mesh sti ness. In this study, the elastic modulus of cohesive element and solid element are the same, and the initial constitutive thickness of cohesive element h oc is one thousandths of the constitutive thickness of solid element h s .

Damage Initiation and Evolution of
Cohesive Element e 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 agged as negative values. As shown in Figure 5, once the traction exceeds the strength, the cohesive element begins to damage, and its sti ness decreases as separation displacement increases. When the cohesive element is completely damaged, its sti ness reduces to zero, indicating forward crack propagation. e traction-separation law can be divided into three stages: elastic stage, damage initiation stage, and damage evolution stage.

Linear Elastic Traction-Separation Equation.
It is assumed that traction-separation behavior in the elastic stage Advances in Civil Engineering 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 where ε n and ε s are the nominal tensile and shear strain components, respectively, and δ n and δ s are the corresponding separation displacement components, respectively. e nominal traction stress vector t consists of two components in 2D problems: t n and t s , 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. e elastic relationship between nominal stresses and nominal strains is represented as where E n and E s are the tensile and shear elastic moduli, respectively.

Damage Initiation.
In deep underground surrounding a SCA borehole, stress eld is a ected by the SCA's expansion pressure and crustal stress eld. Rock shear strength usually increases with con ning pressure; therefore, compressive shear strength characteristics are crucial for crack propagation modelling.

Compressive Shear Damage Criterion (t n < 0).
When t n < 0, the damage initiation of cohesive elements is assumed to obey the Mohr-Coulomb criterion. e damage criterion function is given by where S is the pure shear strength (i.e., the cohesion) of rock and φ is the friction angle. In the t n -t s plane, Equation (6) is plotted as two radials, as illustrated in Figure 6.

Tensile-Shear Damage Criterion (t n ≥ 0).
e cohesive element begins to damage when the nominal stress ratios in the quadratic interaction function below reaches 1. is criterion is represented as where T is the tensile strength of cohesive elements. For rock materials, T is usually smaller than S. erefore, Equation (7) represents a semiellipse in the t n -t s plane, as shown in Figure 6. Equation (7) describes the mixed mode of tensile and shear stresses that causes damage initiation.

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 e ects of all the active damage mechanisms. D is zero initially and then gradually evolves from 0 to 1 upon further loading. e stress components of the traction-separation law are a ected by the damage according to the following equation: 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 sti ness is degraded when t n ≥ 0 but unchanged when t n < 0. e mixed mode happens when cohesive elements undergo normal and shear tractions simultaneously. To describe damage evolution of this mode, an e ective displacement δ m [35] is de ned as where the operator is de ned as erefore, δ m is equal to |δ s | when δ n < 0 (i.e., t n < 0). According to Equations (4) and (5), tractions at the damage initiation point are calculated by where t o n and t o s are the normal and shear stress components at the damage initiation point, respectively, and δ o n and δ o s are normal and shear displacement components at the damage initiation point, respectively (Figure 7). e mode mixity ratio β is de ned as

Advances in Civil Engineering
According to Equation (9), the e ective displacement at the damage initiation point, δ o m , is given by According to Equations (12) and (13), δ o n and δ o s can be represented by β and δ o m as where δ o m can be obtained by substituting Equations (11) and (14) into (6) and (7), and solving for δ o m gives Equation (15) indicates that δ o m 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 where G I and G II are tension and shear energy components, respectively, and G IC and G IIC are I-mode and II-mode fracture energies, respectively. When G I and G II satisfy Equation (16), the damage variable D reaches 1. G I and G II are calculated as where δ f n and δ f s are the tensile and shear displacement components at the complete damage point.
According to Equations (9) and (12), δ f n and δ f s can be represented by β and δ f m : where δ f m is the e ective displacement at the complete damage point. When α 1, δ f m is obtained by substituting Equations (18) and (19) into (16), given by Under compressive shear loading, according to Equation (12), β tends to in nity. us, δ f m can be solved by letting β → ∞ in Equation (20), which is given by For the linear softening evolution, D is expressed as where δ max m is the maximum e ective 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).
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. e damage initiation and evolution laws were written in the Fortran programming language as a UMAT subroutine.

Mechanical Properties of Rock.
e mechanical properties used in this study were obtained from laboratory tests. e 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. e 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 fourpoint shear tests, respectively. Table 1 lists the mechanical properties. Note that the elastic modulus E o , cohesion c, friction angle φ o , and tensile strength s t of the solid elements are equal to the elastic modulus E n , pure shear strength S, friction angle φ, and tensile strength T of the cohesive elements.

SCA Expansion Property.
e variation of SCA expansion pressure over time can be measured in a con ned state. As shown in Figure 10, SCAs are put into a thickwalled steel cylinder tube, and then the SCA expansion property will be re ected 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 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 P in and P ex are the pressures of the inner surface and outside surface.
According to Equation (23), the stress components on the outside surface of the steel cylinder tube are calculated by According to the generalized Hooke law, the tangential strain component is where E p is the elastic modulus of the surrounding materials. By substituting Equation (24) into Equation (25), P in can be obtained as Equation (26) shows that P in is in uenced by both E p and P ex . Borehole expansion pressure changes with P ex , so the SCA expansion property measured under P ex 0 is not suitable for numerical simulation when ground stress eld 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. erefore, the loading rate of expansion pressure has very few e ects 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.

Mesh Model.
In fact, SCA boreholes must have an opening that allows it to expand freely along the axial direction. e free expansion behavior along the axial Advances in Civil Engineering 7 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 Figure 11, the mesh is a 500 mm × 500 mm square. e SCA borehole diameter is 30 mm, and the distance between the two notch tips is 56 mm. e 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. e bottom of the model was xed in the vertical direction. e maximum principal stress σ 1 and the minimum principal stress σ 3 were applied in both the vertical and horizontal directions, respectively. e numerical simulation consists of two steps: rst, σ 1 and σ 3 were applied to the corresponding model boundaries, and the geostatic stress was equilibrated. en, SCA expansion pressure was gradually applied to the borehole wall from 0 MPa to 150 MPa.

Test Schemes.
A series of numerical tests were carried out to study the characteristics of SCA-induced crack initiation and propagation under di erent stress conditions and notch directions. Stress conditions were categorized into σ 1 σ 3 0 MPa (free boundary conditions where the rock sample is not xed and not loaded with any external forces), σ 1 σ 3 ≠ 0 MPa (bidirectional isobaric stress eld), and σ 1 −σ 3 5 MPa (bidirectional unequal stress eld), as shown in Table 2, where θ is the angle between notch direction and σ 1 direction (Figure 2).

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 rst. e numerical results were compared with   Advances in Civil Engineering    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]. e 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]. is 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. e 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. e crack angle difference in Figures 12(a) and 12(b) is due to different mineral grain and natural microfracture distribution. e simulations in this study did not focus on the mesoscopic heterogeneity of rock materials. e 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.

Crack Characteristics in Bidirectional Isobaric Stress
Fields. e numerical results of crack propagation under σ 1 � σ 3 � 10 MPa conditions are illustrated in Figure 13.   Variations of crack initiation pressure P 0 and the slope of P-l linear ts a l as con ning stress increases.  Advances in Civil Engineering ere are three main cracks under the free boundary conditions, but only two main cracks in the bidirectional isobaric stress eld 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.
To apply SCAs in deep underground, the spacing of SCA boreholes, which depends on the crack length l, is crucial, so further analysis of the relationship between P and l under di erent stress conditions should be undertaken. P-l linear t in di erent bidirectional isobaric stress elds 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 eld, SCA-induced crack length increases with its expansion pressure. e SCA crack initiation pressure P 0 is a proper parameter for evaluating SCAs' rock breaking ability (the higher the P 0 , the more di cult it is to break rock). Under the free boundary conditions, P 0 is only approximately 4-6 MPa [11]. Under σ 1 σ 3 5 MPa, P 0 is 16.6236 MPa. However, when σ 1 σ 3 20 MPa, P 0 increases to 37.8456 MPa, which is approximately 6.3-9.5 times that under the free boundary conditions. e slope of the P-l linear t (unit: mm/MPa) a l is the incremental ratio of crack length to SCA expansion pressure. As shown in Figure 14(b), a l decreases as con ning stress increases. e lower the a l , the more di cult it is to break  e reduction rate of a l gradually slows down as con ning 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 eld does not signi cantly a ect the crack direction, but can signi cantly a ect crack initiation pressure P 0 and crack length l.
According to Equation (26), the maximum SCA expansion pressure depends on the rock elastic modulus and geostress eld, 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.

Crack Characteristics in Bidirectional Unequal Stress
Fields.
e 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. eir linear relationship is 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.  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 P 0 , while both the stress ratio σ 3 /σ 1 and notch direction control crack direction. e corresponding P-l linear fit equations are listed in Table 4. Figure 18 shows the variations of P 0 and a l versus σ 3 when σ 1 − σ 3 is kept constant at 5 MPa. According to Figure 18, P 0 increases with σ 3 under all notch direction conditions; P 0 increases with θ under all stress conditions. e 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 P 0 -σ 3 curves are approximately linear, but when θ � 90°, the trend flatlines as σ 3 > 12.5 MPa. Under all notch conditions, a l is inversely proportional to σ 3 . e results suggest that crack generation becomes more difficult as σ 3 and θ increase.

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 tractionseparation law for describing rock compressive shear strength was proposed. Crack length and crack direction in bidirectional isobaric and unequal stress fields were analyzed