A method for predicting crack morphology considering interlayer weak planes

Abstract

The present invention discloses a method for predicting fracture morphology considering interlayer weak planes, comprising the following steps: establishing a three-dimensional geological model according to the formation parameters of the actual interactive strata obtained; establishing a rock seepage-stress coupling model based on the stress balance equation, the continuity equation and the boundary conditions; determining the fracture extension and expansion criterion; inputting the fracturing well parameters of the target well area, establishing a three-dimensional full coupling model for fracturing; simulating the influence of different geological factors and engineering factors on fracture morphology. The present invention is based on the finite element method, establishes a full three-dimensional seepage-stress coupling model, uses the Cohesive pore pressure unit to preset two fracture surfaces, considers interlayer weak planes or bedding weak planes, studies the influence of different geological factors and engineering factors on fracture morphology, can quickly predict the fracture morphology during the fracturing of reservoirs containing interlayer weak planes and guide the design of fracturing construction parameters.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas exploration and development, and in particular to a method for predicting crack morphology taking interlayer weak surfaces into consideration. Background Art

[0002] As global oil and gas resource exploration and development gradually shifts to unconventional oil and gas development, hydraulic fracturing technology, as a measure to increase production of unconventional oil and gas wells, has been widely used in my country. However, unconventional oil and gas reservoirs are highly heterogeneous in the vertical direction. For example, tight sandstone has interlayers of sand and mudstone and mudstone interlayers, and shale has interlayers of shale-sandstone and shale-carbonate. This interlayer makes the fracture morphology complex and difficult to predict during hydraulic fracturing. How to predict the fracture morphology based on the actual interlayer interaction of the formation to guide the optimization of construction parameters is one of the current difficulties.

[0003] In order to study the extension morphology of hydraulic fractures in alternating strata, many scholars currently mainly establish models from two-dimensional and three-dimensional perspectives for research.

[0004] In terms of two-dimensional model research: The patent document with application number 202010545763.X discloses a multi-cluster fracturing fracture morphology prediction method based on discrete element fluid-solid coupling. This method can effectively consider the influence of interlayer weak surfaces on fracture morphology, but does not consider the flow of liquid in the third dimension. There is a certain gap with the actual three-dimensional formation and it has certain limitations.

[0005] In terms of three-dimensional model research, the patent document with application number 201910470372.3 discloses a method for constructing a cross-layer fracturing map. This method establishes a three-dimensional model and studies the height extension law of fractures in sandstone-mudstone alternating reservoirs, but does not consider the influence of interlayer weak surfaces on fracture morphology.

[0006] The existing methods that can consider the fracture morphology of interlayer weak surfaces are limited to two-dimensional models, which is somewhat different from the actual three-dimensional strata. The three-dimensional model is difficult to consider the influence of interlayer weak surfaces on fracture morphology. Summary of the invention

[0007] The purpose of the present invention is to address the deficiencies in the prior art and provide a method for predicting fracture morphology by taking into account interlayer weak surfaces. By establishing a full three-dimensional seepage-stress coupling model and considering the influence of interlayer weak surfaces on hydraulic fracturing, fracture morphology prediction is performed.

[0008] The technical solution adopted by the present invention to solve the above technical problems includes the following contents:

[0009] A method for predicting crack morphology considering interlayer weak planes comprises the following steps:

[0010] S1. Establishing a three-dimensional geological model based on the stratigraphic parameters of the actual interactive strata obtained;

[0011] S2. Based on the stress equilibrium equation, continuity equation and boundary conditions, a rock seepage-stress coupling model is established;

[0012] S3. Determine the crack extension and propagation criterion;

[0013] S4. Input the parameters of the target well area fracturing well and establish a three-dimensional fully coupled fracturing model;

[0014] S5. Simulate the influence of different geological factors and engineering factors on fracture morphology;

[0015] The method for establishing the rock seepage-stress coupling model in step S2 is:

[0016] (1) Based on porous media theory and the principle of virtual work, the rock stress equilibrium equation is established:

[0017]

[0018] In the formula, is the effective stress borne by the solid skeleton, Pa; P nw is the non-wetting phase pressure, Pa; I is [1,1,1,0,0,0] T The unit matrix; δε is the imaginary strain rate matrix, s -1 ; t is the surface force matrix, N / m 3 ; δv is the virtual velocity matrix, m / s; f is the body force matrix, N / m 3 ; dV is the unit cell, m 3 ; S is the integral space surface, m 2 ;

[0019] (2) Based on the principle of conservation of matter, it is assumed that the flow of fluid in rocks satisfies Darcy's law, and the fluid continuity equation is established:

[0020]

[0021] Where J is the volume change rate of the porous medium, dimensionless; n ω is the ratio of the fluid volume to the total volume, dimensionless; ρ ω is the fluid density, kg / m 3 ; x is the space vector, m; dt is the time increment, s; v ω is the seepage velocity, m / s;

[0022] (3) Set the flow boundary conditions, pore pressure boundary conditions and position boundary conditions respectively:

[0023] The flow boundary condition is:

[0024]

[0025] Where n is the unit normal direction of the flow boundary; k is the permeability coefficient tensor; q0 is the total amount of liquid flowing through the boundary per unit time;

[0026] Among them, the pore pressure boundary condition can be expressed as P0=P 0b It means that the pore pressure on the boundary is a certain value P 0b ;

[0027] Among them, the position boundary condition is used to constrain the displacement of the boundary nodes in the X, Y, and Z directions. x =0,U y =0,U z =0;

[0028] (4) By defining the shape function and using the Galerkin method for finite element discretization, the seepage-stress coupling equation can be obtained:

[0029]

[0030] Where: K = ∫ V B T D ep BdV;

[0031] F = ∫ V (▽N p ) T kk r ▽N p dV;

[0032]

[0033] Among them, N u , B is the defined shape function vector matrix; is the node displacement; is the pore pressure at the unit node; N p is the shape function; D ep is the elasticity matrix; m is the transformation matrix; q ob is the fluid flow rate on the boundary;

[0034] The rock seepage-stress coupling model includes a rock physical property parameter dynamic model, and the rock physical property parameter dynamic model is:

[0035] (1) Multi-field coupled dynamic porosity model:

[0036]

[0037] Where ω is the coefficient of thermal expansion, 1 / k; ε v is the volume strain, dimensionless;

[0038] (2) Permeability dynamic model derived based on the Kozeyn-Cannan capillary bundle model:

[0039]

[0040] In the formula, ε v is the volume strain, dimensionless; ω is the thermal expansion coefficient, 1 / k; ΔT is the temperature change, K;

[0041] (3) Dynamic model of pore compressibility:

[0042]

[0043] In the formula, ε v is the volume strain, dimensionless; ΔT is the temperature change, K; Δp is the pressure change, MPa;

[0044] In step S3, the fracture extension criterion includes using a viscoelastic Cohesive pore pressure unit to preset two fracture surfaces to simulate fracture initiation and extension, fracture surface I is perpendicular to the minimum principal stress direction, fracture surface II exists at the junction of the stratum layer, and fracture surface I and fracture II intersect;

[0045] There are four pore pressure points at the unit where the fracture surface I and the fracture surface II intersect. By replacing the four pore pressure nodes with a common pore pressure node, the continuity condition of the fluid pressure at the intersection is met;

[0046] The crack initiation conditions are:

[0047]

[0048] In the formula, σ n is the applied stress in the normal direction of the Colesive unit, MPa; τ s , τ t is the applied stress in two tangential directions of the unit, MPa; is the critical stress in the normal direction when the unit fails, MPa; is the critical stress in two directions when the unit fails in the tangential direction, MPa;

[0049] The viscoelastic Cohesive pore pressure unit will enter the damage evolution stage after the crack damage, that is, the damage will continue to occur until the unit is completely destroyed when the fracture energy is reached, and finally a crack is formed. The damage factor calculation formula is:

[0050]

[0051] In the formula, is the maximum displacement reached by the unit during loading, m; is the displacement when the unit is completely destroyed, m; is the displacement of the unit at the initial damage, m;

[0052] The damage factor D is used to determine the damage of the unit: when the damage has not yet occurred, the D value is 0; when the unit is completely damaged, the D value is 1;

[0053] The injected liquid during fracturing is a Newtonian fluid. The flow model simulating the fracturing fluid in the fracture surface includes the tangential flow along the fracture surface and the normal flow perpendicular to the fracture surface:

[0054] (1) The tangential flow along the fracture surface is defined by Newtonian flow pressure conduction as:

[0055] Q=-k t Δp

[0056] Where Q is the displacement of fracturing fluid, m 3 / s;k t is the flow coefficient; p is the flow pressure, MPa;

[0057] Among them, the flow coefficient k t It can be expressed as:

[0058]

[0059] Where, d is the crack opening width; μ is the fracturing fluid viscosity, Pa.s;

[0060] (2) The normal flow perpendicular to the fracture surface simulates the seepage loss in the normal direction of the fracture surface and is characterized by defining the loss coefficient. The normal flow is expressed as:

[0061]

[0062] In the formula, q t ,q b are the seepage flow rates of the fluid on the upper and lower surfaces of the fracture, m 3 / s;c t 、c b are the fluid loss coefficients on the upper and lower surfaces of the fracture respectively; p t 、p b are the pore pressures of the body on the upper and lower surfaces of the crack, MPa; p i is the fluid pressure at the middle surface of the Colesive fracture unit, MPa.

[0063] Preferably, in step S5, the geological factors include any one or more of reservoir thickness, barrier thickness, number of layers, tensile strength of layer interfaces, rock physical parameters, rock mechanical parameters and stress difference; the engineering factors include any one or more of construction displacement, perforation orientation, fracturing fluid viscosity, proppant concentration and construction scale.

[0064] The beneficial effects of the present invention are:

[0065] The present invention is based on the finite element method, establishes a full three-dimensional seepage-stress coupling model, uses the Cohesive pore pressure unit to preset two types of fracture surfaces, considers the interlayer weak surface or the bedding weak surface, studies the influence of different geological factors and engineering factors on the fracture morphology, and can guide the design of fracturing construction parameters. It has preliminarily solved the current problem that it is difficult to study the influence of the interlayer weak surface or the bedding weak surface on the fracture extension from a three-dimensional perspective. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 It is a schematic diagram of the implementation flow of a method for predicting crack morphology considering interlayer weak planes according to an embodiment of the present invention;

[0067] Figure 2 is a schematic diagram of a three-dimensional geological model of shale-sandstone alternating strata established according to an embodiment of the present invention;

[0068] Figure 3 is a schematic diagram of intersecting fracture planes I and II preset according to an embodiment of the present invention;

[0069] Figure 4 This is a schematic diagram of the merging process of the intersection nodes of the crack surface I and the crack surface II;

[0070] Figure 5 is a schematic diagram of a three-dimensional fully coupled fracturing model established according to an embodiment of the present invention;

[0071] Figure 6 is a schematic diagram of crack morphology calculated according to an embodiment of the present invention;

[0072] Figure 7 is a schematic diagram of the relationship between the seam height and the tensile strength of the layer interface calculated according to an embodiment of the present invention;

[0073] Figure 8 is a schematic diagram of the relationship between the joint height and the construction displacement calculated according to an embodiment of the present invention;

[0074] Fig. 9 Schematic diagram of the relationship between fracture height and fracturing fluid viscosity calculated according to an embodiment of the present invention. DETAILED DESCRIPTION

[0075] The following will describe the implementation methods of the present invention in detail with reference to the accompanying drawings and embodiments, so that the implementation process of how the present invention applies technical means to solve technical problems and achieve technical effects can be fully understood and implemented accordingly. It should be noted that as long as there is no conflict, the various embodiments of the present invention and the various features in the embodiments can be combined with each other, and the technical solutions formed are all within the protection scope of the present invention.

[0076] Meanwhile, in the following description, many specific details are set forth for the purpose of explanation to provide a thorough understanding of the embodiments of the present invention. However, it is obvious to those skilled in the art that the present invention can be implemented without the specific details or the specific manner described herein.

[0077] In addition, the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be performed in an order different from that shown here.

[0078] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0079] A method for predicting crack morphology considering interlayer weak planes comprises the following steps:

[0080] S1. Establishing a three-dimensional geological model based on the stratigraphic parameters of the actual interactive strata obtained;

[0081] S2. Based on the stress equilibrium equation, continuity equation and boundary conditions, a rock seepage-stress coupling model is established;

[0082] S3. Determine the crack extension and propagation criterion;

[0083] S4. Input the parameters of the target well area fracturing well and establish a three-dimensional fully coupled fracturing model;

[0084] S5. Simulate the impact of different geological and engineering factors on crack morphology.

[0085] Figure 1 A schematic diagram of the implementation flow of the crack morphology prediction method considering interlayer weak planes provided in this embodiment is shown.

[0086] In this embodiment, the Abaqus finite element simulation calculation tool is preferably used as a platform to establish a three-dimensional geological model, a rock seepage-stress coupling model, and two types of fracture surfaces are preset using the Cohesive pore pressure unit, and the fracture initiation and expansion are simulated according to the corresponding fracture extension rules. Of course, in other embodiments of the present invention, according to actual needs, the method can also use other reasonable tools or methods to realize the fracture morphology prediction of the reservoir containing interlayer weak planes.

[0087] In this embodiment, the method preferably first establishes a three-dimensional geological model of shale-sandstone alternating strata in step S1. The geological characteristics of the shale-sandstone alternating strata are: 6 shale layers, 5 sandstone layers, and the perforation section is located in the middle shale layer. The three-dimensional geological model is established, from top to bottom, shale (4m), sandstone (4m), shale (10m), sandstone (4m), shale (8m), sandstone (4m), shale (4m), sandstone (4m), shale (2m), sandstone (4m), shale (2m), Figure 2 shown.

[0088] The method for establishing the rock seepage-stress coupling model in step S2 is:

[0089] (1) Based on porous media theory and the principle of virtual work, the rock stress equilibrium equation is established:

[0090]

[0091] In the formula, is the effective stress borne by the solid skeleton, Pa; P nw is the non-wetting phase pressure, Pa; I is [1,1,1,0,0,0] T The unit matrix; δε is the imaginary strain rate matrix, s -1 ; t is the surface force matrix, N / m 3 ; δv is the virtual velocity matrix, m / s; f is the body force matrix, N / m 3 ; dV is the unit cell, m 3 ; S is the integral space surface, m 2 ;

[0092] (2) Based on the principle of conservation of matter, it is assumed that the flow of fluid in rocks satisfies Darcy's law, and the fluid continuity equation is established:

[0093]

[0094] Where J is the volume change rate of the porous medium, dimensionless; n ω is the ratio of the fluid volume to the total volume, dimensionless; ρ ω is the fluid density, kg / m 3 ; x is the space vector, m; dt is the time increment, s; v ω is the seepage velocity, m / s;

[0095] (3) Set the flow boundary conditions, pore pressure boundary conditions and position boundary conditions respectively:

[0096] The flow boundary condition is:

[0097]

[0098] Where n is the unit normal direction of the flow boundary; k is the permeability coefficient tensor; q0 is the total amount of liquid flowing through the boundary per unit time;

[0099] Among them, the pore pressure boundary condition can be expressed as P0=P 0b It means that the pore pressure on the boundary is a certain value P 0b ;

[0100] Among them, the position boundary condition is used to constrain the displacement of the boundary nodes in the X, Y, and Z directions. x =0,U y =0,U z =0;

[0101] (4) By defining the shape function and using the Galerkin method for finite element discretization, the seepage-stress coupling equation can be obtained:

[0102]

[0103] Where: K = ∫ V B T D ep BdV;

[0104] F = ∫ V (▽N p ) T kk r ▽N p dV;

[0105]

[0106] Among them, N u , B is the defined shape function vector matrix; is the node displacement; is the pore pressure at the unit node; N p is the shape function; D ep is the elasticity matrix; m is the transformation matrix; q ob is the fluid flow rate on the boundary.

[0107] As the fracturing fluid is injected, the change in fluid pressure will cause the rock physical parameters to change. It is necessary to establish a dynamic model of rock physical parameters to update the rock physical parameter data in the rock seepage-stress coupling model. The dynamic model of rock physical parameters is:

[0108] (1) Multi-field coupled dynamic porosity model:

[0109]

[0110] Where ω is the coefficient of thermal expansion, 1 / k; εv is the volume strain, dimensionless;

[0111] (2) Permeability dynamic model derived based on the Kozeyn-Cannan capillary bundle model:

[0112]

[0113] In the formula, ε v is the volume strain, dimensionless; ω is the thermal expansion coefficient, 1 / k; ΔT is the temperature change, K;

[0114] (3) Dynamic model of pore compressibility:

[0115]

[0116] In the formula, ε v is the volume strain, dimensionless; ΔT is the temperature change, K; Δp is the pressure change, MPa.

[0117] In step S3, the fracture extension criterion preferably uses a viscoelastic Cohesive pore pressure unit to preset two fracture surfaces to simulate fracture initiation and extension, fracture surface I is perpendicular to the minimum principal stress direction, fracture surface II exists at the junction of the stratum layer, and fracture surface I and fracture II intersect, such as Figure 3 shown.

[0118] like Figure 4 As shown in the figure, there are four pore pressure points at the unit where the fracture surface I and the fracture surface II intersect. By replacing the four pore pressure nodes with a common pore pressure node, that is, replacing the P1, P2, P3, and P4 pore pressure nodes in the intersecting unit with the common pore pressure node P0, the continuity condition of the fluid pressure at the intersection is met.

[0119] The prefabricated Cohesive pore pressure unit simulates the cracks, and the crack initiation conditions are:

[0120]

[0121] In the formula, σ n is the applied stress in the normal direction of the Colesive unit, MPa; τ s , τ t is the applied stress in two tangential directions of the unit, MPa; is the critical stress in the normal direction when the unit fails, MPa; is the critical stress in two directions when the unit fails tangentially, MPa.

[0122] When the Cohesive unit is damaged by cracking, it will enter the damage evolution stage, that is, the damage will continue to occur until the crack fracture energy is reached, and the unit will be completely destroyed, and finally a crack will be formed. The damage factor calculation formula is:

[0123]

[0124] In the formula, is the maximum displacement reached by the unit during loading, m; is the displacement when the unit is completely destroyed, m; is the displacement of the unit at the initial damage, m;

[0125] The damage factor D is used to determine the damage of the unit: when the damage has not yet occurred, the D value is 0; when the unit is completely damaged, the D value is 1.

[0126] In this embodiment, the fracturing fluid is regarded as a Newtonian fluid, and the flow model simulating the fracturing fluid in the fracture surface includes a tangential flow along the fracture surface and a normal flow perpendicular to the fracture surface:

[0127] (1) The tangential flow along the fracture surface is defined by Newtonian flow pressure conduction as:

[0128] Q=-k t Δp

[0129] Where Q is the displacement of fracturing fluid, m 3 / s;k t is the flow coefficient; p is the flow pressure, MPa;

[0130] Among them, the flow coefficient k t It can be expressed as:

[0131]

[0132] Where, d is the crack opening width; μ is the fracturing fluid viscosity, Pa.s;

[0133] (2) The normal flow perpendicular to the fracture surface simulates the seepage loss in the normal direction of the fracture surface and is characterized by defining the loss coefficient. The normal flow is expressed as:

[0134]

[0135] In the formula, q t ,q b are the seepage flow rates of the fluid on the upper and lower surfaces of the fracture, m 3 / s;c t 、c b are the fluid loss coefficients on the upper and lower surfaces of the fracture respectively; p t 、p b are the pore pressures of the body on the upper and lower surfaces of the crack, MPa; pi is the fluid pressure at the middle surface of the Colesive fracture unit, MPa.

[0136] Step S4, input the fracturing well parameters of the target well area, and establish a three-dimensional fully coupled fracturing model.

[0137] like Figure 5 As shown, the statistical parameters of the target well area fracturing well are input, including rock elastic modulus, Poisson's ratio, rock body tensile strength, rock body shear strength, prefabricated fracture surface tensile strength, ground stress, injection displacement, fracturing fluid viscosity, initial permeability, initial porosity, and a three-dimensional fully coupled fracturing model is established.

[0138] Table 1 Parameters of fracturing wells in target well area

[0139]

[0140]

[0141] The Abaqus implicit solver is used to solve the three-dimensional fully coupled fracturing model. Figure 6 As shown, cracks started to occur at the weak surface between the layers, and an "I-shaped crack" was clearly formed. The predicted half-crack length was 44m and the crack height was 10m.

[0142] Step S5, simulating the influence of different geological factors and engineering factors on the fracture morphology.

[0143] (1) By changing the geological factors and the tensile strength of the shale-sandstone interface, the crack extension was simulated under tensile strengths of 0.5, 1.0, 1.5, 2.0, and 3.0 MPa. The results are as follows: Figure 7 As shown:

[0144] The simulation results show that the smaller the tensile strength of the layer interface, the easier it is for the layer interface to break, and the easier it is for the crack to extend along the layer interface. When the tensile strength of the layer interface exceeds 2MPa, the layer interface will break in a small range, and the crack height will easily penetrate the layer interface.

[0145] (2) Changing the engineering factors and construction displacement, the construction displacements of 8, 10, 12, 14, and 16 m were simulated respectively. 3 / min crack expansion, the results are as follows Figure 8 As shown:

[0146] The simulation results show that as the construction displacement increases, the cracks continue to communicate with the weak surfaces between layers and the crack height increases; when the displacement increases by 1 times, the crack height increases by 25%. The displacement has a certain impact on the crack height.

[0147] (3) By changing the engineering factors and the viscosity of the fracturing fluid, the crack propagation was simulated under the conditions of fracturing fluid viscosity of 1, 10, 20, 50, 100, and 150 mPa.s. The results are shown in Figure 2. Fig. 9 shown.

[0148] The simulation results show that as the viscosity increases, the fracture height increases. If the viscosity is too high, the fracture height will break through the layer interface and the layer interface will rupture. The viscosity of the fracturing fluid has a great influence on the fracture morphology.

[0149] The "one embodiment" or "embodiment" mentioned in the specification means that the specific features, structures or characteristics described in conjunction with the embodiment are included in at least one embodiment of the present invention. Therefore, the phrases "one embodiment" or "embodiment" appearing in various places throughout the specification do not necessarily refer to the same embodiment.

[0150] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment as above, it is not used to limit the present invention. Any technician familiar with this profession can make some changes or modifications to equivalent embodiments of equivalent changes using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention still falls within the scope of the technical solution of the present invention.

Claims

1. A method for predicting crack morphology considering interlayer weak planes, characterized in that: The following steps are involved: S1. Establishing a three-dimensional geological model based on the stratigraphic parameters of the actual interactive strata obtained; S2. Based on the stress equilibrium equation, continuity equation and boundary conditions, a rock seepage-stress coupling model is established; S3. Determine the crack extension and propagation criterion; S4. Input the parameters of the target well area fracturing well and establish a three-dimensional fully coupled fracturing model; S5. Simulate the influence of different geological factors and engineering factors on fracture morphology; The method for establishing the rock seepage-stress coupling model in step S2 is: (1) Based on porous media theory and the principle of virtual work, the rock stress equilibrium equation is established: In the formula, is the effective stress borne by the solid skeleton, Pa; P nw is the non-wetting phase pressure, Pa; I is [1,1,1,0,0,0] T Identity matrix; δε is the virtual strain rate matrix, s -1 ; t is the surface force matrix, N / m 3 ;δv is the virtual velocity matrix, m / s; f is the physical strength matrix, N / m 3 ; dV is the unit cell, m 3 ; S is the integral space surface, m 2 ; (2) Based on the principle of conservation of matter, it is assumed that the flow of fluid in rocks satisfies Darcy's law, and the fluid continuity equation is established: Where J is the volume change rate of the porous medium, dimensionless; n ω is the ratio of the fluid volume to the total volume, dimensionless; ρ ω is the fluid density, kg / m 3 ; x is the space vector, m; dt is the time increment, s; v ω is the seepage velocity, m / s; (3) Set the flow boundary conditions, pore pressure boundary conditions and position boundary conditions respectively: The flow boundary condition is: Where n is the unit normal direction of the flow boundary; k is the permeability coefficient tensor; q0 is the total amount of liquid flowing through the boundary per unit time; Among them, the pore pressure boundary condition can be expressed as P0=P 0b It means that the pore pressure on the boundary is a certain value P 0b ; Among them, the position boundary condition is used to constrain the displacement of the boundary nodes in the X, Y, and Z directions. Let U x =0,U y =0,U z =0; (4) By defining the shape function and using the Galerkin method for finite element discretization, the seepage-stress coupling equation can be obtained: Where: K = ∫ V B T D ep BdV; Among them, N u , B is the defined shape function vector matrix; is the node displacement; is the pore pressure at the unit node; N p is the shape function; D ep is the elasticity matrix; m is the transformation matrix; q ob is the fluid flow rate on the boundary; The rock seepage-stress coupling model includes a rock property parameter dynamic model, and the rock property parameter dynamic model is: (1) Multi-field coupled dynamic porosity model: Where ω is the coefficient of thermal expansion, 1 / k; ε v is the volume strain, dimensionless; (2) Permeability dynamic model derived based on the Kozeyn-Cannan capillary bundle model: In the formula, ε v is the volume strain, dimensionless; ω is the thermal expansion coefficient, 1 / k; ΔT is the temperature change, K; (3) Dynamic model of pore compressibility: In the formula, ε v is the volume strain, dimensionless; ΔT is the temperature change, K; Δp is the pressure change, MPa; In step S3, the fracture extension criterion includes using a viscoelastic Cohesive pore pressure unit to preset two fracture surfaces to simulate fracture initiation and extension, fracture surface I is perpendicular to the minimum principal stress direction, fracture surface II exists at the junction of the stratum layer, and fracture surface I and fracture II intersect; There are four pore pressure points at the unit where the fracture surface I and the fracture surface II intersect. By replacing the four pore pressure nodes with a common pore pressure node, the continuity condition of the fluid pressure at the intersection is met; The crack initiation conditions are: In the formula, σ n is the applied stress in the normal direction of the Colesive unit, MPa; τ s , τ t is the applied stress in two tangential directions of the unit, MPa; is the critical stress in the normal direction when the unit fails, MPa; is the critical stress in two directions when the unit fails in the tangential direction, MPa; The viscoelastic Cohesive pore pressure unit will enter the damage evolution stage after the crack damage, that is, the damage will continue to occur until the unit is completely destroyed when the fracture energy is reached, and finally a crack is formed. The damage factor calculation formula is: In the formula, is the maximum displacement reached by the unit during loading, m; is the displacement when the unit is completely destroyed, m; is the displacement of the unit at the initial damage, m; The damage factor D is used to determine the damage of the unit: when the damage has not yet occurred, the D value is 0; when the unit is completely damaged, the D value is 1; The injected liquid during fracturing is a Newtonian fluid. The flow model simulating the fracturing fluid in the fracture surface includes the tangential flow along the fracture surface and the normal flow perpendicular to the fracture surface: (1) The tangential flow along the fracture surface is defined by Newtonian flow pressure conduction as: Q=-k t Δp Where Q is the displacement of fracturing fluid, m 3 / s;k t is the flow coefficient; p is the flow pressure, MPa; Among them, the flow coefficient k t It can be expressed as: Where, d is the crack opening width; μ is the fracturing fluid viscosity, Pa.s; (2) The normal flow perpendicular to the fracture surface simulates the seepage loss in the normal direction of the fracture surface and is characterized by defining the loss coefficient. The normal flow is expressed as: In the formula, q t ,q b are the seepage flow rates of the fluid on the upper and lower surfaces of the fracture, m 3 / s;c t 、c b are the fluid loss coefficients on the upper and lower surfaces of the fracture respectively; p t 、p b are the pore pressures of the body on the upper and lower surfaces of the crack, MPa; p i is the fluid pressure at the middle surface of the Colesive fracture unit, MPa.

2. A method for predicting crack morphology considering interlayer weak planes according to claim 1, characterized in that: In step S5, the geological factors include any one or more of reservoir thickness, barrier thickness, number of layers, tensile strength of layer interfaces, rock physical parameters, rock mechanical parameters and stress difference; the engineering factors include any one or more of construction displacement, perforation orientation, fracturing fluid viscosity, proppant concentration and construction scale.

Images