The anisotropic mechanical characteristics of layered rocks under numerical simulation

Layered rocks pose the challenge of wellbore stability in drilling engineering because of the anisotropic mechanical properties caused by the distinct weak planes. To understand the significant anisotropy of layered rocks in real formation condition, true triaxial compression tests are conducted by numerical simulation in this study. It is revealed that the mechanical responses of layered rocks are either controlled by the rock matrix or dominated by the weak plane and exhibit three different types associated with the orientations of the weak plane (including the dip direction α and dip angle β). When the orientations of the weak plane are α = 0°–90° and β = 0°, 60°–90°, the failure and strength properties of layered rocks depend entirely on the rock matrix, classified to the first type. Whereas the layered rocks with angle α ≤ 45° and β = 15°–45° fail by slipping failure along the weak plane, the relationship curves of rock strength versus the intermediate principal stress (σ2) are downward convex parabolas. In the last type, the mechanical behaviors of layered rocks with α > 45° and β = 15°–45°, involved in the changes of failure mode and the strength curve, are complex. Besides, the limitation of the simulation is discussed, and further studies on layered rocks are essential.


Introduction
The mechanical characteristics of rock are the basis of engineering design and construction in oil and gas exploration. As a special rock material containing a (or a set of parallel) weak plane(s), the layered rocks exhibit anisotropic mechanical characteristics related to the orientation of weak planes, resulting in the challenge of mechanically induced wellbore instability in drilling process. Thus, it is essential to thoroughly understand the mechanical responses of layered rocks. As early as 1959, Jaeger (1959) had carried out compression tests on layered rocks with artificial joint surfaces. Subsequently, great efforts have been made to investigate the mechanical properties of various layered rocks by using conventional triaxial tests and other methods. Based on experimental results, there are well comprehensions of anisotropic mechanical properties under laboratory conditions for the natural layered rocks, such as shale (Donath 1961;Chenevert and Gatlin 1965;McLamore and Gray 1967;Kwaśniewski and Nguyen 1988;Niandou et al. 1997;Ajalloeian and Lashkaripour 2000;Cho et al. 2012;Geng et al. 2016), schist (McCabe andKoerner 1975;Behrestaghi et al. 1996;Nasseri et al. 1997;Zhang et al. 2011) and phyllite (Donath 1972;Ramamurthy 1993;Fereidooni et al. 2016).
The anisotropy of mechanical characteristics of layered rocks consists of the rock strength, failure mode and deformation, associated with the dip angle β between the axial loading direction and the normal of the weak plane. The anisotropy of strength is usually the primary concern among them. Based on the results of conventional triaxial compression tests on shale and slate, Donath (1961) revealed that the rock strength curve of layered rocks with respect to the dip angle is an upward concave parabola. As the angle β increases from 0° to 90°, the rock strength decreases first and then increases for natural and artificial layered rocks. Consequently, it is widely recognized that the strength curve shows a U-shaped line, where the maximum strength occurs at β = 0° or 90° and the minimum probably corresponds to the angle β of 30° (Fereidooni et al. 2016). Further, Behrestaghi et al. (1996) defined the shape of the relationship 1 3 curve between the UCS (uniaxial compressive strength) and β as "the type of anisotropy," with the classifications of "Shoulder type," "U-shaped type" and "Wave type." The failure and deformation of layered rocks are also vitally important for the engineering projects. According to the results of compression tests on layered rocks, there are two types of failure modes: the failure of the rock matrix and the failure of the weak plane (Tien et al. 2006). When the angle β belongs to the intermediate angles (around 15° ~ 45°), the sliding failure occurs along weak planes, which results in a relatively lower rock strength. The failure modes of layered rocks at other dip angles are mainly shear failure across weak planes, occurring at a stress determined by the matrix. The anisotropy of deformation is also related to the dip angle, characterized by the differences of elastic modulus and Poisson's ratio at different angles.
Since around the same time, a number of true triaxial tests have been conducted to investigate the effect of principal stresses (especially the intermediate principal stress σ 2 ) on rock mechanical characteristics (Mogi 2006). These experiments made by Mogi (1967Mogi ( , 1971, Takahashi and Koide (1989), and others studied a variety of homogeneous rocks, and the results revealed the strengthening effect of the intermediate principal stress on rock strength. With the discovery of it, the σ 2 effect on layered rock has also attracted the attention of scholars. Akai (1971) carried out triaxial compression tests on the chlorite schist under a constant average principal stress. Then, Kwaśniewski and Mogi (2000) investigated mechanical characteristics of the Chichibu greenschist under true triaxial condition with the consideration of the dip direction (α) of the weak plane. Lately, Lin et al. (2013) explored the σ 2 effects on the low-strength layered rocks at the angle α = 0° or 90° by using numerical simulation. The previous studies have proved that the orientations of weak planes (including the dip direction α and dip angle β) have significant influences on the strength and failure of layered rocks. However, with a few experiments on layered rocks, there are still many unknowns about their mechanical properties under true triaxial compression, such as the coupled effect of the angles α and β (Mogi 2006).
In this paper, the numerical simulation method is used to explore the anisotropic mechanical properties of layered rocks under true triaxial compression. With the full consideration of weak plane orientations and stress conditions, the corresponding research on layered rocks has been systematically completed. Thus, the further understanding of the mechanical responses of layered rocks is obtained from the study, which can promote the development of the related experimental and theoretical studies. Future research on layered rocks, subjected to the stress states of real formation, will make it possible to solve engineering problems such as the mechanically induced wellbore instability in shale formation.

Modeling
In this study, the sample models of layered rocks were built by the same method adopted by Lin et al. (2013) and then imported into FLAC3D for calculation and analysis. At present, cubic and rectangular prismatic samples are commonly used in true triaxial compression tests (Kwaśniewski et al. 2012). The rock models adopted in this study are standard rectangular prisms with the size of 50 mm × 50 mm × 100 mm. As mentioned before, the orientations of weak planes consist of the dip angle β and dip direction α under true triaxial condition. The dip direction α is the angle between the trace line of the weak surface in the horizontal plane (the plane formed by σ 2 and σ 3 ) and the intermediate principal stress σ 2 , and the dip angle β is the angle between the weak surface and the maximum principal stress σ 1 , as shown in Fig. 1. Both α and β vary from 0° to 90°, with the change gradient of 15°. Therefore, there are seven rock models with different angle α at the same angle β, taking the specimens with the angle β = 0° as examples, as shown in Fig. 2a-g. However, the weak plane is parallel to the horizontal plane at β = 90°, resulting in only one rock model, as shown in Fig. 2h. Thus, there is an assumption that layered rocks with the angle β = 90° and different dip directions exist, which are same as each other. In the rock models, the weak planes are represented by the white grid surfaces with no thickness, which divide the rock matrix into blocks (the equal-thickness blue part) with a thickness of 10 mm.
The Drucker-Prager plasticity model in FLAC3D, characterized by the Drucker-Prager criterion and a tensile limit, was selected as the constitutive model of the rock matrix (Itasca Consulting Group 2002). For this model, the strength where σ are τ are the octahedral normal stress and shear stress (Alejano and Bobet 2012), respectively; q and k are constants related to the material, determined by cohesion and internal friction angle using the following equations: (1) where c is the cohesion of the rock matrix, MPa; ϕ is the internal friction angle of the rock matrix, °. The tensile failure is defined by the tensile failure criterion as where σ t is the tensile strength of the rock matrix, MPa.
In addition, this model is an ideal plastic model and follows the basic assumption of plasticity theory, that is, the total strain increment is composed of elastic and plastic parts, but only the elastic part will contribute to the stress increment. In elastic deformation stage, stress increment and strain increment satisfy Hooke's law: where Δγ e and Δε e are the increments of shear strain and volume strain, Δτ and Δσ are the increments of tangential stress and normal force, respectively, G is the shear modulus, and K is the bulk modulus. (4) On the other hand, the weak planes of layered rock are controlled by the interfaces, which are characterized by the Coulomb shear strength criterion (Itasca Consulting Group 2002). If the shear yield of the weak plane happens, the required shear force is limited by the following relation: where c w is the cohesion of the interface, ϕ w is the friction angle of the interface, F n is the normal force of the interface, A is the interface area, and p is the pore pressure (not considered in this study).
The properties of the layered rock are cited from the results of conventional triaxial compression tests on shale: c = 28.00 MPa, ϕ = 38.50°; c w = 15.90 MPa, ϕ w = 30.00°. Following the theory and background of the FLAC3D user manual (Itasca Consulting Group 2002), the parameters of the rock matrix and weak plane are shown in Table 1 and  Table 2, respectively. The intermediate principal stress σ 2 and the minimum principal stress σ 3 , as stress loads, are both uniformly applied to the corresponding side surfaces of the rock sample, as shown in Fig. 1 and Fig. 2. And the two rock surfaces perpendicular to the Z-axis direction are simultaneously subjected to the displacement loads with the same magnitude and opposite directions. During the simulation process, the stress loads that ignore the influence of stress increments are first performed, and then the displacement loads with the magnitude of 4.5 × 10 -5 mm/step run for predetermined steps.

True triaxial tests of σ 2 = σ 3
At first, the true triaxial compression experiments with the stress state of σ 2 = σ 3 were performed on all the rock samples. Due to the same magnitude, the stress loads (σ 2 and σ 3 ) were referred to as the confining pressure, which increases from 0 to 30 MPa with the gradient of 10 MPa. Based on the basic assumptions of plasticity theory, the stress-strain curve of rock sample in the numerical simulation turns into two phases: The first is the linear elastic phase, where the stress increases linearly with the increasing strain, whereas the change of stress in the second phase (called as the plastic phase) is relatively small, as shown in Fig. 3. The intersection of the two phases is the critical point of deformation with the rock failure. Thus, the peak strength can be obtained from the stress-strain curve when the Z-axis strain of rock reaches this point.
For layered rocks with dip direction α = 0°, the relationship curves between rock strength and dip angle β are in the shape of "Shoulder type" under different confining pressures, as shown in Fig. 4. As β increases from 0° to 90°, the rock strength decreases first and then increases, which is consistent with the result of conventional triaxial experiment. Taking two rock samples (one with α = 0° and β = 0° and another one with α = 0° and β = 30°) as examples, the errors between the simulated and experimental value of rock strength are small, and the curves of rock strength versus confining pressure all follow a linear law, as shown in Fig. 5. These strength characteristics are also applicable to layered rocks with other dip directions.
The angle β plays a significant role in the strength of layered rocks under the condition of σ 2 = σ 3 , whereas the variations of rock strength with angle α are almost negligible for  . 3 The stress-strain curve of some layered rocks; a α = 0°, β = 0°, and σ 2 = σ 3 = 10 MPa, b α = 0°, β = 30°, and σ 2 = σ 3 = 10 MPa (Z-axis: the Z-axis strain of specimen, X-axis: the X-axis strain of specimen, Y-axis: the Y-axis strain of specimen. The latter two overlap each other) 1 3 layered rocks with the same dip angle. Taking the two series of rock samples with β = 0° and β = 30° as examples, the rock strength curves with respect to the angle α are nearly horizontal straight lines at different confining pressures, as shown in Fig. 6. Based on the assumption about the rock models, it is also suitable for layered rocks with β = 90°. In terms of the failure characteristics of layered rocks, the simulation results are in agreement with the experiment results. The failure modes of layered rocks in numerical simulation can also be classified into two types, including the shear failure of the rock matrix and slipping failure along the weak plane. In this paper, the rock failure is illustrated by the Z-direction displacement diagram after simulated due to the limit of software. When α = 0° ~ 90° and β = 0°, 60° ~ 90°, the layered rocks fail by the shear failure of matrix, characterized by the uniform distribution of displacement in the Z-axis direction, as shown in Fig. 7. And the slipping failure along the weak plane occurs at α = 0° ~ 90° and β = 15° ~ 45°. As a result of it, the rock sample disintegrates into two parts along a weak plane, which are equal in magnitude but opposite in direction, as shown in Fig. 8.

True triaxial tests of σ 2 > σ 3
Based on the true triaxial tests mentioned in Sect. 3.1, a series of true triaxial tests with the stress load of σ 2 > σ 3 are further carried out on layered rocks. The minimum principal stress σ 3 , same as before, can be one  1 3 σ 1 -σ 2 curves of these rock samples are roughly the same as each other at the same stress σ 3 . Taking some rock samples as examples, the rock strength of each specimen increases almost linearly with the σ 2 , as shown in Fig. 9. Correspondingly, these layered rocks all have been broken by the shear failure of the rock matrix and characterized by the similar Z-direction displacement in Fig. 7. Therefore, these mechanical characteristics unassociated with the weak plane are summarized as the first type.
For the layered rocks with α = 0°-45° and β = 15°-45°, the mechanical responses are mainly affected by the weak plane. With the increase in the σ 2 , the rock strength first increases and then decreases, as shown in Fig. 10. The strength curves of these rocks are convex parabolas with opening downward. Under a certain stress condition, the initial slope and peak value of σ 1 -σ 2 curve increase with the increase in the angle α, that is, the σ 2 effect on rock strength of these layered rocks strengthens gradually as the increasing angle α in the range of 0°-45°. Thus, the strength curves of layered rocks with α = 0° are almost straight lines as a result of the negligible effect of σ 2 . In terms of failure characteristics, these layered rocks are dominated by sliding failure Fig. 7 The Z-direction displacement diagram of the shear failure of the rock matrix; a α = 0°, β = 0°, and σ 2 = σ 3 = 0 MPa, b α = 45°, β = 60°, and σ 2 = σ 3 = 10 MPa Fig. 8 The Z-direction displacement diagram of the slipping failure along the weak plane; a α = 0°, β = 30°, and σ 2 = σ 3 = 0 MPa, b α = 45°, β = 30°, and σ 2 = σ 3 = 10 MPa 1 3 along the weak plane, as shown in Fig. 11. These mechanical characteristics are classified into the second type under the coupled effect of the σ 2 and the weak plane orientations.
In the third type, the mechanical characteristics of layered rocks with α > 45° and β = 15°-45° are complex, involved in the changes of the strength curve and failure mode. Due to the increase in angle α, the rock strength curves are characterized by two linear portions, as shown in Fig. 12. The slope of the previous linear portion changes from the minimum (the slope of strength curve at α = 60°) to the maximum (the slope of the critical line with σ 2 = σ 3 ) as angle α increases from 60° to 90°, whereas the slope of the next linear portion is mostly a constant for layered rocks (expect for these with α = 60°) at any stress σ 3 . In fact, the slope of strength curve decreases suddenly to a lower value as a result of the change of the failure mode. Taking the layered rock with β = 30° and α = 90° as an example, the failure modes corresponding to the previous and next portion of strength curve are sliding failure along the weak plane and shear failure of matrix, respectively, as presented in Fig. 13. Because of the lack of this change, there is only the previous portion for the layered rocks with α = 60°. It is evident that the domination of mechanical mechanisms (either rock matrix or weak planes) is determined by the coupled effect of the stress state and the weak plane orientations. Thus, the last type in mechanical characteristics of layered rocks can be regarded as a combination or transition of the former two types.

Discussion
It has been proved that the mechanical characteristics of layered rocks are determined by either the rock matrix or the weak plane based on the results of conventional triaxial tests. As mentioned above, the anisotropy of mechanical behaviors is closely related to the dip angle of weak planes (β). In the simulation results of this study, it is found that the anisotropy is also associated with the dip direction of weak planes (α) and the stress state (σ 2 and σ 3 ) under the true triaxial conditions. However, there are a very few existing true triaxial tests on layered rocks to identify this discovery. By comparing with the results of Chichibu green schist (Kwaśniewski and Mogi 2000), the results of the numerical simulation on layered rocks are discussed below.
In terms of the first-type property, the similar results are obtained in experiments and simulations. With the increase Fig. 9 The relationship curves between rock strength and intermediate principal stress in the first type; a α = 0°, β = 0°, b α = 45°, β = 0°, c α = 90°, β = 0°, d β = 90° 1 3 Fig. 10 The relationship curves between rock strength and intermediate principal stress at the angle α = 0° ~ 45° and β = 30°; a α = 0°, β = 30°, b α = 15°, β = 30°, c α = 30°, β = 30°, d α = 45°, β = 30°F ig. 11 The Z-direction displacement diagrams of layered rock with α = 0° and β = 30° at different stress states; a σ 3 = 0 MPa, σ 2 = 30 MPa, b σ 3 = 0 MPa, σ 2 = 60 MPa 1 3 Fig. 12 The curves of rock strength versus σ 2 at the angle α = 60° ~ 90°and β = 30°; a α = 60°, β = 30°, b α = 75°, β = 30°, c α = 90°, β = 30°F ig. 13 The Z-direction displacement diagrams of layered rock with α = 90°and β = 30° at different stress states; a σ 3 = 0 MPa, σ 2 = 30 MPa, b σ 3 = 0 MPa, σ 2 = 90 MPa 1 3 in σ 2 , the rock strength of the green schist with β = 90° increases gradually to an approximate constant; see Fig. 8 in the literature by Kwaśniewski and Mogi 2000. The failure patterns of these schists at different stress states all are the shear failure of the rock matrix, presented in the form of a schematic diagram as shown in Fig. 14a. Obviously, the mechanical mechanism of layered rocks entirely depends on rock matrix for both numerical simulation and experiment, whereas there are some differences between them. With the application of the Drucker-Prager constitutive model, the slope of the strength curve in the simulation remains constant, without the gradual decline as green schist. Thus, due to the overestimation of rock strength by the Drucker-Prager criterion (Colmenares and Zoback 2002), the use of it should be careful in engineering evaluation. Moreover, the failure surface is uncertain in the simulation as a result of the limitation of the constitutive model. Based on the whole experimental and simulated studies on layered rocks, it can be inferred that the mechanical characteristics of some layered rocks are related to the rock matrix, and similar to each other. The range of the weak plane orientations for these rocks (β = 0°, 60° ~ 90° and α = 0° ~ 90°) obtained from the simulation can provide reference for further study on layered rocks.
On the other hand, the mechanical characteristics of green schist with angle β = 30° are significantly affected by the weak plane. The σ 1 -σ 2 curves of these schists are in shapes of downward convex parabola; see Fig. 8 in the literature by Kwaśniewski and Mogi 2000. And the failure modes of these rocks exhibit the slipping failure along the weak plane, as shown in Fig. 14(b -c). When the angle α = 0° and β = 30°, the mechanical behaviors of green schist are roughly the same as those of layered rocks in the simulation due to the failure surface parallel to the σ 2 . At α = 45° and β = 30°, the failure mechanisms of layered rocks in simulation and experiment all are associated with the weak plane, but it is obvious that the rock strength obtained from simulation is overestimated at high σ 2 . Therefore, the failure criterion of the weak plane tends to overpredict its strength at some conditions. Besides, there is a prediction proposed by Mogi (2006) that a larger σ 2 may cause a change in the failure mode of the green schist with α = 90° and β = 30°, confirmed by the results of numerical simulation. But due to the limitations of theoretical models, these results still need to be identified by real tests.
Overall, the mechanical behaviors of layered rocks under the true triaxial stress state are complex, and the accuracy of numerical simulation results is limited as a result of the insufficient research on layered rocks. This study shows that the mechanical mechanism of layered rocks is strongly affected by the stress state and the weak plane orientations. As discussed above, there are some theoretical and experimental problems about layered rocks. For the demands of engineering construction and theory development, further studies on layered rocks are extremely essential.

Conclusions
By simulating the true triaxial compression tests on layered rocks with stress condition of σ 2 = σ 3 , the results of it are similar to those of real conventional triaxial tests, and it is found that the influence of the weak plane direction α on the rock mechanical characteristics is negligible. Based on it, the true triaxial tests with stress state of σ 2 > σ 3 were further simulated, and the following conclusions were obtained: (1) Under the true triaxial compression, the anisotropic mechanical characteristics of layered rocks are either controlled by the weak plane or dominated by the matrix. The failure mechanism of layered rocks is closely related to the stress state and orientations of weak planes, resulting in three different types of property. In each type, there are unique mechanical characteristics, including failure mode, strength curve and the σ 2 effect on rock strength, etc. 1 3 (2) When weak planes are at angle α = 0° ~ 90° and β = 0°, 60° ~ 90°, the mechanical characteristics of layered rocks are entirely associated with rock matrix, classified to the first type. Under the simulated condition, these specimens fail by the shear failure of the rock matrix, and the σ 1 -σ 2 curves at a certain σ 3 are the same as each other for these rocks. Thus, the σ 2 effect on rock strength only depends on the mechanical properties of the rock matrix. (3) In the second type, the mechanical behaviors of layered rocks with α ≤ 45° and β = 15° ~ 45° are related to the weak plane. The rock failure modes are all the sliding failure along the weak plane. And the σ 1 -σ 2 curves of these specimens are in the shape of a downward convex parabola, which is affected by the orientation of weak planes. As a result, the σ 2 effect on the strength of these layered rocks strengthens gradually as the angle α increases from 0° to 45°. (4) For the layered rocks with angle α > 45° and β = 15° ~ 45°, the mechanical characteristics are complex. Under the coupled effect of the σ 2 and the weak plane orientations, the changes of the strength curve and failure mode are observed in the results of simulation. The mechanical mechanism is determined by this coupled effect. Thus, the last type in mechanical characteristics of layered rocks can be regarded as a combination or transition of the former two types.

Conflicts of interest
The authors declare that they have no conflict of interest. Ethics approval The authors certify that this work is original and has not been published and will not be submitted elsewhere for publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.