1 Introduction
The tensile strength of shale is a critical parameter in shale gas exploitation because tensile failure usually occurs in the process of hydraulic fracturing. The indirect tensile strength can be determined by the Brazilian splitting method as described by ISRM (
ISRM, 1978) and is the basis of extensive studies conducted on the tensile strength and failure mode of isotropic media (
Mellor and Hawkes, 1971;
Saksala et al., 2013;
Yuan and Shen, 2017). Furthermore, the Brazilian splitting method also enables the effects of sample size, boundary conditions, loading configurations and stress distribution on the shale failure behavior to be clarified (
Rocco et al., 1999; Li and Wong, 2013).
However, the results of a large number of studies regarding rock as a kind of homogeneous medium are difficult to directly apply to shale-related engineering because real shale is actually a heterogeneous material containing cracks (
Li et al., 2014;
Liu et al., 2021). To discuss shale anisotropy,
Chen et al. (1998) and
Claesson and Bohloli (2002) established a relationship between the tensile strength of layered rocks and the five isotropic elastic parameters in the transverse direction.
Wang et al. (2016) unfolded the influence of the inclination angle (
θ) on the acoustic emission characteristics and crack initiation.
Tan et al. (2015) studied the tensile strength and failure mode of anisotropic shale by a combination of numerical and experimental methods, suggesting that the anisotropy plays a significant role in the stress distribution, crack initiation and bearing capacity of shale.
Wang et al. (2017) further discussed the tensile strength anisotropy of Longmaxi shale at a strain rate of 10
−5 s
−1 to ~10
−2 s
−1, where more complex fractures tend to be formed under a higher loading rate.
Meanwhile, depending on the Brazilian splitting technique applied, disc specimens containing one or multiple preexisting cracks have been adopted by many researchers who are interested in fracture mechanics and crack initiation, propagation and coalescence (
Guo et al., 1993;
Fowell and Xu, 1994;
Fowell et al., 2006;
Haeri et al., 2014a). It is also noted that
θ and
L play a crucial role in determining the stress intensity factors (SIFs), among the geometrical parameters of a central straight-through crack Brazilian disc (CSCBD) specimen, where the stress intensity factors (SIFs) decrease with decreasing
θ when the ratio of the crack half-length to the radius (
a/
R) is less than 0.3 (
Ghazvinian et al., 2013).
Dong (2008) recommended that the range of the
a/
R ratio is 0.4–0.6 to improve the precision of the mode I and II stress intensity factors (
KI,
KII).
Cai (2013) revealed the initiation and propagation process of the primary wing crack and other secondary cracks based on the FEM/DEM method, indicating that the wing crack initiation positions deviate from the preexisting crack tips with increasing frictional resistance. Moreover, the specimens embedded within double or multiple cracks in the disk center were studied by
Haeri et al. (2014b,
2015), which showed that wing cracks were initiated at the first stage of loading and then propagated toward the compressive loading direction. In addition,
Sarfarazi et al. (2017) showed that crack spacing (
H) and number significantly affect the failure mode and crack initiation and failure stress, according to the PFC
2D-based investigation on the stress distribution and crack initiation, propagation and coalescence in disc specimens with multiple center cracks.
In recent years, the combined effects of preexisting cracks and rock heterogeneity have attracted increasing attention. For example,
Suo et al. (2018) exhibited the fracture toughness and fracture propagation of anisotropic shale rocks based on experimental and numerical attempts. Another example is that
Cai and Kaiser (2004) and
Zhao et al. (2021) simulated the tensile behavior of rock samples with anisotropic properties and preexisting cracks.
Yang and Huang (2014) discussed the failure strength and failure mode of heterogeneous rock with one preexisting crack and accordingly explained the shear failure and tensile failure mechanisms of specimens with different notch angles. Basically, previous works have concluded that rock anisotropy and preexisting cracks are key parameters in the determination of tensile crack initiation and propagation, which is helpful for understanding the tensile mechanical behavior of jointed rocks. Unfortunately, regarding the existing numerical simulations related to shale Brazilian disc specimens, discs of anisotropic shale with multiple cracks have rarely been studied.
In this paper, PFC2D was used to determine the tensile strength and the corresponding mechanism of anisotropic shale specimens embedded with two parallel cracks under Brazilian splitting conditions. First, a set of micro- and macro-parameters was obtained from the intact disc experimental results based on the trial-and-error method for input into PFC2D. Then, a series of numerical operations were conducted to clarify in detail how α, θ, β, L, and H influence the tensile strength and failure mode, and the involved mechanisms were also determined. This numerical work is expected to gain a better understanding of the failure behavior of real shale, which is heterogeneous and usually contains complex natural fractures.
2 Numerical model development
2.1 Intact shale model and its validation
As shown in Fig. 1, an intact disc model (diameter of 50 mm) with different inclination angles
θ (
θ of 0°, 15°, 30°, 45°, 60°, 75°, and 90°) is established in PFC
2D for microparameter calibration and model validation by comparison with experimental test results. According to previous studies, reasonable results can be obtained when the model size is 100 times larger than the average particle radius and the ratio between the maximum particle radius and the minimum particle radius is close to 1.66 (
Potyondy and Cundall , 2004;
Yang et al., 2006).
Zhang et al. (2017) showed that the laminar thickness (
l) of Longmaxi shale, which has the same geological background as the shale studied in this paper, ranges from 0.1 mm to 1.0 mm. Therefore, in this work, the
l value was set to 1.0 mm in the numerical model, and the particle sizes ranged from 0.15 mm to 0.25 mm with a uniform distribution, ensuring the efficiency and accuracy of this numerical simulation. Herein, the smooth joint model (SJM) and the flat joint model (FJM) are used to simulate the shale layers and matrix planes, respectively. More details about SJM and FJM are exhibited by
Wang et al. (2014) and
Wu et al. (2018).
Because it is difficult to directly obtain the microparameters by experimental methods, numerical calculations are usually employed. For example,
Yoon (2007) used regression analysis to determine microparameters. However, the microparameters obtained from mathematical methods are only applicable for situations with simple models with few parameters. Hence, the trial-and-error method is generally used for parameter capture in PFC models (
Park and Min, 2015;
Park et al., 2018).
Xia and Zeng (2018) suggested a macro-/microscopic calibration method for transversely isotropic rocks in which the microparameters of particles and parallel bonds should be calibrated first by simulating the macroscopic mechanical behavior of
θ = 90°. Then, the microparameters of the SJM were calibrated at
θ = 0° and 90°; this method avoids the defects that arise in the
Park et al. (2018) method, which would cause a large error after adding the bedding planes. Herein, Xia’s method was adopted to calibrate the shale microparameters, and the key microparameters are shown in Table 1
To stimulate quasistatic loading, it is necessary to keep the load plate moving extremely slowly during loading. Therefore, a loading rate of 0.2 m/s was adopted in the numerical model, and the time step was 2.36 × 10−8 m/step. That is, the load plate moving 1 mm requires approximately 211864 steps.
To verify the reliability of the numerical model, samples (with a diameter of 50 mm and height of 25 mm) are obtained from a Longmaxi shale outcrop located within the Pengshui zone, Chongqing, China. Three samples are investigated for each θ to reduce experimental dispersion. The indirect tensile strength can be calculated by Eq. (1).
where
F is the maximum splitting load,
D is the diameter of the disc specimen, and
t is the thickness of the sample. Although Eq. (1) is established for linear elastic, homogeneous and isotropic materials and does not truly represent the tensile strength of transversely isotropic materials, it is very useful to compare the effect of the inclination angle (
Tavallali and Vervoot, 2010).
Figure 2 shows the comparisons of the tensile strength results obtained by numerical and experimental outputs under Brazilian testing for intact shale specimens, where the failure strength is anisotropic and decreases with increasing θ in the case of θ≥15°, and the strengths corresponding to θ = 0° and θ = 15° are almost equal, that is, the simulated tensile strength agrees well with the experimental results under Brazilian test. Meanwhile, Figure 3 shows that as θ ranges from 0° to 90°, the fracture mechanism and failure modes obtained by the numerical method are in good agreement with those from the experimental operations. In summary, the comparisons between numerical and experimental outcomes exhibit the reliability of the numerical model proposed in this work. Thus, the parameters used in the numerical model are applicable to the analysis of shale Brazilian disc specimens.
2.2 Model development for Brazilian discs containing two parallel cracks
To discuss the influence of preexisting cracks on the mechanical behavior of transversely rocks, such as shale, numerical specimens containing two parallel notches are created by weaking the contact parameters in the middle region of the disc specimen (Fig. 4), and the key microparameters of the notches are shown in Table 2. The other parameters of the numerical model are the same as those of the intact shale model. Herein, the ligament angle α (the angle between the ligament direction and the loading direction), notch angle β, notch length L and notch spacing H were considered in a comprehensive manner. Note that the β value varies from 0° to 90°, with an interval of 15°.
3 Results and discussion
3.1 Effect of ligament angle (α) on tensile strength of Brazilian disc specimens with two parallel cracks
As previously mentioned, Eq. (1) is used for comparison between only intact and preexisting crack samples. To discuss the effect of the ligament direction, which is perpendicular to the loading direction (LVL) or parallel to the loading direction (LPL), and the β value on the tensile strength, the crack length and crack spacing remain constant, where L = 10 mm and H = 25 mm.
In this work, the relative strength metric is defined as the ratio of the tensile strength of specimens containing preexisting cracks to that of intact specimens. That is, the relative strength is a parameter that is introduced here to judge how the preexisting cracks influence the failure behavior of shale specimens so that the impact of the notch angle on the tensile strength can be discussed. Herein, the label TiBj (i, j = 0, 15, 30, 45, 60, 75, 90) means that the modeling specimen has an inclination angle of i degrees and contains two parallel cracks with a notch angle of j degrees.
For the LPL case, as shown in Fig. 5, all the specimens containing cracks (except T0B0) have a lower indirect tensile strength than the intact specimen. When θ remains constant, the relative strengths of the specimens roughly form a U-shaped curve with increasing β. For all the specimens, the tensile strength decreases gradually as β increases from 0° to 45°, where the decrease is low in the β = 0°–15° phase and tends to be high in the β = 15°–45° phase, as revealed by the curve slope in Fig. 5. For θ = 0°–90°, changing at an interval of 15°, the minimum relative strengths are observed in the cases of β = 60°, 45°, 45°, 45°, 75°, 60°, and 45°. The minimum relative strengths are 38.2%, 41.1%, 47.1%, 42.0%, 50.5%, 55.6%, and 47.0%. Then, the relative strength increases slowly as β increases. Among all the specimens, in the case of θ = 75°, the preexisting cracks have the least influence on the tensile strength of the samples.
According to the numerical outcomes under LVL conditions (Fig. 6), the β value has an insignificant effect on the tensile strength of the specimen when θ is 75° or 90°, where the greatest reduction rate is just 5.9% (T75B60) compared to the intact specimen. However, for θ≤60°, the specimen tensile strength exhibits a U-shaped distribution with increasing β value, among which the β value induces the greatest effect on the tensile strength of the specimen when θ is 45°. As a result, the relative strength of a specimen is the smallest at θ = 45°, and its maximum decrease is 38.2% at T45B45 compared to the intact specimen. In the case of T0B0, T0B15, T15B0, T30B15, T45B0, T75B0, T75B45, T75B75, T75B90, and T90B90, the tensile strength is increased by 0.3%–5.2% compared to the intact specimen.
In summary, the ligament angle can play an important role in determining the tensile strength of Brazilian disc specimens with two parallel cracks. Therein, preexisting cracks have a significant and insignificant influence on the tensile strength of the specimen at α = 0° and at α = 90°, respectively.
3.2 Effect of crack length (L) on tensile strength
To discuss the effect of crack length on indirect tensile strength, the constant parameters of α = 0°, β = 45° and H = 25 mm are set. For most of the crack lengths (L = 5, 12, 15, and 20 mm), the tensile strength of a specimen peaks at an inclination angle of 15° (Fig. 7), while the greatest value of tensile strength is obtained at an inclination angle of 0° when L is 10 mm. Basically, the tensile strength decreases with increasing θ, where the ratio of the maximum tensile strength to the minimum tensile strength is 1.94, 1.82, 1.74, 1.54, and 1.45 when the crack length L is 0, 5, 10, 12, and 15 mm, respectively. Hence, the degree of specimen anisotropy is reduced because of the presence of the prefabricated cracks and even decreases more at a greater L value. However, when the crack length is 20 mm, this ratio becomes 1.66, which is higher than that in the situation in which the crack length is 12 or 15 mm. This phenomenon suggests that the anisotropy caused by L and θ becomes more obvious when the crack length reaches a certain value.
According to Fig. 8, the tensile strength decreases as the crack length increases, where the rate of decrease is high at first (L<12 mm) and then becomes lower (L>12 mm). The percentage of the tensile strength of the specimen at L = 12 mm to that at L = 0 mm (intact specimen) varied from 32.6% to 49.7%, suggesting that the tensile strength was reduced by 51.3%–67.4% when L increased from 0 to 12 mm. Similarly, it could be calculated that the tensile strength was reduced by only 1.12%–9.81% from L = 12 mm to L = 20 mm. Accordingly, it can be reasonably speculated that the influence of a lower ratio of crack length to sample diameter (L/D) to the tensile strength of shale specimens is greater than that of a higher ratio.
3.3 Influence of the crack spacing (H) on tensile strength
To clarify how the H value affects the indirect tensile strength, constant parameters of α = 0°, β = 45°, and L = 10 mm are adopted. Generally, as shown in Fig. 9, the parameter H has little effect on tensile strength. The ratio between the maximum and minimum tensile strength of the specimens is 1.49, 1.49, 1.55, 1.60, and 1.74 when the crack spacing H is 5, 10, 15, 20, and 25 mm, respectively, while the corresponding ratio for the intact specimen is 1.94. This phenomenon indicates that the prefabricated cracks can reduce the anisotropy degree of modeling shale specimens, where the anisotropy becomes more unobvious at a smaller crack spacing.
4 Failure mode recognition
As previously mentioned, the comprehensive action of α, θ, and β has a great impact on the tensile strength of shale under Brazilian splitting. To further discuss the effect of the comprehensive action of α, θ, and β on the shale mechanical behavior, the failure modes and micromechanism were analyzed with different α, θ, and β values.
4.1 Failure mode of LPL
Figure 10 shows the numerical load-displacement curve for specimen T60B45 under LPL conditions. As shown in Fig. 10, point a represents the early loading stage in which the microcracks first initiate from the tip of the notch. Point b represents the stage when the peak force was reached, and point c represents the stage when the loading was stopped. To investigate the deformation progress of shale containing two parallel cracks, three interesting zones are distinguished by kl (k, l = a, b, c), where the first letter k represents one of the three loading stages in Fig. 10 and the second letter l represents one of the areas in the samples (the zone between the upper loading point and upper notch, the zone between the two notches, and the zone between the lower loading point and lower notch, respectively). At point a of Fig. 10, before the peak load, the maximum contact force between particles occurs near the loading point, where the particles move significantly faster than those in the specimen center (aa and ab in Fig. 11(b)). However, the particle displacement near the loading point is mainly toward the sample center. In addition, the particle displacement vectors at the center of the sample begin to diverge from each other (ab in Fig. 11(b)). Meanwhile, the particles on the left side of the notch move toward the left, while those on the right side of the notch move toward the right. As a result, some microcracks first initiate from the tips of the notch (ac in Fig. 11(b)).
When the load force increases to the peak force (point b in Fig. 10), as shown in Fig. 12, the particle displacement at the notch tips continues to increase. Several macrocracks that form by microcrack coalescence are initiated from the tips of the notch and loading point (Fig. 12). At this moment, the dispersion phenomenon of particles in the specimen center is more obvious. The macrocracks near the notch tip and loading point gradually approach each other and finally merge (ba and bc in Fig. 12(b)). Similarly, the cracks near the tips of the two notches gradually spread toward each other (bb in Fig. 12(b)). Meanwhile, the contact force of the particles at the notch tip still remains high, indicating that the macrocracks will continue to propagate forward.
When the load decreases to point c of Fig. 10, the macrocracks penetrate the entire specimen between the two loading points, as shown in Fig. 13. The wing cracks between the notches and those between the notch and the loading point are mainly characterized as tensile failures because the displacement vectors of the particles diverge from each other. However, for some secondary cracks at the nearby loading points, two mechanisms are recognized. One is the different displacements of the particles passing through the bedding layer in the same direction, leading to tensile failure along the bedding planes (black segments of ca in Fig. 13(b)). Another is that the displacements of the particles on the left and right sides of the bedding share the same direction as those in the bedding layer; however, the displacement magnitudes are different, which leads to shear failure (pink segments of Fig. 13(b)).
Basically, there are four typical failure modes for specimens under LPL conditions induced by the variables θ and β, as shown in Fig. 14. In the case of β = 0°, the macrocracks propagate in a straight path (approximately), parallel to the loading direction (mode I), as shown in Appendix I. In the cases of β = 15°, 30°, and 45°, some macrocracks initiate from the tips of the notch away from the loading point (f1, f2) and continue growing parallel to the loading direction, ultimately coalescing in the bridge area. Here, these cracks form a closed quadrilateral or triangular area at the bridge area (i.e., the area between the two preexisting cracks), revealed by mode III in Fig. 14. Note that when β = 15° and θ = 45° and 60°, the crack propagates rapidly from f1 to n2, while the crack propagation from f2 to n1 is suppressed due to the bedding effect, by which failure mode II occurs. In the case of β≥60°, the inclination angle has almost no effect on the failure mode, while the notch angle plays an important role in crack propagation. Therefore, notch tips f1 and f2 are connected, and tips n1 and n2 are connected to the lower and upper loading points, respectively, where failure mode IV occurs.
4.2 Failure mode under LVL conditions
Similar to the LPL case, three interesting zones are distinguished by kl (k, l = A, B, C), where the first letter k represents one of the three loading stages in Fig. 15 and the second letter l represents one of the areas in the samples (the zone near the upper loading, the zone between the two notches, and the zone near lower loading, respectively).
Figure 16 shows the contact force and displacement vector distribution corresponding to the LVL case at point A in Fig. 15 (before the peak load). Different from the LPL case, the bedding plane shear cracks are initiated from loading points because 1) the particle displacement along the bedding direction is large near the loading points (AA and AC in Fig. 16(b)) and 2) the displacement direction of the particles around the notches is approximately the same; thus, their relative displacement is small (AB in Fig. 16(b)). In addition, the maximum contact force between particles occurs near the loading point. The displacement of particles in the middle of the disc specimen is directed diagonally upward and downward along the bedding planes, suggesting that the bedding planes significantly affects particle displacement.
As the load force increases to the peak force (point B in Fig. 15), as shown in Fig. 17, both the particle displacement and contact force continue to increase. Some particles near the bedding planes move along the bedding direction, while others move in the direction normal to bedding. Comparatively, the displacement difference among the particles moving along the bedding direction is greater, which mainly induces shear cracks (pink segments), accompanied by a small number of tensile cracks (black segments) along the bedding planes (BA and BC in Fig. 17(b)). Meanwhile, there are some matrix tensile cracks (green segments) and bedding plane shear cracks around the notch tips (BB in Fig. 17(b)). Although the particles in the middle of the disc specimen always displace in opposing directions (BB in Fig. 17(b)), the displacement of the particles at the center of the specimen is much smaller than the displacement at the loading points. Hence, it is difficult to first initiate cracks at the center of the specimen.
When the load decreases to point C (Fig. 15), as shown in Fig. 18, the macrocracks distribute throughout the specimen between the two loading points, mainly caused by bedding plane shear and tensile failure. Meanwhile, the cracks at the tips of the notch, resulting from mainly matrix tensile failure and bedding plane shear failure, continue to propagate toward the loading points.
There are three representative failure modes for specimens under LVL conditions according to the variables θ and β (Fig. 19). For failure mode A, as shown in Appendix II, a straight or arch line connects the loading points, leading to specimen splitting, which mainly occurs under the condition that β is small or θ is large. When β is 45°–75° and θ is smaller than 45°, the most complex failure mode arises (mode B), which not only includes the cracks connecting the loading points but also includes the cracks propagating from the notch tips to the loading points. Note that in the case of β = 75° and θ = 0°/15°, the cracks at the notch tips are suppressed by the layering and thus less likely to propagate forward; therefore, failure mode A occurs. In the case of 90°≥β≥30° and θ = 30°/45°/60°, the wing cracks propagate in a curved path and continue growing to the cracks initiated at the loading points (failure mode C).
5 Conclusions
Aiming at enhancing the knowledge of shale failure behavior, the indirect tensile mechanical behavior of anisotropic disc specimens, with two central parallel cracks, was analyzed based on modeling Brazilian splitting tests with PFC2D. The effects of α, θ, β, H, and L on the tensile strength, failure mode and mechanism were clarified, yielding the following conclusions:
1) Preexisting cracks have a significant and insignificant influence on the tensile strength of shale at α = 0° and 90°, respectively. In the case of α = 0°, the maximum tensile strength of shale with different θ and β values is reduced by up to 61.8% (corresponding to T0B60) compared to intact specimens. However, in the case of α = 90°, except for θ = 45°, the tensile strength is reduced by 6%–38.2%; under the other θ and β conditions, the tensile strength is only reduced by 0.6%–16.9%. For some specimens, the tensile strength is even increased by 0.3%–5.2% compared to intact specimens.
2) The anisotropy degree of shale specimens will be reduced with the inclusion of preexisting cracks, and in the case of an L/D ratio of 0.3, the anisotropy is characterized as the weakest. Meanwhile, a smaller crack spacing usually stimulates a greater reduction in the anisotropy degree.
3) Under LPL conditions, there are four typical failure modes at different notch angles and inclination angles. As the notch angle increases, the influence of the inclination angle on the failure mode gradually decreases. In addition, at β≥60°, the failure mode is mainly affected by the notch angle.
4) Under LVL conditions, three typical failure modes arise. In the case of 45°≤β≤75° and θ≤30°, the most complex failure mode will occur.