Abstract
The effect of the bedding structure on the mechanical properties of layered shale was studied by means of experiment and numerical simulation. Based on continuum damage theory and discrete fracture network modeling method (D-DFN), a finite element model describing structural discontinuity and mechanical anisotropy of shale is established. In this model, the degradation process of stiffness and strength of shale after failure is described based on the stress-displacement relationship of elements. In order to distinguish the mechanical properties between the bedding and the matrix, a nonzero initial damage variable is set in bedding elements to show initial lower elastic modulus and strength of bedding elements compared with initially nondamaged matrix elements. The calibration of model parameters is discussed, and the simulation results are compared with the experimental results. The results show that the D-DFN method can effectively simulate the anisotropic characteristics of shale deformation and strength, which verifies the effectiveness of the method.
1. Introduction
As an energy source, oil and gas resources are superior to coal in many aspects such as transportation, heating value, and environmental protection. Therefore, since the 20th century, the proportion of oil and gas in the world energy structure has gradually increased. It is usually used as a strategic material to evaluate, plan, and manage and to formulate special policies and strategies for its development and utilization. At the same time, the problem of resource recovery and reducing environmental pollution during the mining process [1] has been the focus of attention by scholars, and in-depth research has been conducted on the problem of improving recovery, such as carbon dioxide storage [2, 3] and adding nanoparticles to the polymer [4, 5]. In addition, an in-depth understanding of the mechanical properties of rocks also plays an important role in improving the recovery of oil and gas resources in tight reservoirs. Most rock materials show different degrees of anisotropy [6, 7]. During long-term diagenesis, shales develop a large number of discontinuities such as bedding, joints, and cracks through deposition and compaction, leading to structural discontinuity and mechanical anisotropy [8–10]. The results of conventional triaxial compression tests show that shale strength is related not only to confining pressure but also to the angle between maximum principal stress and bedding planes [11]. A reasonable model to describe the structural discontinuity and mechanical anisotropy of shale is meaningful for the design of safety drilling and reservoir treatment in shale oil/gas development. A lot of experiments and theoretical studies have been done by former researchers [12, 13].
[14] carried out compression tests on different kinds of rocks. Tests results showed the effect of bedding plane orientation on the test values of elastic parameters and the yield strengths. [15, 16] found that the anisotropy of shale elastic parameters is significantly affected by the amount of clay and organic content as well as the shale fabric. [17] firstly deduced the analytical solution for the elastic problem of surface loading on transversely isotropic bodies. [18] applied the anisotropic mechanics model to the wellbore stability analysis. The results show that the anisotropy of rock causes an obviously different stress distribution around the wellbore from that calculated by the isotropic model. [19] established the anisotropic Mises-Schleicher criterion (AMS) by introducing a four-order anisotropic tensor, to describe the failure characteristics of transversely anisotropic rocks under compression or tensile conditions. [20] improved the Cam Clay model to describe the deformation characteristics of shale, as well as the strength weakening process after shale failure, by introducing an orthotropic elasticity and an orthotropic pressure-dependent yield surface. Based on the anisotropic strength criterion of McLamore and Gray for rock and soil, [21, 22] improved the isotropic Drucker-Prager strength criterion as an anisotropic elastoplastic constitutive model for layered rock, which can effectively reflect the deformation and strength characteristics of layered rocks. [23] presented a modified Drucker-Prager yield criterion for the transversely isotropic geomaterials by considering the anisotropy of the friction angle and the dilation angle. [24] introduced an anisotropic parameter into the Hoek-Brown failure criterion and experimentally studied the relation between the anisotropy parameter and the degree of anisotropy. There are many other methods to describe the anisotropic characteristics of material deformation or strength by introducing anisotropic parameters [25–27], which will not be detailed here. With the development of computer technology, a numerical simulation method provides a new way to simulate anisotropic materials. Sainsbury used the 3DEC discrete element code to simulate the anisotropic characteristics of the real rock mass by including the joint elements, and the effects of joints on the rock deformation and strength characteristics were also investigated [28]. [29] established a discrete fracture model to characterize the heterogeneity of shale reservoirs with interlaced distributed natural fractures. In addition to describing the fluid flow in complex fracture networks, the model can also deal with the nonuniform distribution of stress and the anisotropy of the strength caused by the opening and shearing of natural fractures. Therefore, the discrete fracture model provides an alternative way of characterization of shale mechanics. [30] used bonded-particle discrete element modeling with embedded smooth joints to simulate the mechanical behavior of transversely isotropic rock and demonstrate the effectiveness of the new method in modeling the equivalent anisotropic medium.
The above research puts forward different methods to describe anisotropy materials, which have their own characteristics and advantages. The failure criterion established by [19] can reflect the strength characteristics of the transversely anisotropic material. The defect is that massive parameters need to be determined and the rock structural discontinuity cannot be characterized. Crook indirectly characterized strength anisotropy of materials using the method of the equivalent stress, which also has the defect of too many material parameters, while the discrete fracture model provides an effective method that can effectively represent the effect of beddings on rock mechanical and physical properties [29]. In this paper, a continuum damage-based discrete fracture network method (D-DFN) is proposed to describe the structural discontinuity and mechanical anisotropy of laminated shales, where the stiffness and strength evolution of shale after failure are described based on the continuum damage theory, and the discrete fracture network modeling method is used to describe the bedding structure of shales. The model parameters are calibrated by experiment, and the validity of the D-DFN method is proved by comparison between the simulation results and experiment results.
2. Compression Tests on Shales
In this paper, the experimental investigation is performed on Longmaxi shales. This is a silicic, fine-grained, black shale with an average of 21% clay content. During diagenesis of shales, the external force (such as tectonic movement) can cause the opening and propagation of shale joints, forming natural fracture networks [31]. Shown in Figure 1 is the distribution of natural fractures in Longmaxi shales with different scales. At least on the size of Figure 1(b), shales should be described by a discrete fracture network (DFN) model [32]. That is, shales should be considered a combination of the natural fractures and the shale matrices cut by fractures, due to the large differences in mechanical and physical properties between natural fractures and matrices. Even if the rock sample size is reduced to the size of the standard core while the core contains weak interfaces (Figure 1(c)), the results of mechanical tests carried out on this core reflect the combined effect of the interfaces and the matrices. Therefore, a new method is needed to explain conventional mechanical test results carried out on shales.
(a) Shale
(b) A shale sample for physical experiment
(c) A standard core with beddings
Standard cores with diameters of 2.5 cm and lengths of 5 cm (Figure 1(c)) are used in the tests. In order to reduce the discreteness of the test results, all the standard cores are taken from the same shale block. The bedding orientation in triaxial tests is denoted by the angle between the loading direction and the bedding planes (Figure 2). Two sets of core compression tests have been performed to investigate the effects of confining pressure and bedding orientations on rock mechanical properties: (1)A series of drained triaxial compression tests with the bedding planes orientated normal to the axis of the specimen () and confining pressures of 0, 10, 20, 30, 40, and 50 MPa(2)A series of drained triaxial compression tests with confining pressures of 40 MPa and different bedding plane orientations (, 30°, 45°, 60°, and 90°)
The conventional triaxial compression tests with different confining pressures show (Figures 3 and 4) that (1)the cores mainly show shear failure except at lower confining pressures (0 MPa and 10 MPa) (Figure 3)(2)with the increase of confining pressure, the elastic modulus of shale is not changed much (average 31 GPa), but the strength of shale is increased evidently (Figure 4)
(a) Deviatoric stress-volumetric strain
(b) Deviatoric stress-axial/radial strain
The conventional triaxial compression tests with different bedding orientations (Figures 5–8) show that (1)the cores mainly show shear failure along bedding planes when equals 30 or 45. When is 60, besides the shear failure along the bedding plane, shear failure also occurs to shale matrices. When equals 0 or 90, mainly shear failure occurs to the shale matrix(2)as increases from 0 to 90, the tested elastic modulus decreased (Figure 7). The anisotropic degree of the modulus of elasticity () can be calculated as 1.10. According to Worotnicki’s statistics of 200 sets of core compression test results [34], the shale samples of this experiment can be classified as a kind of rock with smaller elastic anisotropy(3)with the increase of from 0 to 90, the ultimate deviatoric stress () first decreases and then increases, and the minimum value is obtained at 45 (Figure 8). The strength anisotropy () is calculated to be 1.27. It shows that the strength anisotropy of shale in this experiment is slightly higher than that of elastic modulus due to the bedding structure
(a) Deviatoric stress-volumetric strain
(b) Deviatoric stress-axial/radial strain
3. Simulation Method
In this paper, the mechanical property of shales is studied by combining experimental and numerical methods. Based on current laboratory testing conditions, the following assumptions were made in numerical simulation, as well as for simplified calculation. (1)The elastic deformation stage before shale failure is described by the anisotropic elastic constitutive model. The damage theory is used to explain the stiffness and strength degradation after shale failure, and the damage evolution process is assumed to be isotropic(2)The shale bedding structure is described by the discrete fracture network (DFN) modeling method. The mechanical and physical difference between shale beddings and masses was taken into account by setting different initial damage variables to the bedding elements and the matrix elements (Figure 9), where a nonzero initial damage is set to bedding elements to show initial lower elastic modulus and strength of beddings compared with initially nondamaged matrix elements
(a)
(b)
3.1. Isotropic Elastic Constitutive Model
Because of the geological deposition, shales pose a distinct layered structure, and the transversely isotropic assumption is usually carried out to simplify the calculation. According to the theory of elasticity, the expression of the strain-stress relationship of the orthotropic body is as follows [17]:
In Equation (1), there are nine independent parameters in the orthotropic model, including Young’s modulus, , , and in three orthogonal directions, three Poisson ratios, , , and , and three shear modulus, , , and .
Assuming that the 1-2 plane is an isotropic plane, then , , and , in which and represent the transverse and normal direction, respectively, and the strain-stress expression of the transversely isotropic body is expressed as follows [18–20]: where . According to the Saint-Venant principle [35], can be approximated as . Therefore, four parameters are needed for describing the transversely isotropic model of shale, which are , , , and .
3.2. Hyperbolic Drucker-Prager Plastic Model
The Drucker-Prager model is widely used to characterize the deformation and strength of geotechnical materials. The yield surface equation can be expressed as a linear relationship between the equivalent stress and the mean stress [36], as shown in Equation (3). In the three-dimensional principal stress space, the shape of the yield surface is an open circular cone. Drucker-Prager of the linear form itself has shortcomings in predicting failure in low confining pressures and tensile failure, which needs to be improved according to specific engineering problems: where is the yield function, is the equivalent stress, is the mean stress, and and are the friction angle and cohesion in space, respectively.
Here, to predict failure of shales under both compression and tension stress, the modified Drucker-Prager criterion of hyperbolic form was employed, which can be expressed as [37] where is the initial hydrostatic tension strength. Figure 10 shows the comparison of the linear form and hyperbolic form of the Drucker-Prager criterion. Later, it can be found that the curve of hyperbolic form fits the experimental data better.
In the plastic stage, the plastic potential function of shales is defined as where is the plastic potential function, , and is the dilation angle-related parameter.
The evolution of the plastic flow is defined by a nonassociative flow rule (): where is the plastic strain and is the plastic multiplier.
3.3. Damage Constitutive Model
When the stress reaches the peak strength, the shale will soften, exhibiting a transformation from a homogeneous strain field to a heterogeneous strain field with the localized region of intense strain due to the stiffness degradation in the fractured zone. Damage mechanics is often used to describe the failure process of concrete or rock materials [38]. A damage variable related to plastic deformation is implied in the model to describe the softening stage. To eliminate the effect of element size on the calculation results, the damage variable is stated by the stress-displacement relationship rather than the stress-strain relationship [39]. Considering the brittle failure characteristics of shale, it is assumed that the relationship between the damage variable and the postfailure equivalent plastic deformation is suited to the first-order exponential decay function (Equation (7)), and the related parameter can be calibrated by results of the triaxial tests. In order to facilitate programming, the relationship between damage variable and equivalent plastic strain is represented by a piecewise linear function, as shown in Figure 11: where is the base of the natural logarithm, is the postfailure equivalent plastic displacement, and is the material parameter, which reflects the evolution rate of with the element deformation.
Since shale is a brittle material, the elastic damage constitutive relation of the element under uniaxial compressive stress and tensile stress can be further simplified, as shown in Figure 12. When the stress of the element satisfies the strength criterion (such as the Drucker and Prager criterion, expressed in Equation (4)), the element begins to fail. In elastic damage mechanics, with the development of damage, the elastic modulus of the element may gradually degrade. The elastic modulus of the damaged element is defined as follows [40]: where is the damage variable, which is a scaler in this paper; is the initial elastic modulus; and is the elastic modulus when the element is damaged.
Tang et al. summarized the elastic damage constitutive model, where the element will be destroyed when the strain of the element exceeds a certain value under tension or compression conditions [40]. Since the description of damage is based on the stress-strain relation, Tang et al.’s modeling results are sensitive to the element size. Therefore, the size of all elements in Tang et al.’s models is uniform and close to the size of rock crystal particles [40]. Hillerborg et al. took the energy required for fracture opening as a material parameter, based on the theory of brittle fracture mechanics. In such a method, the fracture behavior of elements is described by the stress-displacement relation rather than the stress-strain relation, which can reduce the influence of mesh size on the calculation results [39].
In this paper, the definition of damage variable is based on the stress-displacement relation. Figure 12 shows the stress-displacement constitutive relation of an element under uniaxial tension and compression conditions.
When an element is subjected to tensile failure, taking the uniaxial tension condition as an example, the expression of the damage variable is as follows: where is the elongation of the element under tension stress, which is the product of element strain and element characteristic length ; is the element elongation when the peak stress is reached; and is the element elongation at the residual stress . When , the element is completely damaged and a smeared fracture will be formed in the element.
When a compressive failure occurs to an element, taking the uniaxial compression condition as an example, the expression of the damage variable is as follows: where is the change of element length under compression stress. Other parameters in the equation are defined in Figure 12. Compared with tension failure, the damage cannot reach 1 under compression conditions. That is because rock will shear slide along the broken surface after compression failure, thus maintaining certain residual stress.
4. Calibration of the Numerical Model
The elastic properties can be determined directly from the linear part of the stress-strain relation curves of compression tests on cores with different bedding orientations (Figure 6). The shale studied in the paper is assumed transverse isotropic, and the calibrated elastic parameters are summarized in Table 1.
Compression tests on cores with bedding planes orientated normal to the axis of the specimen () and different confining pressures (0, 10, 20, 30, 40, and 50 MPa) (Figure 4) are selected to calibrate the parameters in the hyperbolic Drucker-Prager model (Equation (4)). Also, the tension strength (about 14 MPa, obtained from Brazil split tests) is needed to determine the intersection of the failure curve and the horizontal axis (Figure 10). The mean pressures and respective deviated stresses at the failure points under different test conditions are collected in Table 2. The yield function [19] is used to fit the data (Figure 13), and the parameters, , , and , are determined as shown in Table 3.
After yielding, there is a hardening stage before shale reaching the peak strength in each test (Figure 4). The hardening behavior can be described by one of the stress-strain curves in Figure 4. The principle is that the hardening curve can represent the average trend of the relationship between the equivalent stress increment and the equivalent plastic strain increment under all confining pressure conditions. Here, we choose the plastic section of the stress-strain curve with confining pressure of 30 MPa to determine the hardening stage.
The simulation model is built using ABAQUS as a platform, and a user subroutine is implemented in the model to relate rock element parameters (e.g., elastic modulus and strength) to the damage variable [37]. Figure 14 shows the algorithm flowchart for the th calculation step.
According to the geometry of the standard rock core, a 3D finite element model discretized by hexahedral elements was developed (Figure 9). In order to take account of computational efficiency and accuracy, the average size of cells is optimized to 1.5 mm, and C3D8R is chosen as the element type. The displacement loading method is adopted (Figure 9(a)).
In order to simulate the influence of weak beddings on shale strength, the finite element model is divided into matrix element sets and bedding element sets (Figure 9(b)). As mentioned before, the mechanical difference between shale beddings and masses was taken into account by setting a nonzero initial damage variable to bedding elements to show initial lower elastic modulus and strength compared with initially nondamaged matrix elements. However, it is difficult to determine the initial damage strength of bedding only by experiment. Here, the paper combines numerical simulation and experiments to solve the problem. Specifically, establishing a core model with of 30° to simulate, different initial damage variables (such as 0.1, 0.3, 0.5, 0.7, and 0.9) are assigned to the bedding elements every time, and the corresponding stress-strain curves are obtained. Then, the initial damage variable of bedding can be determined by comparing the simulation strength of core with of 30 to compression test results. Finally, good agreement can be found when the initial damage variable of beddings is set to 0.3. Therefore, the initial damage variable of bedding in all numerical models is set as 0.3 for the shale samples used in this paper. All material parameters are summarized in Table 4.
5. Modeling Results
Figure 15 shows the initial and postfailure damage distribution of models with different bedding orientations. It can be found that when is 30 or 45, the main failure mode is shearing along bedding planes (Figures 15(b) and 15(c)). When is 0 or 90, the shear failure of the matrices is occurring (Figures 15(a) and 15(e)). With changing from 45 to 60, the normal stress acting on beddings will increase, resulting in the greater shear strength of bedding planes. When is 60°, the shear failure will occur both in matrices and along the bedding planes (Figure 15(d)). The numerical simulation results are in good agreement with the experimental results (Figure 5).
Figure 16 shows the stress-strain curves obtained from numerical simulations with different bedding orientations. The variation trends of elastic modulus and strength of shale over are also compared well with the experimental results (Figures 17 and 18). It is proved that the D-DFN model is effective in simulating the mechanical properties of the Longmaxi black shale.
(a) Simulation results
(b) The contrast between simulation and test results
6. Conclusion
In view of the bedding structure of the Longmaxi black shale in Sichuan basin, a finite element model based on continuous damage theory and discrete fracture network modeling method (D-DFN method) is developed to describe the shale’s structural discontinuity and the mechanical anisotropy.
In the model, it is assumed that the elastic deformation of shales satisfies the transversely isotropic deformation law. After failure, the isotropic damage model is used to describe the strain softening trend. The discrete fracture modeling method is used to characterize the shale matrices and beddings by setting different initial damage variables to bedding elements and matrix elements. By comparing with the triaxial test results, the validity of the D-DFN method in simulating the deformation and strength anisotropy of the shale is verified.
The model assumes that the strength of the shale matrix is isotropic. Therefore, the D-DFN method proposed in this work is suitable for shale masses with low matrix strength anisotropy. A further study is needed for the characterization of strength anisotropy of shale matrix, to expand the application of the model.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
This work is supported by the Shandong Provincial Natural Science Foundation (Nos. ZR2019ZD14 and ZR2020QE120).