Cross-scale bedding shale elastic mechanical property calculation method

Abstract

The invention discloses a cross-scale bedding shale elastic mechanical property calculation method. According to the method, a geometric model is established on a microscopic scale by using a Thiessen polygon method, the geometric model contains various minerals and bedding planes, then tensile and compressive loads are applied to solve mechanical characteristics, and the mechanical characteristics are compared with indoor experimental results for verification; on the macroscopic scale, the mechanical property of the bedding shale is predicted through periodic boundary conditions and a mechanical equivalent theory, cross-scale solving is achieved, numerical simulation of a uniaxial compression test is conducted on the bedding shale on the basis of a finite element method, and macroscopic elastic mechanical parameters of the bedding shale are obtained. According to the method, stress-strain data of the bedding-containing shale can be determined, mechanical parameters such as the elastic modulus and the Poisson's ratio are calculated on the basis, and a reliable basis is provided for shale fracturing transformation.

Description

Technical Field

[0001] The present invention relates to the research field of oilfield rock mechanical properties, and specifically to a method for calculating the elastic mechanical properties of cross-scale layered shale based on the combination of microscopic and macroscopic finite element simulations. Background Art

[0002] With the deepening of oil and gas resource exploration and development, conventional oil and gas resources alone can no longer meet my country's growing energy demand. my country's shale oil and gas resources are rich in reserves and have huge development potential. The effective development of shale oil and gas reservoirs will provide an important way to alleviate my country's energy shortage. Shale is a layered sedimentary rock with significant heterogeneity, ultra-low porosity and low permeability. Horizontal well multi-cluster fracturing technology is an effective means to develop shale oil and gas. The key to the success of this technology lies in whether effective extended cracks can be formed in shale reservoirs after hydraulic fracturing. However, the actual expansion of hydraulic fractures in rocks is affected by formation conditions and the mechanical properties of reservoir rocks, and is highly complex. Therefore, in-depth research on the mechanical properties of shale is of great significance for shale oil and gas development.

[0003] Shale can be regarded as a typical heterogeneous natural material, composed of multiple minerals, and exhibits significant anisotropy and multi-scale characteristics. Bedding is a common structural feature in shale, and bedding structure can affect the fracturing effect. It is difficult for hydraulic fractures to penetrate bedding vertically, which hinders the communication between adjacent shale oil reservoirs and makes it difficult to hydraulically fracture horizontal wells. At present, commonly used rock mechanical parameters, such as tensile strength, Young's modulus and Poisson's ratio, are mainly obtained through indoor experiments and logging data. For example, triaxial stress tests and Brazilian splitting tests can measure stress-strain curves, and then deduce mechanical parameters such as elastic modulus, Poisson's ratio and compressive strength. The results obtained by this method have high reliability and accuracy. However, in the process of obtaining shale mechanical parameters, there are the following difficulties: first, shale sampling is difficult, and indoor experiments are difficult to fully simulate formation conditions; second, the high brittleness of shale makes it easy to break during the test, which increases the complexity of indoor experiments. Therefore, there are certain limitations in obtaining shale mechanical parameters using indoor experiments. Compared with indoor experiments, numerical simulation experiments are more efficient, controllable, flexible and safer. They can study the mechanical properties of shale more comprehensively and in-depth, and provide more effective theoretical guidance and technical support for shale oil and gas extraction.

[0004] The composition, layout, location and mechanical properties of different minerals in shale directly affect the macroscopic mechanical properties of shale. If the elastic mechanical properties of layered shale can be solved and verified at the mesoscale, and the elastic mechanical properties can be solved at the macroscale through mechanical equivalence theory and the imposition of periodic boundary conditions, this will comprehensively characterize the mechanical properties of layered shale from multiple angles and facilitate the study of the influence of single factors on the mechanical properties of rock. Summary of the invention

[0005] In view of the above problems and research status, the present invention first establishes a geometric model at a microscopic scale using the Thiessen polygon method, which contains a variety of minerals and bedding planes, and then applies tensile and compressive loads to solve the mechanical characteristics, and compares and verifies with indoor experimental results; then, at a macroscopic scale, the mechanical properties of the bedding shale are predicted through periodic boundary conditions and mechanical equivalent theory, and cross-scale solutions are achieved. Based on the finite element method, a uniaxial compression test of the bedding shale is numerically simulated to obtain the macroscopic elastic mechanical parameters of the bedding shale.

[0006] The technical solution provided by the present invention is: a cross-scale method for calculating the elastic mechanical properties of bedding shale, comprising the following steps: 1) Establishment of Thiessen polygon model: The cell shape of the three-dimensional Thiessen polygon Voronoi shows good similarity with the shape of mineral particles. Therefore, the multi-cell of the Voronoi diagram is used to represent the structure of the internal particles of the rock. The Thiessen polygon is a set of continuous polygons composed of perpendicular bisectors connecting two adjacent point segments. The distance from any point in a Thiessen polygon to the polygon is less than the distance to the control points of other polygons. According to the Thiessen polygon coordinates, the Python script is used in the abaqus software to generate the Thiessen polygon model; 2) Assign different mineral grain volumes to the polygonal model according to the mineral ratio of the bedding shale, and calculate the grain volume of each Thiessen polygon; 3) Use Python script to randomly place mineral particles at any position of the polygonal model according to the calculated volume, and assign different attributes to the polygonal model containing mixed mineral particles; 4) The polygonal model is divided into tetrahedral meshes as a verification model, with a global mesh size of 3 and a mesh attribute of C3D4. By adding reference points and applying stress surface coupling, loads are applied to the rock, i.e., tension / compression, to create fixed boundary conditions. One end of the model is fixed, i.e., the degree of freedom of its end face in the X direction is set to 0, and the degrees of freedom of the three corner points in the XYZ direction, XY direction, and XZ direction are also set to 0. The other end is compressed, i.e., a displacement of 1 mm is added to its end face along the X direction, and the degrees of freedom of the three corner points in the YZ direction, Y direction, and Z direction are also set to 0; 5) Select the static general mode to compress the polygonal model at the microscopic scale, and set the total compression time to 1s; 6) Create a new geometric model of the same size as the polygonal model, i.e., 50×50×100. Set four reference points (0,0,20), (50,0,20), (0,05,80), and (50,50,80) respectively. Perform three-point segmentation on the model and insert cohesive units as bedding planes on the segmentation planes. Specify the unit type as COH3D8. The bedding plane forms an inclination angle with the geometric model, which can be determined according to the bedding angle of the actual shale. 7) Create a repeated mesh pattern for the geometric model in a specific area or direction to simulate periodic behavior; first, perform tetrahedral meshing on the geometric model, then offset and copy the mesh faces of the geometric model so that the meshes of the opposite faces are the same, then divide the mesh volume to achieve the periodicity of the geometric model mesh, define the global mesh size as 3, and specify the model mesh attribute as C3D4; 8) Map the periodic grid divided in the geometric model to the polygonal model. According to whether the center coordinates of the unit are in the corresponding geometric area, the grid is mapped by coordinate comparison. The mapped model is a representative volume element RVE model. RVE can provide sufficient accuracy to represent the larger scale of the material. 9) In order to eliminate the boundary effect caused by the scale limitation of the simulation system and make the stress and displacement at the boundary continuous, periodic boundary conditions are created for the geometric model; Python script is used to add periodic boundary conditions: First, six boundary surface sets are established and read for the polygonal model, and the node data on the six surfaces are saved; second, the node number and coordinates on each surface are read; third, the nearest points on each pair of relative surfaces are calculated, and the node pair numbers are saved; fourth, the eight diagonal points are numbered and divided into two groups, namely the lower left, lower right, upper left, and upper right; fifth, equation constraints are established for the node pairs on the left and right surfaces, the upper and lower surfaces, and the front and back surfaces, to constrain the degrees of freedom in the three directions of X, Y, and Z; finally, equation constraints are established for the eight diagonal points; 10) Perform compression and tensile tests on polygonal models at the microscopic scale to solve elastic mechanical properties such as stress-strain curve, compressive strength, elastic modulus and Poisson's ratio: 11) On a macro scale, geometric models of bedding in stratified shales with bedding angles of 0°, 30°, and 60° can be established respectively, and then solutions and predictions can be made based on the model building and calculation processes in steps 6) to 10).

[0007] In the above step 2), the volume is calculated using the area and elevation difference of Thiessen polygons. First, the area of ​​each Thiessen polygon is calculated using GIS software. Then, the elevation difference between the highest point and the lowest point in each Thiessen polygon is calculated, such as formula 1-formula 7 to calculate the volume of each Thiessen polygon. ; Where l, w, h are the length, width and height of the model respectively; V and Vquartz, Vcalcite and Villite are the total volume of the polygonal model and the volume of quartz, calcite and illite mineral particles in the model respectively; wquartz is the volume fraction of quartz; v is the volume of any particle in the polygonal model; s is the area of ​​any particle in the polygonal model; Δh is the elevation difference between the highest point and the lowest point in any particle in the polygonal model.

[0008] In the above step 9): Surface node constraint equations: On the opposite surface perpendicular to the X axis: ; On the opposite surface perpendicular to the Y axis: ; On the opposite surface perpendicular to the Z axis: ; Constraint equations between edge nodes: Three sets of linear constraint equations between the four edges parallel to the Z axis: ; Constraint equations between corner nodes: ; In the formula, ε 0 x , ε 0 y , ε 0 z , γ 0 xy , γ 0 xz , γ 0 yz are the strains in the X, Y, Z, XY, XZ, and YZ directions respectively; x =W x , y =W y , z = hThe three planes are called principal planes, and the planes parallel to the principal planes are called subordinate planes. x , W y , h They are the length, width and height of the cube polygonal model respectively; u , v , w They are x , y , z The displacement in the direction, the subscripts represent different corner nodes (D, E, F, G), edge nodes (EA, HD, FB, GC), surface nodes ( x =0, x =W x , y =0, y =W y , z =0, z = h ) position, that is u ∣ x=w for x direction x =W x The displacement at the surface node position, v EA for y The displacement at the edge node in the EA direction, w F for z Direction F displacement at the angle node.

[0009] In the above step 10), firstly, the reaction forces received by all the collection points on the tensile surface are summed up using formula (17), and the stress-strain curve of the entire representative volume element is plotted using formulas (18) and (19). Then, according to the displacement deformation of the model in the X direction and the Y direction (or the Z direction), the mechanical parameters of the bedding shale are obtained using formulas (20-23), and compared with the results of the indoor rock compression test, the compression cloud map and stress-strain curve of the bedding shale and black shale.

[0010] The sum of the reaction forces at all points on the extruded surface: ; Stress and strain in the X direction: ; Strain in Y and Z directions: ; Poisson's ratio: ; Modulus of elasticity: ; Where: rfi is the first i The reaction force of the point, RF x is the sum of all reaction forces on the plane; A x is the area of ​​the plane perpendicular to the X direction; u x , u y , u z are the displacements in the X, Y, and Z directions respectively, l x , l y , l z They are the maximum dimensions in the X, Y, and Z directions respectively; v is Poisson's ratio, σ x is the stress in the X direction, ε x , ε y , ε z are strains in the X, Y, and Z directions respectively; E is the elastic modulus, σ x1 , σ x2 , ε x1 , ε x2 are the stress and strain at any two points before the highest point of the stress-strain curve in the X direction.

[0011] The beneficial effects of the present invention are as follows: the present invention establishes a Thiessen polygon model of layered shale, adopts a stretching and compression method, enables the model to perform operations in a target loading environment, and can clearly describe the process of stretching the shale containing layers until the structure is destroyed by considering the mineral type and content through the output cloud map, clarify the stress and strain data of the layered shale, and calculate the elastic modulus, Poisson's ratio and other mechanical parameters on this basis, providing a reliable basis for shale fracturing transformation. The specific advantages are as follows: a) Considering the bedding shale as a heterogeneous composite material composed of various mineral particles and bedding, the elastic mechanical properties of the bedding shale are calculated by combining the random arrangement of mineral particles in the stratum and the effect of bedding on the shale, taking into account the heterogeneity of the shale. The results are compared and verified with the indoor experimental results, and the elastic mechanical properties of the shale with beddings of different dip angles can be predicted.

[0012] b) As a novel mechanical characterization and solution method, this method establishes representative volume elements based on Thiessen polygons and creates periodic boundary conditions to calculate the elastic mechanical properties of bedding shale. Unlike the previously disclosed patent "A method for solving the elastic mechanical properties of shale reservoir rocks from micro to macro", this patent uses a planar Thiessen polygon model to establish representative volume elements of various minerals randomly distributed and embedded in bedding surfaces, which is more in line with the actual situation of bedding shale and provides a more accurate solution.

[0013] c) Compared with indoor experimental methods, this method can simulate actual formation conditions and is suitable for solving the elastic mechanical parameters of layered shales with different mineral contents. It can provide mechanical parameters for a variety of fracturing fracture extension simulation schemes under various formation conditions, and can reveal the influence mechanism of various factors such as rock physical properties, mineral proportions, arrangement, clay mineral hydration, etc. on the elastic mechanical properties of shale reservoir rocks. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Attached Figure 1 : Thiessen polygon graph; Attached Figure 2 : Generate Thiessen polygon flow chart; Attached Figure 3 :Mineral proportions of bedding shales in northeastern Sichuan; Attached Figure 4 : attribute model diagram; Attached Figure 5 : Displacement load boundary diagram; Attached Figure 6 : Bedding plane diagram; Attached Figure 7 : Periodic grid generation flow chart; Attached Figure 8 : periodic grid graph; Attached Fig. 9 : RVE schematic diagram; Attached Fig.10 : Mapping grid map; Attached Fig.11 : Equation boundary condition diagram; Attached Fig.12 : represents the boundary condition diagram of the volume element model; Attached Fig.13 : Layered shale compression cloud map; Attached Fig.14 : Black shale compression cloud map; Attached Fig.15 : Stress-strain curves of layered shale and black shale; Attached Fig.16 : Technical roadmap for cross-scale calculation methods of elastic mechanical properties of layered shales. DETAILED DESCRIPTION

[0015] The present method is not limited to the following examples. In the detailed description of the present method below, some specific details are described in detail. However, for the part that is not described in detail, those skilled in the art can also fully understand.

[0016] A cross-scale calculation method for elastic mechanical properties of bedding shale is proposed. The specific process is as follows: 1. Generation of Thiessen polygon model: Thiessen polygon, also known as Voronoi diagram, is a set of continuous polygons composed of perpendicular bisectors connecting two adjacent line segments. The distance from any point in a Thiessen polygon to the polygon is less than the distance to the control points of other polygons. According to the Thiessen polygon coordinates in Appendix 1, a Python script (code in Table 2) is used in abaqus software to generate a Thiessen polygon model (polygon model) (see Appendix Figure 1 ): First define the point set in Appendix 1, then use the Delaunay function in the Python library to build a Delaunay triangulation, find the boundary of the Thiessen polygon, and finally draw the Thiessen polygon. See the flowchart in the attached Figure 2 , Figure 2 In: The Delaunay triangulation is a collection of connected but non-overlapping triangles whose circumscribed circles do not contain any other points in the region.

[0017] Appendix 1 Thiessen polygon coordinates RVE x y z Starting 0 0 0 end 50 50 100

[0018] 2. According to the mineral proportions of the bedding shale in northeastern Sichuan (see Appendix Figure 3 ), the volume of different mineral particles is assigned to the polygonal model. The mineral content of the bedding shale and black rock is shown in Appendix 3 (black shale is used for comparison with the bedding shale). The calculation process of the particle volume of each Thiessen polygon is as follows: The volume is calculated using the area and elevation difference of Thiessen polygons. First, use GIS software to calculate the area of ​​each Thiessen polygon, and then calculate the elevation difference between the highest point and the lowest point in each Thiessen polygon, as shown in formula (1-7) to calculate the volume of each Thiessen polygon.

[0019] ; Where l, w, h are the length, width and height of the model respectively; V and Vquartz, Vcalcite and Villite are the total volume of the polygonal model and the volume of quartz, calcite and illite mineral particles in the model respectively; wquartz is the volume fraction of quartz; v is the volume of any particle in the polygonal model; s is the area of ​​any particle in the polygonal model; Δh is the elevation difference between the highest point and the lowest point in any particle in the polygonal model.

[0020] 3. Use Python script (part of the code is shown in Appendix 5) to randomly place the mineral particles at any position of the polygonal model according to the calculated volume. The polygonal model containing mixed mineral particles is the attribute model, and different attributes are assigned to it. The attribute table is shown in Appendix 4. (The effect diagram is shown in Appendix Figure 4 ).

[0021] Appendix 3 Shale Bedding Rock Mineral Content Table Lithology Quartz content / % Calcite content / % Illite content / % Stratified shale 33 22 45 Black shale 35 5 60 Appendix 4 Mineral particle properties mineral Young's modulus / MPa Poisson's ratio Compressive strength / MPa quartz 85330 0.31 250 Calcite 90153 0.073 208 Illite 7900 0.313 / Schedule 5 Allocation of mineral particle codes

[0022] 4. Divide the polygonal model into tetrahedral meshes as the verification model, with a global mesh size of 3 and mesh attributes of C3D4. Add reference points and apply stress surface coupling to apply loads (tension / compression) to the rock, create fixed boundary conditions, fix one end of the model, that is, set the X-direction degree of freedom of its end face to 0, and set the degrees of freedom of the three corner points in the XYZ direction, XY direction, and XZ direction to 0 respectively. Compress the other end, that is, add a displacement of 1 mm to its end face along the X direction, and set the degrees of freedom of the three corner points in the YZ direction, Y direction, and Z direction to 0 respectively (see Appendix Figure 5 ).

[0023] 5. Select the static general mode to compress the polygonal model at the microscopic scale, and set the total compression time to 1s.

[0024] 6. Create a new geometric model of the same size as the polygonal model (50×50×100), set four reference points (0,0,20), (50,0,20), (0,05,80), (50,50,80) respectively, perform three-point meshing on the model, insert cohesive units as bedding planes on the meshing plane, and specify the unit type as COH3D8 (see Appendix Figure 6 ). The bedding plane forms an inclination angle with the geometric model, which can be determined according to the bedding angle of the actual shale.

[0025] 7. Create a repeated mesh pattern in a specific area or direction for the geometric model to simulate periodic behavior. First, perform tetrahedral meshing on the geometric model, then offset and copy the mesh surface of the geometric model (see the attached flowchart). Figure 7 ), so that the meshes of the opposite faces are the same, then the mesh volume is divided to achieve the periodicity of the geometric model mesh, and the global mesh size is defined as 3 (see Appendix Figure 8), specify the model mesh attribute as C3D4.

[0026] 8. Map the periodic grid divided in the geometric model to the polygonal model. According to whether the center coordinates of the unit are in the corresponding geometric area, use coordinate comparison to map the grid (see Appendix Fig. 9 ), the mapped model is the representative volume element (RVE) model, which can provide sufficient accuracy to represent the larger scale of the material (see Appendix Fig.10 ).

[0027] 9. In order to eliminate the boundary effect caused by the scale limitation of the simulation system and make the stress and displacement continuous at the boundary, periodic boundary conditions are created for the geometric model. Use python script to add periodic boundary conditions: first, establish and read six boundary surface sets for the polygonal model, and save the node data on the six surfaces; second, read the node number and coordinates on each surface; third, calculate the nearest points on each pair of relative surfaces and save the node pair numbers; fourth, number the eight diagonal points and divide them into two groups, front and back, lower left, lower right, upper left, and upper right; fifth, establish Equation constraints for the node pairs on the left and right surfaces, upper and lower surfaces, and front and back surfaces to constrain the degrees of freedom in the X, Y, and Z directions; finally, establish Equation constraints for the eight diagonal points. Fig.11 The unit shown is for reference. Add periodic boundary conditions to achieve the corresponding load conditions through the following constraint equations. Taking the X direction as an example, after adding periodic boundaries, as shown in the attached figure Fig.12 .

[0028] Surface node constraint equations: On the opposite surface perpendicular to the X axis: ; On the opposite surface perpendicular to the Y axis: ; On the opposite surface perpendicular to the Z axis: ; Constraint equations between edge nodes: Three sets of linear constraint equations between the four edges parallel to the Z axis: ; Constraint equations between corner nodes: ; In the formula, ε 0 x , ε 0 y , ε 0z , γ 0 xy , γ 0 xz , γ 0 yz are the strains in the X, Y, Z, XY, XZ, and YZ directions respectively; x =W x , y =W y , z = h The three planes are called principal planes, and the planes parallel to the principal planes are called subordinate planes. x , W y , h They are the length, width and height of the cube polygonal model respectively; u , v , w They are x , y , z The displacement in the direction, the subscripts represent different corner nodes (D, E, F, G), edge nodes (EA, HD, FB, GC), surface nodes ( x =0, x =W x , y =0, y =W y , z =0, z = h ) position, that is u ∣ x=w for x direction x =W x The displacement at the surface node position, v EA for y The displacement at the edge node in the EA direction, w F for z Direction F displacement at the angle node.

[0029] 10. Compression and tension experiments were carried out on the polygonal model at the microscopic scale to solve the elastic mechanical properties such as stress-strain curve, compressive strength, elastic modulus and Poisson's ratio: First, the reaction forces received by all the collection points on the tensile surface were summed up using formula (17), and the stress-strain curve of the entire representative volume element was drawn using formulas (18) and (19). Then, according to the displacement deformation of the model in the X direction and the Y direction (or the Z direction), the mechanical parameters of the layered shale were obtained using formulas (20-23), and compared with the results of the indoor rock compression experiment (see Appendix 6), the compression cloud map of the layered shale and the black shale (see Appendix 10). Figure 13-14 ) and stress-strain curve (see Appendix Fig.15 ).

[0030] The sum of the reaction forces at all points on the extruded surface: ; Stress and strain in the X direction: ; Strain in Y and Z directions: ; Poisson's ratio: ; Modulus of elasticity: ; Where: rf i is the first i The reaction force of the point, RF x is the sum of all reaction forces on the plane; A x is the area of ​​the plane perpendicular to the X direction; u x , u y , u z are the displacements in the X, Y, and Z directions respectively, l x , l y , l z They are the maximum dimensions in the X, Y, and Z directions respectively; v is Poisson's ratio, σ x is the stress in the X direction, ε x , ε y , ε z are strains in the X, Y, and Z directions respectively; E is the elastic modulus, σ x1 , σ x2 , ε x1 , ε x2 are the stress and strain at any two points before the highest point of the stress-strain curve in the X direction.

[0031] Appendix 6 Comparison and verification of mechanical parameters of layered shale

[0032] 11. By the attached Figure 12-13 It can be seen that in the model of bedding shale under compressive load, the stress values ​​are mainly concentrated around the bedding plane, and shear failure first occurs along the bedding plane, that is, the area with relatively weak mechanical properties. In addition, the elastic modulus and tensile strength of bedding shale are also affected by its mineral composition and bedding mechanical properties. The mineral particles are not as tightly bound on the bedding plane as inside the layer, so the bedding has a smaller elastic modulus and tensile strength. This is because the bedding structure is relatively weak. These characteristics make bedding shale easier to break. Black shale can maintain a high ability to resist stretching when subjected to external forces. Appendix 6 lists the simulation results of rock mechanical parameters corresponding to the mineral content in Appendix 4 and the comparison with the experimental results. It can be seen that for black shale, the tensile strength of the experimental sample was measured to be 25.02MPa, and the elastic modulus was 13.28GPa. Fig.14 From the stress-strain curve, it can be seen that the stress required for bedding shale and black shale first increases and then slowly decreases, and the elastic modulus and compressive strength of bedding shale are both lower than those of black shale. In summary, the error of the mechanical parameters obtained by comparing the experiment and the simulation is within 5%. Therefore, this method, as a means of solving mechanical properties by combining micro-macro, meets the accuracy requirements in engineering. The effects of the physical properties, mineral proportions, arrangement patterns, and clay mineral hydration of bedding shale on the elastic mechanical properties of bedding shale can be further studied as needed.

[0033] 12. On a macro scale, geometric models with bedding angles of 0°, 30°, and 60° can be established in the bedding shale, and then the prediction can be solved according to the model establishment and calculation process in steps (6-10). The technical route of specific implementation methods 1-12 is shown in the attached Fig.16 shown.

Claims

1. A cross-scale calculation method for elastic mechanical properties of layered shale, comprising the following steps: 1) Establishment of Thiessen polygon model: The cell shape of the three-dimensional Thiessen polygon Voronoi shows good similarity with the shape of mineral particles. Therefore, the multi-cell of the Voronoi diagram is used to represent the structure of the internal particles of the rock. The Thiessen polygon is a set of continuous polygons composed of perpendicular bisectors connecting two adjacent point segments. The distance from any point in a Thiessen polygon to the polygon is less than the distance to the control points of other polygons. According to the Thiessen polygon coordinates, the Python script is used in the abaqus software to generate the Thiessen polygon model; 2) Assign different mineral grain volumes to the polygonal model according to the mineral ratio of the bedding shale, and calculate the grain volume of each Thiessen polygon; 3) Use Python script to randomly place mineral particles at any position of the polygonal model according to the calculated volume, and assign different attributes to the polygonal model containing mixed mineral particles; 4) The polygonal model is divided into tetrahedral meshes as a verification model, with a global mesh size of 3 and a mesh attribute of C3D4. By adding reference points and applying stress surface coupling, loads are applied to the rock, i.e., tension / compression, to create fixed boundary conditions. One end of the model is fixed, i.e., the degree of freedom of its end face in the X direction is set to 0, and the degrees of freedom of the three corner points in the XYZ direction, XY direction, and XZ direction are also set to 0. The other end is compressed, i.e., a displacement of 1 mm is added to its end face along the X direction, and the degrees of freedom of the three corner points in the YZ direction, Y direction, and Z direction are also set to 0; 5) Select the static general mode to compress the polygonal model at the microscopic scale, and set the total compression time to 1s; 6) Create a new geometric model of the same size as the polygonal model, i.e., 50×50×100. Set four reference points (0,0,20), (50,0,20), (0,05,80), and (50,50,80) respectively. Perform three-point segmentation on the model and insert cohesive units as bedding planes on the segmentation planes. Specify the unit type as COH3D8. The bedding plane forms an inclination angle with the geometric model, which can be determined according to the bedding angle of the actual shale. 7) Create a repeated mesh pattern for the geometric model in a specific area or direction to simulate periodic behavior; first, perform tetrahedral meshing on the geometric model, then offset and copy the mesh faces of the geometric model so that the meshes of the opposite faces are the same, then divide the mesh volume to achieve the periodicity of the geometric model mesh, define the global mesh size as 3, and specify the model mesh attribute as C3D4; 8) Map the periodic grid divided in the geometric model to the polygonal model. According to whether the center coordinates of the unit are in the corresponding geometric area, the grid is mapped by coordinate comparison. The mapped model is a representative volume element RVE model. RVE can provide sufficient accuracy to represent the larger scale of the material. 9) In order to eliminate the boundary effect caused by the scale limitation of the simulation system and make the stress and displacement at the boundary continuous, periodic boundary conditions are created for the geometric model; Python script is used to add periodic boundary conditions: First, six boundary surface sets are established and read for the polygonal model, and the node data on the six surfaces are saved; second, the node number and coordinates on each surface are read; third, the nearest points on each pair of relative surfaces are calculated, and the node pair numbers are saved; fourth, the eight diagonal points are numbered and divided into two groups, namely the lower left, lower right, upper left, and upper right; fifth, equation constraints are established for the node pairs on the left and right surfaces, the upper and lower surfaces, and the front and back surfaces, to constrain the degrees of freedom in the three directions of X, Y, and Z; finally, equation constraints are established for the eight diagonal points; 10) Perform compression and tensile tests on polygonal models at the microscopic scale to solve elastic mechanical properties such as stress-strain curve, compressive strength, elastic modulus and Poisson's ratio: 11) On a macro scale, geometric models of bedding in stratified shales with bedding angles of 0°, 30°, and 60° can be established respectively, and then solutions and predictions can be made based on the model building and calculation processes in steps 6) to 10).

2. The method for calculating elastic mechanical properties of cross-scale layered shale according to claim 1, characterized in that: In step 2), the volume is calculated using the Thiessen polygon area and elevation difference. First, the area of ​​each Thiessen polygon is calculated using GIS software. Then, the elevation difference between the highest point and the lowest point in each Thiessen polygon is calculated, such as formula 1-formula 7 to calculate the volume of each Thiessen polygon. ； Where l, w, h are the length, width and height of the model respectively; V and Vquartz, Vcalcite and Villite are the total volume of the polygonal model and the volume of quartz, calcite and illite mineral particles in the model respectively; wquartz is the volume fraction of quartz; v is the volume of any particle in the polygonal model; s is the area of ​​any particle in the polygonal model; Δh is the elevation difference between the highest point and the lowest point in any particle in the polygonal model.

3. The method for calculating elastic mechanical properties of bedding shale across scales according to claim 1, characterized in that: In step 9): Surface node constraint equations: On the opposite surface perpendicular to the X axis: ； On the opposite surface perpendicular to the Y axis: ； On the opposite surface perpendicular to the Z axis: ； Constraint equations between edge nodes: Three sets of linear constraint equations between the four edges parallel to the Z axis: ； Constraint equations between corner nodes: ； In the formula, ε 0 x , ε 0 y , ε 0 z , γ 0 xy , γ 0 xz , γ 0 yz are the strains in the X, Y, Z, XY, XZ, and YZ directions respectively; x =W x , y =W y , z = h The three planes of the system are called principal planes, and the planes parallel to the principal planes are called subordinate planes. x , W y , h They are the length, width and height of the cube polygonal model respectively; u , v , w They are x , y , z The displacement in the direction, the subscripts represent different corner nodes (D, E, F, G), edge nodes (EA, HD, FB, GC), surface nodes ( x =0, x =W x , y =0, y =W y , z =0, z = h ) position, that is u ∣ x=w for x direction x =W x The displacement at the surface node position, v EA for y The displacement at the edge node in the EA direction, w F for z Direction F displacement at the angle node.

4. The method for calculating elastic mechanical properties of layered shale across scales according to claim 1, characterized in that: In step 10), firstly, the reaction forces received by all the collection points on the tensile surface are summed using formula (17), and the stress-strain curve of the entire representative volume element is drawn using formulas (18) and (19). Then, according to the displacement deformation of the model in the X direction, Y direction or Z direction, the mechanical parameters of the bedding shale are obtained using formulas (20-23), and compared with the results of the indoor rock compression test, the compression cloud map and stress-strain curve of the bedding shale and black shale; The sum of the reaction forces at all points on the extruded surface: ； Stress and strain in the X direction: ； Strain in Y and Z directions: ； Poisson's ratio: ； Modulus of elasticity: ; Where: rf i is the first i The reaction force of the point, RF x is the sum of all reaction forces on the plane; A x is the area of ​​the plane perpendicular to the X direction; u x , u y , u z are the displacements in the X, Y, and Z directions respectively, l x , l y , l z They are the maximum dimensions in the X, Y, and Z directions respectively; v is Poisson's ratio, σ x is the stress in the X direction, ε x , ε y , ε z are strains in the X, Y, and Z directions respectively; E is the elastic modulus, σ x1 , σ x2 , ε x1 , ε x2 are the stress and strain at any two points before the highest point of the stress-strain curve in the X direction.

Images