Numerical simulation study on the mechanical behavior of shale core and reservoir rock mass

After the hydraulic fracturing of the shale reservoir, the stress field changes during the mining process. The stress sensitivity of the reservoir affects its development capabilities. Moreover, the stress sensitivity of shale reservoirs obtained by various researchers using different experimental methods is not uniform. This study analyzes the physical background differences between reservoir rock mass and experimental core. In this study, we establish the “support structure model” of the reservoir and show that the reservoir rock body has a structural overall effect, such that the supporting effect of the surrounding rock on the overlying strata is very strong. At the same time, we also establish the transverse isotropic numerical model of the shale reservoir on the basis of the real physical background of the reservoir and the experimental physical core model on the basis of the real physical background of the core test. This study shows that the equivalent mechanical behavior of the core test in the holder is pure core compression without the support and pull of the whole rock mass structure. In the real reservoir, the reservoir rock mass has significant structural effect. The actual compression change of real reservoir mining is less than the pure compression displacement of the experimental test core. The core stress sensitivity of the conventional experimental test does not directly describe the true stress‐strain effect of the macroscopic reservoir rock mass as a structural whole. A laboratory stress‐sensitive test apparatus and method that is equivalent to the physical background of a real reservoir needs to be improved or designed.

The permeability obtained by the method of constant pore fluid pressure and changing confining pressure is sensitive to the confining pressure and is exponentially related to the effective stress. [2][3][4][5] Jiao et al 6 applied the double effective stress theory to evaluate the core stress sensitivity using the method of constant pore fluid pressure and changing confining pressure. The result showed that there is no strong stress sensitivity in the ultralow-permeability reservoir. Gao et al 7 used the method of constant confining pressure and changing pore fluid pressure to study the permeability sensitivity of ultralow-permeability reservoirs. The result showed that permeability is less sensitive to pore fluid pressure. Zhu et al 8 conducted a permeability stress sensitivity test by using both the methods simultaneously. The result indicated that the smaller the effective stress (Boit) coefficient obtained by the experiment is, the stronger is the exponential relationship between the permeability and the effective stress, and the stronger is the stress sensitivity of the shale. Shen et al 9 found that the permeability sensitivity of shale fractures under stress was very strong, and the exponential relationship better described the pressure dependency of the permeability for the tested shale samples. For high-pressure deep shale gas reservoirs, which ware usually highly stress-sensitive, the impact of strong stress sensitivity should be fully considered. 10 These studies on stress sensitivity are abundant, but the results are inconsistent. Researchers have also tried to find a reasonable explanation from the experimental process and principles, but there is still no consensus.
In terms of pore compression stress sensitivity, Luo et al 11 found that the low-permeability and ultralow-permeability reservoirs are basically comprised of small pores, and hence, the capillary effect is strong, so there is a large threshold pressure gradient. Therefore, it has stronger stress sensitivity compared with the medium-to-high permeability reservoirs. Luo et al 12 combined rock mineral analysis data and mercury intrusion test to study the microporous structure of reservoirs. It was believed that skeleton particle content, pore throat radius, and clay mineral content affect the stress sensitivity of reservoirs. Therefore, the effects of stress sensitivity should be considered during the development of ultralow-porosity and -permeability reservoirs. Tian et al 13 showed, while analyzing the relationship between shale stress sensitivity and shale productivity, that the stress sensitivity curve of shale with fracture was L-shaped. Pseudoplastic deformation was predominant in the first stage, whereas the elastic compression deformation of the rock skeleton was mainly exhibited in the second stage. As a consequence, the stress sensitivity of fractured reservoirs was relatively weak. On the basis of the research of reservoir rock compressibility, Li 14 believed that the compression of the skeleton volume led to the compression of rock pore volume. The low-permeability compact rock had a low pore compression coefficient, so the stress sensitivity of low-permeability reservoirs was weak. Furthermore, Wang et al 15 found that the stress sensitivity coefficient of permeability was positively correlated with the aspect ratio of pores and the elastic modulus of matrix. With the increase of organic matter content, the stress sensitivity coefficient of permeability decreased under the combined effect of high quartz content and low aspect ratio organic pores. These pore compression stress-sensitive studies have mainly two different results. One is that shale reservoirs mainly develop a large number of micropores and nanopores, resulting in a large threshold pressure gradient and strong stress sensitivity. The other is that shale reservoirs have obvious high-pressure characteristics, and ultradense and low compressibility factors will lead to weak stress sensitivity.
In short, the current research status of shale stress sensitivity can be summarized as shown in Table 1.
Through the above survey, we can find the stress sensitivity of shale reservoirs directly reflected by the stress sensitivity of shale core permeability and pore compression is not uniform. Most researches mainly concentrated on experimental methods and theoretical analysis, but whether the experiment itself is similar to real reservoirs has not attracted much attention. The most basic and important step for simplifying an engineering problem to a physical model is similarity transformation and mechanical modeling. In this paper, the difference in physical background and mechanical behavior between shale reservoir rock mass and experimental core is studied. The mismatch between the experimental model itself and the real reservoir was found. It reveals the root cause of the inconsistent results of core stress sensitivity.

CHARACTERISTICS AND DEVELOPMENT BACKGROUND OF SHALE RESERVOIR
Rocks are aggregates of minerals that make up the crust, which is different from the rock mass. Rock mass refers to a composite geological body formed during the geological history process, with certain rock components, specific structural types, and certain environmental conditions. The two basic properties of rock mass are as follows: material composition and structure. 17 The constituent elements of the structure have two basic units: the structural plane and the structural body. The structural plane refers to a two-dimensional planar geological interface with a certain direction, a large extension, and a small thickness. It includes material dissimilar surfaces and discontinuous surfaces, such as faults, joints, bedding, schistosity, weak interlayer, unconformity surface, and so on. The structural body refers to rock blocks that are cut by structural planes. Table 2 shows the mechanism of rock mass deformation. 18 The layered rock mass structure is layered or thin-layered. The shale rock structure type belongs to the thin-layered structure. The geological background is a thin-layered rock mass. Under tectonic action, it exhibits relatively strong folds and interlayer movement. The genetic type of the structural plane is the sedimentary structural plane of the original structural plane. The appearance is generally consistent with the occurrence of the rock layer, which is the interlayer structure surface. The distribution of such structural planes in marine rocks is stable. The deposition compaction is strong, and the structure is stable. Therefore, the overall deformation modulus and bearing capacity of this type of rock mass are relatively high. Rock mass deformation is controlled by rock layer combination and structural plane.
The stratigraphic rocks, a research unit, belong to a porous solid medium. Its structure consists of three parts: rock matrix (skeleton), pores (fissures), and fluid, and the structural properties of these three components are related to the macromechanical properties of the rock manifestation. With the continuous exploitation of oil and gas, the pore pressure of oil and gas reservoirs continues to decrease. The effective stress of the rock skeleton of the reservoir changes, which changes the mechanical properties of the rock and deforms the rock skeleton. This leads to changes in the pore volume in the rock, which in turn leads to changes in the permeability of the pore fluid and the state of oil, gas, water, and water (saturation) of the rock. When the rock is not under pressure, the pores in the rock coexist with the throat. When stressed, the throat in the rock is closed first, but the pores are not closed. As the pressure increases, there are fewer and fewer unclosed throats, and most are throats that are not easy to close, resulting in a reduction in the amount of rock compression after compression. The cores belong to drilled rock samples, and current stress-sensitive studies are basically based on the compressive nature of the rock, ignoring the structural integrity of the rock. In fact, the macroscopic real stratigraphy is comparable to continuous solid materials in terms of strength, deformation, and external loading relationships due to the joint action of rock structure and material composition. During the mining process after the reservoir fracture, the overlying strata will be supported by a certain degree of surrounding rock due to the structure of the rock mass. It transforms and decomposes a certain degree of stress changes in the transformation area and slows the compression of the roaring channels and pores in the rocks in the transformation area. By analyzing the mechanical mechanism of real shale reservoir mining and experimental physical simulation testing, as shown in Figures 1 and 2, it is pointed out that the physical background of the experimental physical simulation testing is quite different from the real reservoir. The stress sensitivity of the experimental physical simulation test cannot directly describe the real effect of the reservoir rock mass as a whole. Considering the macroscopic continuity and structural integrity of the reservoir rock mass, in this study, we analyzed the supporting degree of the surrounding rock to the overlying strata. The simulation of rock sample compression and the core held by the holder were carried out at the full-diameter core scale, and the equivalent mechanical behavior of the core held by the holder was revealed. A numerical model of large-scale reservoir was also developed. Moreover, the stress and strain values of the reservoir elements in the numerical model were compared with those in the core held by the holder model at the core scale, indicating the difference between the core held by the holder and the real reservoir stress-strain characteristics.

| Elastic constitutive model
It was assumed that the deformation of reservoirs was in the elastic range. For the natural foundation formed by sediments such as shale reservoirs, the vertical physicochemical properties of rocks were quite different, and the horizontal nature difference was relatively small. Therefore, a transversely isotropic elastic-plastic constitutive model was introduced to describe the shale reservoir state. 19 The basic equations describe the axisymmetric space problem of transversely isotropic elastomers in Cartesian coordinates, including the equilibrium equations, geometric equations, and the physical equations described by Hooke's law (constitutive relations). 20 The stress component represents the equilibrium differential equations, as written in the following equation: The geometric equation is The physical equation can be written as where E v is the elastic modulus perpendicular to the isotropic plane, E h is the elastic modulus parallel to the isotropic plane, v vh is Poisson's ratio of the horizontal strain caused by the application of vertical strain, v hh is Poisson's ratio in isotropic plane, G vh is the shear modulus perpendicular to the isotropic plane, Δ x , Δ y , and Δ z are the increment of strain, Δ x , Δ y , and Δ z are the increment of normal stress, Δ yz , Δ zx , and Δ xy are the increment of shear strain, and Δ yz , Δ xz , and Δ xy are the increment of shear stress.

| Plastic constitutive model
The plastic constitutive model is composed of three parts: yield criterion, flow rule, and hardening criterion. The bedding of shale is mainly horizontal. The bedding plane is isotropic. The transverse isotropic plastic yield function of rock derived by Wang et al 21 is shown in Equation (4).
where g ( ) is the yield function, F and H are model parameters, is the dilatancy angle, D is the cohesive force, I 1 is the first invariant of the stress tensor, and 11 , 22 , and 33 are the first, second, and third principal stresses, respectively.
The flow rule is shown in Equation (5).
The potential function f of the geotechnical materials is expressed as shown in Equation (6).
where d is the plasticity factor; f is potential function; = sin √ 3 √ 3+sin 2 , and is internal friction angle.
In this paper, referring to the study of Wang et al, 21 the strain hardening criterion is used to describe the properties of geotechnical materials under complex loading conditions.

| Rubber constitutive model in core holder
The elastic constitutive models describing rubber materials include a phenomenological model based on the strain energy function of the continuous medium and a statistical model considering the statistical thermodynamic molecular chain network. It is well known that the rubber material is a super-elastic material, and even if it is quite deformed, it can be treated similarly as an isotropically uniformly deformed elastic material. The deformation process of the elastic material is reversible, and the elastic constitutive equation is mainly expressed by the strain energy function. In this study, we neglect the microstructure and molecular nature of the elastomer, and the Mooney-Rivlin phenomenological model, which considered the elastomer as a continuous medium, was used to describe the mechanical behavior of rubber. 22 Mooney-Rivlin model is a commonly used model, which can simulate the mechanical behavior of almost all rubber materials. The strain energy density function model was given by Equation (7).
For incompressible materials, the typical binomial third-order expansion is as Equation (8).
where N, C ij , and d k are material constants determined experimentally.

MODEL
In this section, we studied the overlying strata to demonstrate the effect of surrounding rock caused by the integrality of rock mass structure on the overlying strata. The effect of the overall rock mass structure on the supporting and pulling degree of the overlying strata is analyzed by comparing an extremely hypothetical support structure model (the box area of different transformation sizes is killed by the life and death element technology) with the reservoir solid compression model.

| Equal-overlying strata thickness support structure model with different spans
A solid model with length, width, and height of 200 m × 60 m × 40 m is established using the finite element software ANSYS. The vertical wellbore direction in the horizontal plane is in the x-axis, the vertical reservoir direction is in the y-axis, and the direction of the wellbore is in the z-axis. The origin of the coordinates is at the center of the model. Furthermore, a series of equal-overlying strata thickness support structure models with different spans were created (in the center of the solid model, the elements in the boxes of size 100 m × 40 m × 10 m, 50 m × 40 m × 10 m, 20 m × 40 m × 10 m, and 10 m × 40 m × 10 m were killed). The effect of surrounding rock with different spans on overlying strata is studied by comparing the results of equal-overlying strata thickness support structure models with those of different spans and solid model.
The solid element SOLID64 is used to simulate rock, and the shale transverse isotropic material parameters are introduced, which can be described by the parameters shown in Tables 3 and 4.
The grid size in this model is controlled in element size, and the total number of grids is 480 000 of the solid model. A fixed constraint was applied to the bottom and sides of the model, and a uniform pressure of 10 MPa was applied to the top. It should be noted that the prototype of this model is stratum. The bottom and sides of the stratum are approximately fixed. Therefore, when exploring the influence of a certain part of the model's stress changes on the whole, we applied fixed constraints on the bottom and sides of the model under the condition of satisfying Saint-Venant's principle. In addition, there is an interaction between the overlying strata and the lower reservoir when the reservoir is mined. There will be a pressure difference between the overlying strata and the lower reservoir. With the increase of mining time, this pressure difference can reach more than 10 6 Pa. Here, we choose to study the mechanical behavior of various equivalent reservoir models when the pressure difference is 10 6 Pa. Figure 3 shows the stress and displacement nephograms of the solid model and the equal-overlying strata thickness support structure model with different spans. Figure 4 shows the stress and displacement distribution curves of the long axis corresponding to the bottom plane of the overlying strata.
The stress distribution curve indicates that the stress of the surrounding rock increases in the solid support area of the overlying strata. The stress at the bottom of the overlying strata decreases initially from the plane center to the surrounding rock support interface and then it increases. It shows that the pressure transmitted by the overlying strata is partially decomposed and assumed by the surrounding rock of the reservoir when it reaches the lower mining reservoir. In this example, when the span is 10 m, the stress in the bottom region of the overlying strata is less than that of the solid model, indicating that the supporting effect of the small-span surrounding rock on the overlying strata is excellent. The displacement distribution curve reflects that the surrounding rock span has a great influence on the deformation of the overlying strata supported by the surrounding rock, and the maximum displacement does not change uniformly with the span. In this case, for the convenience of comparison, the degree of support of the rock mass is defined as Equation (9). where sd is the degree of support (%), ∆ max d the vertical maximum relative displacement (m), and s the supported span (m). The support degrees of supported spans of 100, 50, 20, and 10 m are 88.715%, 86.226%, 90.895%, and 93.405%, respectively. The result also indicates that the support effect of the small-span surrounding rock on the overlying strata is extremely significant, which must also be considered. The degree of support changes a little with increasing span to a larger range, and it is in the hardening stage which has a risk of damage. When the reservoir is mined, the pore pressure drops and the rock mass are compressed, with varying degrees of solid bedrock and surrounding rock support in both the mining zone and the overlying strata. The pressure transmitted from the overlying strata to the reservoir will be redistributed within the rock mass according to the structure support and bedrock deformation and then it reaches an equilibrium state.

| Equal span support structure model with different overlying strata thickness
On the basis of the 4.1.1 solid model, a series of equal span support structure models with different overlying strata thickness of 5, 10, 15, and 20 m (the elements in the boxes of size 10 m × 40 m × 10 m in the center of the solid model were killed) were established. They are used to study the supporting effect of surrounding rock on overlying strata of different thickness. Figure 5 shows the stress and displacement nephograms of the equal span support structure model with different overlying strata thicknesses obtained by finite element solution. Figure 6 shows the stress and displacement distribution curves of the central long axis corresponding to the bottom plane of the overlying strata with different thicknesses.
We infer from Figure 6 that the greater the thickness of the overlying strata on the top of the support structure is, the smaller is the stress value at the bottom of the overlying strata, and the stronger is the surrounding rock supports. The supporting extents of the overlying strata of 5, 10, 15, and 20 m are 78.16%, 93.36%, 96.35%, and 97.5%, respectively. It shows that the thicker the overlying rock layer is, the stronger is the ability to decompose and transform the propagation pressure. That is to say, the thicker the overlying strata are, the less the influence of the pressure transmitted by the overlying strata on the stress of the reservoir mining area is. The results also indicate that the overall support effect of the rock mass structure is significant and cannot be ignored.

| Core test numerical model
The simulation of rock sample compression and the core held by the holder were carried out at the full-diameter core scale, and the equivalent mechanical behavior of the core held by the holder was revealed.
A full-diameter core compression simulation model was established with a diameter of 0.1 m and a length of 0.4 m. The core simulation element and parameters were the same as the transverse shale isotropic solid model, as seen in Section 4.1. This model is controlled by a mesh size of 0.005 m. According to our prototype of physical simulation experiment, 25 a fixed support is applied to both ends, and an effective uniform stress of 10 6 Pa is applied to the side.
In addition, a simulation model of the full-diameter core held by the holder was developed. The core in the holder is held by the steel blocks at both ends, and the confining pressure is loaded by the rubber sleeve. 25 The core diameter is 0.1 m, and the length is 0.15 m, 26 and the simulation element and parameters are the same as the shale transverse isotropic solid model, as seen in Section 4.1. In the core holder part, the diameter of the steel block clamped at both ends is 0.1 m, the length is 0.125 m, the inner diameter of the rubber sleeve is 0.1 m, the outer diameter is 0.11 m, and the length is 0.4 m.  Both rubber and steel blocks are simulated using SOLID185 elements that can comprehensively simulate superelasticity and small deformation. Tables 5 and 6 present, respectively, the material parameters of rubber and steel blocks. Also, the contact pairs between the contact faces of different materials were built by the unique contact guide set in ANSYS, which aims to transmit the contact pressure (positive pressure and friction). Specifically, the rigidity of the core is significantly greater than that of the rubber, so the contact between the core and the rubber is the contact between the rigid body and the flexible body. The two are in symmetrical contact. The core surface is the target surface, and the rubber surface is the contact surface. The friction coefficient is set to 0.2, and the normal contact stiffness is taken as 3 (the contact stiffness can be updated automatically). The rigidity of the core and the steel block is similar. The two are in symmetrical contact. The friction coefficient is set to 0.55, and the normal contact stiffness is taken as 0.1 (the contact stiffness can be updated automatically). Other options select the default settings for the software. The boundary conditions are set to be the same as the full-diameter core compression simulation model. Figures 7 and 8 present, respectively, the displacement nephograms of the full-diameter core compression model and the full-diameter core held by the holder model.   Figure 9 compares the distribution curves of the stress and displacement on the cylinder generatrix corresponding to different radii of the two models.
We understand from Figure 9 that in the full-diameter core compression model, the stress distribution inside the core is basically the same, and the displacement decreases uniformly. It also shows that the core is purely evenly compressed. In the full-diameter core held by the holder model, Figure 10 reflects the friction hindering effect of the core and the clamping steel blocks at both ends, but it is far from enough to simulate the displacement limit of the core shear resistance. The rubber casing has a uniform distribution of confining pressure. The stress and displacement curves of the core in the core held by the holder model are basically consistent with the core compression model. The results indicate that the equivalent physical background of the core held by the holder model is simple core compression, without the overall structural support and pulling of the rock mass.

| Comparative analysis of core and reservoir models
A large-scale shale reservoir model is established with the coordinate origin as the center, with the dimension of the reservoir model is 400 m length × 240 m width × 40 m height. The length, width, and height of the reservoir mining zone are 100 m × 60 m × 10 m, located in the center of the reservoir model. The coordinate system, element type, material parameters, and meshing used in the model are exactly the same as the shale transverse isotropic solid model, as seen in Section 4.1. The bottom and sides are fully constrained, and an effective uniform stress of 10 6 Pa is applied to the model. Figures 10 and 11 present, respectively, the stress nephograms of the core held by the holder model and the reservoir model obtained by the finite element method. Figure 12 shows the stress changes for the reservoir model, and Figure 13 compares the stress and displacement distribution curves.
It can be seen from Figure 13 that in the reservoir model, the stress and displacement of the bottom area of the overlying strata fluctuate, and the supporting effect on the mining area is relatively balanced. The results indicate that the structural integrity of the reservoir rock mass plays a significant role. At the same time, the strength and deformation characteristics of the rock mass are comparable to the continuous solid material. This is also consistent with the generally deep burial, high-pressure, ultradense, and high compressive shear strength characteristics of shale reservoirs.
During the mining process, after the reservoir fractures, the pressure transferred from the overlying strata to the reservoir mining zone will be redistributed within the rock mass, as the support degree of the reservoir changes until a new balance is reached. In the process of dynamic equilibrium, the stress and displacement transferred from the overlying strata to the mining area will be transformed and decomposed to a certain extent by the supporting effect of surrounding rock and bedrock. Since the forces are mutual, the compression of the channels and pores in the rock in the mining area will also be weakened. However, the experimental core test cannot simulate the structural function of the rock mass. The pure