Abstract
The cracking behavior of precracked rocks under mechanical-hydraulic loading is of great significance in underground openings or petroleum engineering. In this study, an advanced in-house finite element code PANDAS proved to be effective in simulating coupled fracturing processes under complex geological conditions was used to simulate the cracking propagation of the precracked rocks under mechanical loading and mechanical-hydraulic loading with different strength parameters. The simulation results demonstrated that (1) the cracks initiate by the induced stresses, and multiple types of tensile cracks originate from the preexisting flaws; (2) crack propagation patterns under mechanical-hydraulic loading were studied with different strength parameters, and the multiple patterns of pure tensile, main tensile, tensile-shear, main shear, and pure shear were observed; and (3) the timing of hydraulic loading has a significant impact on the fracturing process: when hydraulic loading was carried out in the phase of main crack propagation, the tensile fracture was promoted and the shear fracture was inhibited; when hydraulic loading was carried out in the phase of shear crack propagation, the shear fracture and tensile fracture were stimulated. The numerical simulation results are in good agreement with the experimental results by previous studies. The research on the cracking behavior of precracked rocks under mechanical and hydraulic loading will expand the application prospect in the fields of coal seam gas reservoir and tunnel water inrush.
1. Introduction
In the process of underground engineering construction such as traffic tunnel excavation, water conservancy project, and coal and oil exploitation, the geological environment of high ground stress and strong seepage pressure is becoming more and more complex. Under the joint action of in situ stress and seepage water pressure, the primary cracks in engineering rock mass are easy to expand and evolve, which leads to the change of macroscopic mechanical behavior and fracture mode of rock mass [1, 2]. Preexisting cracks or fractures determine the initiation of new cracking and dominate the formation of macroscopic failure planes [3]. When a load is applied, new cracks originate from preexisting cracks and propagate by the influence of major principal stress, sometimes coalescing with other cracks. The cracking of rock masses affected by the coupling effect of mechanical-hydraulic loading is increasingly investigated in engineering.
In past decades, a large number of experimental studies were conducted to characterize the cracking behaviors and the law of fracture in rock samples containing preexisting flaws, where the term flaws refer to those preexisting artificial cracks [4–6]. The fracturing processes and crack coalescence patterns under compression were characterized and estimated on different rock materials [4]. Some studies suggested that tensile cracks are initially originating from flaws under compressive loading, and the tensile cracks are subsequently followed by shear cracks [5]. Another study investigated the failure mode and cracking process of marble specimens with different precracks by the uniaxial compression test. Lee and Jeon [6] carried out uniaxial compression tests on three different materials and obtained the initiation, propagation, and coalescence laws of presingle crack and predouble crack. Yin et al. [7] conducted experimental investigations into hydraulic properties of 3D rough-walled fractures during shearing. At the same time, the high temperature and included angle of preexisting cracks had been experimentally studied to observe macroscopic cracks emanating from the flaws [8, 9].
Moreover, several numerical studies were conducted to explain better the fracturing processes in precracked models, including the finite element methods (FEM) [10–12], boundary element methods (BEM), and displacement discontinuity methods (DDM) [13]. To model the crack behaviors reasonably, appropriate cracking criteria need to be incorporated [14]. A numerical simulation code (RFPA) can be used to analyze the crack propagation in a rock while simulating both global failure of rock and local cracking at flaw tics [15]. Also, cracking behavior in rock and rock-like materials has been modeled by using particle flow code (PFC) [16].
Numerical models are not only useful for analyzing and interpreting the crack types observed in physical experiments but also useful tools for observing the dynamic images of the stress field and seepage field evolution in the process of fracture propagation [17, 18]. However, many studies have only focused on the analysis of preset single crack growth, which cannot completely simulate the whole process of crack initiation, development, and interpenetration under different strength parameters and hydraulic coupling environment [19, 20]. Therefore, some more complex conditions of the actual work need to be further considered to simulate the cracking behavior and propagation mechanism of precracked rock, such as different strength parameters and hydraulic coupling.
To address this knowledge gap, this study was based on the porous mechanical-hydraulic coupled FEM and the continuum damage theory, and an in-house finite element code PANDAS was used to simulate the cracking behavior of precracked rock [21–23]. Moreover, the cracking behaviors that originate and propagate from the flaws with different strength parameters in rock samples were also investigated in this study. On this basis, how the crack behaviors propagate under the mechanical-hydraulic loading and the influence of hydraulic action time sequence was examined. The numerical simulation strategy considering the influence of different rock strength parameters and hydraulic time series on crack growth is innovative. This research is expected to enhance the awareness of the instability process of crack damage evolution under complex stress environments, as well as the mechanism of crack formation and propagation.
2. Methodologies
Many scholars have used the numerical simulation method (RFPA, PFC, and FLAC3D) to study the fracture mechanism of the fissured rock mass, which helps to understand the stress state and provides a powerful criterion for the initiation and propagation of cracks [10–20].
In this study, an in-house code PANDAS (Parallel Adaptive Nonlinear Deformation Analysis Software) was used to simulate the crack propagation patterns of a precracked model under mechanical-hydraulic loading [21]. PANDAS is an advanced computing program in which the in-house finite element code has implemented the advanced unstructured mesh generation of extremely heterogeneous fractured porous media by using a nonlinear coupled finite element solver [22, 23]. It is a finite element code developed and verified for simulating coupled fracturing processes in a complex geological condition, which were set in 3D subjected to the affecting hydraulic and thermal factors for simulating coupled crustal dynamics [24, 25]. The related equations used in this study are briefly listed below [5, 11, 26–29].
The stress balance equation is defined aswhere is the component of the Cauchy stress tensor and is the body force in the th direction.
Jaumann rate is introduced for the stress:where is the Jaumann rate of the Cauchy stress tensor and is the elastoplastic material matrix. The geometrical equation is defined aswhere is the partial derivative of the th velocity at the th direction and is the velocity gradient tensor. and are the symmetric and antisymmetric parts of , respectively.
The fluid and stress coupling equation is defined aswhere is the pore pressure, is the coefficient of pore pressure, and is the Kronecker constant.
The seepage equation is expressed as follows:where is the permeability, is the Biot constant, and is the volume strain.
The permeability follows the exponential law of effective stress:where is the initial permeability and and are the material constants. Once the element is failed, relevant stresses are released, and the element is granted with the properties of air material (i.e., Young’s modulus is reduced to a minimal value, and the permeability is corresponding to a fractured state).
Maximum principal stress criterion and Mohr-Coulomb theory are used as the cracking criteria of the poroelastic-brittle rock sample based on the effective stress concept. Tensile cracks occur when the maximum tensile stress exceeds the tensile strength:
Shear cracks are determined by the Mohr-Coulomb criterion:where is minimum principal stress, is maximum principal stress, is the friction angle, and is the shear strength.
3. Numerical Simulation
3.1. Model
To study the fracture behavior of precracked rock under mechanical hydraulic load, the calculation model is set up in this paper as shown in Figure 1. The precracked rock sample is designed at 50 mm (width) and 100 mm (height) in two dimensions (Figure 1(a)). Three loading stages are designed to include the development of crack initiation (stage 1), crack propagation (stage 2), and fracturing (stage 3). Mechanical loading is performed at stages 1–3 with applied axial displacement on the top boundary and fixed displacement at the bottom boundary (Figure 1(b)). Hydraulic loading is performed at stage 2 at both the top and bottom boundaries (Figure 1(c)). Besides, five cases of preexisting flaws with diverse angles (in the range 15°–75° and 15° in the interval) are designed at the center of the rock samples. The crack initiation and crack propagation patterns under uniaxial compression and mechanical-hydraulic loading are discussed in Sections 3.2 and 3.3, respectively.
(a)
(b)
(c)
3.2. Crack Behaviors of Preexisting Flaw under Uniaxial Compression
3.2.1. Stress Field around the Flaw and Crack Initiation
Five cases of rock samples were performed at stage 1 with a displacement up to 0.2 mm, wherein cracks were originating, and crack-induced stresses were formed. Figure 2 shows the stress field around a 60-degree flaw. The maximum and minimum principal stresses and maximum shear stress were recorded. The contour mapping of the minimum principal stress exhibits the compressive and tensile stress areas (Figure 2(a)). There are two concentrated areas of tensile stress (which resemble an auriform area) on each tip of the flaw that was recorded in the indicators of “MIN1” and “MIN2.” The contour mapping of the maximum principal stress also exhibits an auriform concentrated area on each tip of the flaw (Figure 2(b)), and the indicator “MAX” was used to record the value. The contour mapping of the maximum shear stress exhibits the orthogonal areas of the stress concentration around the flaw (Figure 2(c)), and the indicator “MAXSH” was used to record the value.
(a)
(b)
(c)
Figures 3 and 4 show the evolution of the stresses (MIN1, MIN2, MAX, and MAXSH) with the variation of flaw angles (15°, 30°, 45°, 60°, and 75°). Stress MIN1 decreases after increasing to the maximum value in the 45-degree flaw. Stress MIN2 increases continuously; stresses MAX (absolute value) and MAXSH smoothly decrease continuously.
Diverse tensile cracks originate from the flaws under compressive loading (Figure 5). Type I and II wing cracks are both emerging in the flaws at 15°, 30°, and 45° because of the stress concentrations of MIN1 and MIN2, respectively; moreover, type I crack is followed by type II crack in time history by the reason that stress MIN2 is generally smaller than stress MIN1. Type III crack emerges along the same direction of the flaw at 60°. Type IV crack occurs in the flaw at 75°, wherein cracks are simultaneously emerging from the directions of MIN1 and MIN2 by the reason that stresses of MIN1 and MIN2 are much closer.
(a)
(b)
(c)
(d)
3.2.2. Crack Propagation Patterns Affected by Strength Parameters
During stages 2 and 3 of the mechanical loading, crack propagations are highly related to the specific strength parameters in the rock samples. As shown in Table 1, the tensile strength gradually increases in the five schemes of strength parameter (i.e., SP1, SP2, SP3, SP4, and SP5), which is equivalent to the decrease of the ratio between shear strength to tensile strength. SP1 and SP5 are the minimum and maximum tensile strengths, respectively. From SP2 to SP4, the tensile strengths gradually increase.
Crack propagation patterns (e.g., in the flaw at 60°) affected by the five schemes of strength parameters are shown in Figure 6. A single main fracture forms in the schemes SP1 and SP2. The main (first) fracture and the second fracture approximately perpendicular to each other are forming in the schemes SP3, SP4, and SP5. The crack propagation patterns of the five schemes are pure tensile, main tensile, tensile-shear, main shear, and pure shear, respectively. From SP1 to SP5, the proportion of shear cracks gradually increases in the main fracture, which induces the decrease of the angle (concerning the horizontal direction) of the main fracture plane. The appearance of the second fracture is affected by the increase of the shear crack proportion in the main fracture plane—as the shear cracks increase in the main fracture from SP3 to SP5, the second fracture appears earlier and develops more significantly. Figure 7 shows the shear stress evolution of the initial crack in the second fracture. From SP2 to SP4, the shear stress develops more quickly, and the earliest stress drop (shear crack) appears in SP4. The three moments of the shear stress evolution (A, B, and C) are exhibited by the contour mapping of the maximum shear stress, wherein the shear stress (red area) develops and concentrates in the orthogonal directions leading to the formation of the fractures.
(a)
(b)
(c)
(d)
(e)
The relationships between the angles of the main fracture planes and the schemes of strength parameters are shown in Figure 8. From SP2 to SP4, the angles of the main fracture planes decrease. From 15° to 75° of the flaws, angles of the main fracture planes increase.
3.3. Crack Behaviors under Mechanical-Hydraulic Loading
To explore the hydraulic effect acting in propagated cracks, hydraulic pressures were imposed at both ends of the rock sample after stage 2 of the mechanical loading was performed. Two schemes of hydraulic loading (SH1 and SH2) were designed in the simulation (Figure 9). Under SH1, hydraulic pressures were applied after the main fracture had been developing for a specific period wherein only tensile cracks occurred (including the schemes SP2, SP3, and SP4). Under SH2, hydraulic pressures were applied at a period when the shear cracks were occurring after tensile cracking (including the schemes SP3 and SP4). Table 2 shows the process of mechanical-hydraulic loading in the rock samples (e.g., in the flaw at 60°). The three schemes of strength parameters SP2, SP3, and SP4, are included to study the conditions of nonhydraulic loading (NonHy), hydraulic loading scheme 1 (SH1), and hydraulic loading scheme 2 (SH2). Hydraulic loading was applied after the axial displacement of the rock sample reached stage 2 (crack propagation), and the sample was subsequently continually loading to stage 3 (fracturing). SH2 is not considered in SP2 because of the absence of shear cracks.
(a)
(b)
3.3.1. Mechanical-Hydraulic Action under SH1 Scheme
The contour mapping of the minimum principal stress (tensile) before and after hydraulic loading and the seepage field is shown in Figure 10. The related seepage parameters are introduced in Table 3. The initial permeability (isotropic) is the base value independent from the stress field. Coefficient 1 defines the adjustment of the permeability affected by the stress field (anisotropic), and coefficient 2 defines the increase of the permeability when cracking occurs. Figure 11(a) shows the evolution of the hydraulic pressures in SP2 by the three monitoring points MA, MB, and MC. The hydraulic pressures increase to the constant value of 3.5 MPa by following the sequence MA>MB>MC. Figure 11(b) shows the evolution of the incremental tensile stresses , and 0.713 MPa, 0.695 MPa, and 0.688 MPa are reached at the points MA, MC, and MB, respectively (following the sequence MA>MC>MB). The total tensile stresses () 4.01 MPa, 2.98 MPa, and 2.70 MPa are reached at the points MA, MB, and MC, respectively (following the sequence MA>MB>MC).
(a) Before
(b) After
(c) Seepage field
Figure 12(a) shows the evolutions of the incremental tensile stresses with different flaw angles in SP2. The value of increases with flaw angles from 15° to 75°, and the stress values follow the sequence MA>MC>MB. Figure 12(b) shows the evolutions of the total tensile stresses. The value of increases with flaw angles, and stress values follow the sequence MA>MB>MC.
(a)
(b)
Continued mechanical loading (the loading stage 3) is applied after hydraulic loading. Evolutions of the maximum shear stresses (refer to the initial cracks in the second fracture) are exhibited in Figure 13 under the conditions of nonhydraulic and hydraulic actions. Under the hydraulic action, the maximum shear stresses develop slowly, and the stress drops (shear cracks) occur lately. The comparison between SP4 and SP3 shows that the maximum shear stress develops more quickly, and the stress drop occurs earlier. According to the crack propagation patterns in Figure 14, the angles of the main fracture increase because of the hydraulic action.
(a)
(b)
(a) Scheme 2 of the strength parameters
(b) Scheme 3 of the strength parameters
(c) Scheme 4 of the strength parameters
3.3.2. Mechanical-Hydraulic Action under SH2 Scheme
Figure 15(a) shows the evolution of the hydraulic pressures in SP3 monitored by the three points MA, MB, and MC. Under SH2, the hydraulic pressures increase to the constant value 3.5 MPa by following the sequence MA>MB>MC. Figure 15(b) shows the evolution of the incremental tensile stresses (), and 0.736 MPa, 0.713 MPa, and 0.698 MPa are reached at the points MC, MA, and MB, respectively (following the sequence MC>MA>MB). The total tensile stresses () 4.17 MPa, 2.93 MPa, and 2.91 MPa are reached at the points MA, MB, and MC, respectively (following the sequence MA>MB>MC).
Figure 16(a) shows the evolution of incremental tensile stresses with different flaw angles in SP3. The value of increases with the flaw angles from 15° to 75°, and the stress values follow the sequence MC>MA>MB. Figure 16(b) shows the evolutions of the total tensile stresses. The value of increases with the flaw angles, and the stress values follow the sequence MA>MB>MC.
(a)
(b)
Continued mechanical loading was applied after hydraulic loading. The evolutions of the maximum shear stresses (refer to the initial cracks in the second fracture) are shown in Figure 17 under the conditions of nonhydraulic and hydraulic actions. Under the hydraulic action SH2, the maximum shear stresses develop quickly, and the stress drops occur earlier. According to the crack propagation patterns in Figure 14, the cracking area expanded under the hydraulic action.
(a)
(b)
4. Discussion
In this study, we investigated: the stress field around the flaw and the mode of crack initiation under different strength parameters, as well as the role of hydraulic loading in crack growth. This section will further discuss the influence of a crack growth model and hydraulic loading time sequence on the fracturing process.
4.1. On the Propagation Mode of Crack Growth
The results of this study indicate that the initiation of crack types is controlled by the stress and the angle of precrack. Firstly, the stress field around the flaw (Figure 2) reflects the position that the minimum principal stress (Figure 2(a)) and the maximum principal stress (Figure 2(b)) exert on the tip of the defect and the maximum shear stress (Figure 2(c)) exerts on the orthogonal region of the stress concentration around the defect. Secondly, Figures 3 and 4 show the evolution of the stresses with the variation of flaw angles. Types I and II wing cracks sequentially occur in the flaws at 15°, 30°, and 45°. Type III crack (along the same direction of the flaw) generally occurs in the flaw at 60°. Type IV crossed crack occurs in the flaw at 75°.
Furthermore, we identified that with the increase of rock sample strength, the shear force tends to increase, the secondary crack develops and concentrates in the vertical direction of the precrack, and the main fracture angle decreases. The proportion of shear fractures in the main fractures tends to increase, and the occurrence.
time of the second fracture is earlier and more obvious (Figure 6). Shear cracks were first found in SP4, through shear stress statistics (Figure 7).
Our numerical results are in agreement with previous experimental studies [30–32]. Lajtai [30] conducted uniaxial compression loading tests and found the following crack initiation sequence: (a) tensile fractures, (b) normal shear fractures, (c) additional normal shear fractures leading to the formation of shear zones, ans (d) inclined shear fractures (Figure 18). Experimental result shows that there are various initial defects, such as cracks and crack-like in materials. And the microcracks always appear in the defects or local stress concentration areas. Fracture behavior of rock samples should be based on stress and energy [32]. The cracking behavior in specimens containing single flaws under uniaxial compression was evaluated systematically by Wong and Einstein [31]. Three of the crack types are tensile, and three of the crack types are shear. The remaining one is of mixed tensile-shear nature, with shearing occurring adjacent to the flaw tips and simultaneous tensile opening occurring farther away. The experimental study showed that seven crack types with different trajectories and initiation mechanism (tensile/shear) could originate from the preexisting flaws under uniaxial compression (Figure 19).
(a) Type I tensile crack (tensile wing crack)
(b) Type II tensile crack
(c) Type III tensile crack
(d) Mixed tensile-shear crack
(e) Type I shear crack
(f) Type II shear crack
(g) Type III shear crack
From our study, we derived some knowledge on the crack propagation mode. With the increase of the strength parameters, the crack growth mode develops from tension to shear. And the crack growth modes from SP1 to SP5 are pure tension, main tension, tension shear, main shear, and pure shear, respectively.
4.2. Influence of Hydraulic Loading Time Sequence on the Fracturing Process
In this study, the crack propagation patterns were affected significantly by the hydraulic loading time sequence. Under SH1, hydraulic pressures were applied after the main fracture had been developing for a specific period wherein only tensile cracks occurred (including the schemes SP2, SP3, and SP4). The simulation shows that the tensile stress and the angle of the main fracture crack also increase when SP2 is applied with hydraulic loading (Figures 11 and 12). For SP3 and SP4, the stress drop (shear crack) under mechanical-hydraulic loading occurs later than that under mechanical loading (Figures 13(a) and 13(b)). It shows that the growth of tensile crack on the main fracture surface is enhanced and the shear crack is restrained by mechanical-hydraulic loading. The occurrence time of stress drop (shear crack) of SP4 is earlier than that of SP3, indicating that the increase of rock strength is beneficial to the development of shear crack (Figure 13).
Under SH2, hydraulic pressures were applied at a period when the shear cracks were occurring after tensile cracking (including the schemes SP3 and SP4). The simulation shows that the tensile stress increases with hydraulic loading (Figures 15 and 16). With the increase of the strength ratio of the rock sample, the crack propagation develops from tensile to shear. The occurrence time of stress drop (shear crack) under mechanical-hydraulic loading is earlier than that under mechanical loading (Figures 17(a) and 17(b)). The reason is that the maximum shear stress develops quickly and the shear crack occurs earlier; the shear crack growth is enhanced by mechanical-hydraulic loading. Compared with SP3 and SP4, it is also concluded that the increase of rock strength is beneficial to the development of shear crack (Figure 17).
Previous physical experiments by our team have also shown that mechanical-hydraulic loading can stimulate crack growth and fracture [33, 34]. Through the acoustic emission test under the mechanical-hydraulic loading, Chen [33, 34] found that with the increase of water pressure, the number of internal cracks in the rock increases at the time of final instability (Figure 20). These results may be explained by the fact that the splitting effect of water pressure is produced in the process of crack propagation, which leads to the initiation and propagation of new cracks. The propagation and penetration of primary crack and secondary crack may further accelerate the fracture of rock.
In conclusion, the main strength of this investigation is that the formation and propagation modes of the crack in precracked rock mass are simulated under different strength parameters and different mechanical-hydraulic loading time sequence. However, the following shortcomings exist in this study: (1) only a single precrack is simulated herein, and (2) the influence of high geotherm is not considered. Based on the complexity of internal crack defects and multifield coupling environments in the rock mass, the combination of multiple precracking and high geotherm can be considered to further explain cracking index and fracture mechanism.
5. Conclusions
In this paper, an in-house finite element code PANDAS was used to simulate the cracking behavior of precracked rocks under mechanical-hydraulic loading. The influence of different strength parameters and mechanical-hydraulic loading factors on the law of crack propagation is mainly considered. The major conclusions of this study are as follows:(1)The crack types are controlled by the crack-induced stresses and angle of the flaw. Types I and II wing cracks sequentially occur in the flaws at 15°, 30°, and 45°. Type III crack (along the same direction of the flaw) generally occurs in the flaw at 60°. Type IV crossed crack occurs in the flaw at 75°(2)Crack propagation patterns were strongly affected by the different schemes of strength parameters. As the tensile strength gradually increased, the crack propagation patterns corresponding to the five schemes included pure tensile, main tensile, tensile-shear, main shear, and pure tensile, respectively. The larger the strength ratio of the rock mass is, the earlier the second fracture (roughly perpendicular to the main fracture) appears and develops more obviously(3)The applied timing of hydraulic action has a significant impact on the crack propagation patterns. When hydraulic loading is carried out in the main crack propagation stage, the shear fracture appears later and the maximum shear stress (the initial fracture of the second fracture) develops slowly. It shows that the tensile crack is promoted during the rock fracture and the shear crack is restrained by hydraulic loading(4)When hydraulic loading is carried out in the phase of shear crack propagation, the maximum shear stress develops faster and the shear fracture appears earlier. The results show that both shear fracture and tensile fracture are strengthened under the condition of shear fracture. The splitting effect of water pressure is produced in the process of crack propagation, which leads to the initiation and propagation of new cracks and further accelerates the fracture of rock
This study not only helps to understand the mechanism of crack formation and propagation under complex stress environments but also further expands the research of new engineering scenarios, such as coalbed methane logging and tunnel water inrush [15, 35]. Based on this study, further experimental and theoretical research can be carried out by comprehensively considering multiple groups of precrack and high geotherm factors.
Abbreviations
| MIN: | Minimum principal stress |
| MAX: | Maximum principal stress |
| MAXSH: | Maximum shear stress |
| SP (1-5): | Five schemes of strength parameter |
| SH (1, 2): | Two schemes of hydraulic loading. |
Data Availability
All data, models, and code generated or used in the research process have appeared in the manuscript. According to the calculation formula and method in the paper, interested readers can reproduce these results and data and can make improvements to obtain different results.
Conflicts of Interest
The authors declared that they have no conflicts of interest in this work.
Acknowledgments
We would like to thank Dr. Quanshu Li and Dr. Zhiting Han of the University of Queensland for the helpful discussions. The related calculation was carried out by running the PANDAS software on the Savanna Supercomputer at UQ. The work is financially supported by the National Natural Science Foundation of China (Grant Numbers 52074251, 9205820028, and U19A20111) and the open fund (Grant Numbers SKLGP2015K004 and SKLGP2017Z001) from the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, China