# Characterization of Discontinuity and Mechanical Anisotropy of Shale Based on Continuum Damage Mechanics

Authors: Qinglin Shan; Peng Yan; Hengjie Luan; Yujing Jiang; Sunhao Zhang
Journal: Geofluids
 (2021)

---

## 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.


---

## Body

## 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.Figure 1
Distributions of natural fractures in shales [33].
(a)
Shale(b)
A shale sample for physical experiment(c)
A standard core with beddingsStandard cores with diameters of 2.5 cm and lengths of 5 cm (Figure1(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 (β=90°) 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 (β=0°, 30°, 45°, 60°, and 90°)Figure 2
Illustration of coring directions of shale samples.The conventional triaxial compression tests with different confining pressures show (Figures3 and 4) that
(1)
the cores mainly show shear failure except at lower confining pressures (0 MPa and 10 MPa) (Figure3)(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 (Figure4)Figure 3
Broken cores after compression with different confining pressures.Figure 4
Stress-strain curves of cores under varying confining pressures in triaxial compression tests.
(a)
Deviatoric stress-volumetric strain(b)
Deviatoric stress-axial/radial strainThe conventional triaxial compression tests with different bedding orientations (Figures5–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 (Emax/Emin) 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 (σp) first decreases and then increases, and the minimum value is obtained at 45° (Figure 8). The strength anisotropy (σpmax/σpmin) 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 structureFigure 5
Failure modes of cores after experiments.Figure 6
Stress-strain curves of cores with different bedding orientations.
(a)
Deviatoric stress-volumetric strain(b)
Deviatoric stress-axial/radial strainFigure 7
Variation of elastic modulus overβ.Figure 8
Variation of peak strength overβ.


## 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 (Figure9), 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 elementsFigure 9
Finite element model configuration and load method for triaxial tests.
(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]:
(1)ε11ε22ε33γ12γ13γ23=1E1−ν21E2−ν31E3−ν12E11E2−ν32E3−ν13E1−ν23E21E30000000000000000001G120001G130001G23σ11σ22σ33σ12σ13σ23.In Equation (1), there are nine independent parameters in the orthotropic model, including Young’s modulus, E1, E2, and E3 in three orthogonal directions, three Poisson ratios, ν12, ν13, and ν23, and three shear modulus, G12, G13, and G23.Assuming that the 1-2 plane is an isotropic plane, thenE1=E2=Ep, ν13=ν23=νt, and G13=G23=Gt, in which p and t represent the transverse and normal direction, respectively, and the strain-stress expression of the transversely isotropic body is expressed as follows [18–20]:
(2)ε11ε22ε33γ12γ13γ23=1Ep−νpEp−νtEt−νpEp1Ep−νtEt−νtEp−νtEp1Et0000000000000000001Gp0001Gt0001Gtσ11σ22σ33σ12σ13σ23,where Gp=Ep/21+νp. According to the Saint-Venant principle [35], Gt can be approximated as Gt=EpEt/Ep+Et+2Epνt. Therefore, four parameters are needed for describing the transversely isotropic model of shale, which are Ep, Et, νp, and νt.


### 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:
(3)F=q−ptanθ−d=0,where F is the yield function, q is the equivalent stress, p is the mean stress, and θ and d are the friction angle and cohesion in p~q 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]
(4)F=d−pt0tanθ2+q2−ptanθ−d=0,where pt0 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.Figure 10
Comparison of the linear form and hyperbolic form of Drucker-Prager criterion in thep~q plane.In the plastic stage, the plastic potential function of shales is defined as(5)G=ltanψtanθ2+q2−ptanψ,where G is the plastic potential function, l=d−pt0tanθ, and ψ is the dilation angle-related parameter.The evolution of the plastic flow is defined by a nonassociative flow rule (ψ≠θ):
(6)ε̇p=λ̇∂G∂σ,where εp 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 D 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 D 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 a 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:
(7)D=1e−1/a−1e−up/a−1e−1/a−1,where e is the base of the natural logarithm, up is the postfailure equivalent plastic displacement, and a is the material parameter, which reflects the evolution rate of D with the element deformation.Figure 11
Relationship between damage variable and postfailure equivalent plastic displacement.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 Figure12. 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]:
(8)E=1−DE0,where D is the damage variable, which is a scaler in this paper; E0 is the initial elastic modulus; and E is the elastic modulus when the element is damaged.Figure 12
Stress-displacement relationship of an element subject to uniaxial tensile stress (at lower left) and uniaxial compressive stress (at upper right) [33].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 variableD 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:(9)D=0,u3≤ut0,u3−ut0utr−ut0·1−σtrσt0,ut0<u3≤utr,1−utu−u3utu−utr·σtrσt0,utr<u3≤utu,1,u3>utu,where u3 is the elongation of the element under tension stress, which is the product of element strain and element characteristic length l; ut0 is the element elongation when the peak stress σt0 is reached; and utr is the element elongation at the residual stress σtr. When u3>utu, 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:(10)D=0,u1≤uc0,u1−uc0ucr−uc01−σcrσc0,uc0<u1≤ucr,1−σcrσc0,u1>ucr,where u1 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.


## 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]:
(1)ε11ε22ε33γ12γ13γ23=1E1−ν21E2−ν31E3−ν12E11E2−ν32E3−ν13E1−ν23E21E30000000000000000001G120001G130001G23σ11σ22σ33σ12σ13σ23.In Equation (1), there are nine independent parameters in the orthotropic model, including Young’s modulus, E1, E2, and E3 in three orthogonal directions, three Poisson ratios, ν12, ν13, and ν23, and three shear modulus, G12, G13, and G23.Assuming that the 1-2 plane is an isotropic plane, thenE1=E2=Ep, ν13=ν23=νt, and G13=G23=Gt, in which p and t represent the transverse and normal direction, respectively, and the strain-stress expression of the transversely isotropic body is expressed as follows [18–20]:
(2)ε11ε22ε33γ12γ13γ23=1Ep−νpEp−νtEt−νpEp1Ep−νtEt−νtEp−νtEp1Et0000000000000000001Gp0001Gt0001Gtσ11σ22σ33σ12σ13σ23,where Gp=Ep/21+νp. According to the Saint-Venant principle [35], Gt can be approximated as Gt=EpEt/Ep+Et+2Epνt. Therefore, four parameters are needed for describing the transversely isotropic model of shale, which are Ep, Et, νp, and νt.


## 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:
(3)F=q−ptanθ−d=0,where F is the yield function, q is the equivalent stress, p is the mean stress, and θ and d are the friction angle and cohesion in p~q 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]
(4)F=d−pt0tanθ2+q2−ptanθ−d=0,where pt0 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.Figure 10
Comparison of the linear form and hyperbolic form of Drucker-Prager criterion in thep~q plane.In the plastic stage, the plastic potential function of shales is defined as(5)G=ltanψtanθ2+q2−ptanψ,where G is the plastic potential function, l=d−pt0tanθ, and ψ is the dilation angle-related parameter.The evolution of the plastic flow is defined by a nonassociative flow rule (ψ≠θ):
(6)ε̇p=λ̇∂G∂σ,where εp 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 D 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 D 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 a 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:
(7)D=1e−1/a−1e−up/a−1e−1/a−1,where e is the base of the natural logarithm, up is the postfailure equivalent plastic displacement, and a is the material parameter, which reflects the evolution rate of D with the element deformation.Figure 11
Relationship between damage variable and postfailure equivalent plastic displacement.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 Figure12. 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]:
(8)E=1−DE0,where D is the damage variable, which is a scaler in this paper; E0 is the initial elastic modulus; and E is the elastic modulus when the element is damaged.Figure 12
Stress-displacement relationship of an element subject to uniaxial tensile stress (at lower left) and uniaxial compressive stress (at upper right) [33].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 variableD 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:(9)D=0,u3≤ut0,u3−ut0utr−ut0·1−σtrσt0,ut0<u3≤utr,1−utu−u3utu−utr·σtrσt0,utr<u3≤utu,1,u3>utu,where u3 is the elongation of the element under tension stress, which is the product of element strain and element characteristic length l; ut0 is the element elongation when the peak stress σt0 is reached; and utr is the element elongation at the residual stress σtr. When u3>utu, 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:(10)D=0,u1≤uc0,u1−uc0ucr−uc01−σcrσc0,uc0<u1≤ucr,1−σcrσc0,u1>ucr,where u1 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 (Figure6). The shale studied in the paper is assumed transverse isotropic, and the calibrated elastic parameters are summarized in Table 1.Table 1
Elastic parameters of transversely isotropic shale.
Ep (GPa)Et (GPa)Gp (GPa)Gt (GPa)νpνt35.4932.4114.7914.010.210.22Compression tests on cores with bedding planes orientated normal to the axis of the specimen (β=90°) 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 p and respective deviated stresses q 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, pt0, d, and θ, are determined as shown in Table 3.Table 2
p and q under different test conditions.
p (MPa)-4.6749.2377.67101.67133.96147.54166.90q (MPa)14.10146.96202.86245.02311.89322.68350.72Figure 13
Fitting of the yield surface in thep~q plane.Table 3
Parameters defining yield surface inp~q plane.
ParametersValuespt0 (MPa)5.05d (MPa)88.90θ (°)58.98After yielding, there is a hardening stage before shale reaching the peak strength in each test (Figure4). 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 ith calculation step.Figure 14
Flowchart for numerical simulation.According to the geometry of the standard rock core, a 3D finite element model discretized by hexahedral elements was developed (Figure9). 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 (Figure9(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.Table 4
Material parameters of shale used in the simulations.
Parameter itemsValuesThe initial out-of-plane elastic modulus (Et) (GPa)32.41The initial in-plane elastic modulus (Ep) (GPa)35.49The in-plane Poisson ratio (μp)0.21The out-of-plane Poisson ratio (μt)0.22The initial tensile strength (σt0) (MPa)14The residual tensile strength (σtr) (MPa)2The initial compressive strength (σc0) (MPa)150The residual compressive strength (σcr) (MPa)30The fraction angle relative parameter (θ) (°)58.98The initial mean tensile strength (pt0) (MPa)5.05The dilation angle relative parameter (ψ) (°)48The material parameter (a) (mm)0.003The initial bedding damage variable (D0)0.3


## 5. Modeling Results

Figure15 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 15
Simulations of damage distributions of cores under different bedding orientations.Figure16 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.Figure 16
Simulation results of stress-strain relationship of cores with different bedding orientations and the contrast with test results.
(a)
Simulation results(b)
The contrast between simulation and test resultsFigure 17
The contrast between simulation and test results on the variation of elastic modulus overβ.Figure 18
The contrast between simulation and test results on the variation of shale strength overβ.


## 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.



---

*Source: 5714547-2021-08-14.xml*