A shale gas supercritical carbon dioxide fracturing and production integrated simulation method

Abstract

The application discloses a shale gas supercritical carbon dioxide fracturing and production integrated simulation method and belongs to the fracturing technical field of petroleum engineering. The method comprises the following steps: S1, setting a fracture potential path based on natural fractures and in-situ stress, generating a fracture grid from the fracture potential path, and performing value assignment on the fracture grid; S2, simulating fracture deformation and fracture expansion by adopting a plane three-dimensional displacement discontinuity method, combining an embedded discrete fracture model and a component model to represent matrix-fracture fluid exchange, and performing supercritical carbon dioxide fracturing simulation through iterative coupling; and S3, performing production simulation with a pressure field and component distribution after fracturing simulation as initial conditions according to the opening of the fractures after fracturing simulation and the permeability calculated from the opening. The application can more accurately realize supercritical carbon dioxide fracturing and production integrated simulation, can provide technical support for tracking the migration path of carbon dioxide in the fracturing-well killing-production whole cycle and evaluating the feasibility of carbon dioxide storage after fracturing.

Description

Technical Field

[0001] This invention relates to the field of fracturing technology in petroleum engineering, and in particular to an integrated simulation method for supercritical carbon dioxide fracturing production of shale gas. Background Technology

[0002] Waterless or low-water fracturing technology is gaining increasing attention due to its effectiveness in avoiding shale water-locking, reducing reservoir damage, and conserving water resources. Supercritical carbon dioxide fracturing technology within waterless fracturing not only helps reduce fracture pressure and create complex fracture networks but also contributes to achieving dual-carbon targets. Compared to conventional water-based fracturing fluids, supercritical carbon dioxide is highly compressible and has low viscosity. Conventional fracturing simulation software struggles to reflect the changes in carbon dioxide density, viscosity, and other physical properties with pressure during fracturing. Furthermore, due to its low viscosity, supercritical carbon dioxide readily mixes with reservoir fluids, resulting in a significantly higher filtration rate than conventional water-based fracturing fluids. The carbon dioxide filtration into the matrix causes changes in pressure and composition in the near-fracture region, impacting subsequent production. However, conventional integrated fracturing production simulations only import the fracture geometry after simulation into subsequent production modeling, neglecting the changes in matrix pressure and composition in the near-fracture region after fracturing. Summary of the Invention

[0003] To address the aforementioned problems, this invention aims to provide an integrated simulation method for shale gas supercritical carbon dioxide fracturing production.

[0004] The technical solution of the present invention is as follows:

[0005] An integrated simulation method for shale gas supercritical carbon dioxide fracturing production includes the following steps:

[0006] S1: Based on natural cracks and in-situ stress, set potential crack paths, generate crack meshes from the potential crack paths, and assign values ​​to the crack meshes;

[0007] S2: The planar three-dimensional displacement discontinuity method is used to simulate crack deformation and crack propagation. The embedded discrete crack model and component model are combined to characterize the matrix-crack fluid exchange. Supercritical carbon dioxide fracturing simulation is carried out through iterative coupling.

[0008] S3: Based on the fracture opening and the permeability calculated from it after the fracturing simulation, production simulation is carried out using the pressure field and composition distribution after the fracturing simulation as initial conditions.

[0009] Preferably, in step S1, when assigning values ​​to the fracture mesh, the initial fracture mesh pressure and initial fracture aperture are given, the initial permeability of the fracture mesh is made equal to the matrix permeability, and all fracture meshes are saturated with carbon dioxide at the initial moment.

[0010] Preferably, in step S1, after generating the crack mesh, the step further includes a step of refining the matrix mesh containing the crack mesh.

[0011] Preferably, in step S2, when using the planar three-dimensional displacement discontinuity method to simulate crack deformation, the expression for crack deformation is:

[0012] (1)

[0013] In the formula: and Let m represent the tangential displacement discontinuities parallel to and orthogonal to the crack surface. This represents the discontinuity of normal displacement perpendicular to the crack surface, which is equal to the crack mesh aperture. m; superscript , Indicates the crack grid number; and These represent the in-situ stress components acting tangentially to and orthogonal to the crack surface, in MPa; This represents the in-situ stress component acting in the direction perpendicular to the crack surface, expressed in MPa. This represents the crack mesh pressure, in MPa. Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress in the same direction of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress of crack mesh i in the orthogonal direction, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress of crack mesh i in the orthogonal direction, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress in the same direction of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the tangential stress of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the tangential stress in the other direction of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Indicates the total number of open crack meshes;

[0014] When using the planar three-dimensional displacement discontinuity method to simulate crack propagation, the expression for judging crack propagation is:

[0015] (2)

[0016] In the formula: This represents the Type I stress concentration factor of the crack mesh. ; This represents the Young's modulus of the rock, in MPa. Indicates the Poisson's ratio of the rock; The crack mesh half-length is represented in meters (m). Indicates the type I fracture toughness of rock. .

[0017] Preferably, in step S2, the component model can characterize the multi-component competitive adsorption mechanism and porous media diffusion mechanism in shale gas;

[0018] The multi-component competitive adsorption mechanism is as follows:

[0019] (3)

[0020] In the formula: This represents the adsorption capacity of component i, in kg / m³. 3 ; This represents the Langmuir adsorption capacity of component i, in kg / m³. 3 ; Represents the mole fraction of component i; Indicates matrix pore pressure, in MPa; The Langmuir adsorption pressure of component i is expressed in MPa. Indicates the total number of components;

[0021] The diffusion mechanism of the porous medium is as follows:

[0022] (4)

[0023] In the formula: This represents the diffusion flux of component i in phase α, in kg / (m²). 2 ·s); Indicates the porosity of the matrix rock; Indicates the saturation degree of the α phase; Indicates the tortuosity of the matrix rock; m represents the diffusion coefficient of component i in phase α. 2 / s; Represents the gradient operator; The density of the α phase is expressed in kg / m³. 3 ; This represents the mole fraction of component i in the α phase.

[0024] Preferably, step S2, which involves iterative coupling to simulate supercritical carbon dioxide fracturing, specifically includes the following sub-steps:

[0025] S21: Based on the initial crack mesh pressure or the crack mesh pressure calculated in the previous time, give an initial value of crack mesh pressure, and use the planar three-dimensional displacement discontinuity method to update the opening degree of the opened crack mesh based on this initial pressure value.

[0026] S22: The equivalent carbon dioxide injection time Δt is obtained by calculating the crack mesh aperture updated in step S21. f ;

[0027] S23: Import the fracture mesh aperture obtained in step S21 and the fracture mesh permeability calculated therefrom into the reservoir simulator, and use the component model to simulate and obtain dt+Δt. f Crack mesh pressure after time;

[0028] S24: Determine whether the crack pressure converges based on the initial crack mesh pressure given in step S21 and the crack mesh pressure obtained by simulation in step S23.

[0029] If convergence occurs, the crack pressure cycle is complete, and the process proceeds to step S25.

[0030] If convergence is not achieved, update the initial value of the crack mesh pressure and repeat steps S21-S24.

[0031] S25: Determine whether the crack tip is extending forward;

[0032] If so, then activate the corresponding crack mesh and repeat steps S21-S25;

[0033] If not, the crack propagation cycle is completed, and the process proceeds to step S26;

[0034] S26: Determine whether the current fracturing simulation time has reached the fracturing end time;

[0035] If so, then complete the fracturing simulation;

[0036] If not, repeat steps S21-S26.

[0037] Preferably, in step S21, the initial value of the crack mesh pressure is given by the following formula:

[0038] (5)

[0039] In the formula: This represents the initial pressure value of the given crack mesh, in MPa; This represents the initial fracture mesh pressure or the fracture mesh pressure calculated at the previous time step, in MPa; Indicates the parameters of crack pressure disturbance;

[0040] In step S24, the initial value of the crack mesh pressure is updated using the following formula:

[0041] (6)

[0042] In the formula: This indicates the updated crack mesh pressure; This represents the iteration parameter, whose value is between 0 and 1; This represents the crack mesh pressure obtained from step S23 of the simulation, in MPa.

[0043] Preferably, in step S22, the equivalent carbon dioxide injection time Δt is... f The calculation is performed using the following formula:

[0044] (7)

[0045] In the formula: The volume of carbon dioxide injected into the ground per unit time is expressed in m. 3 / d; This represents the ground density of carbon dioxide, in kg / m³. 3 ; This represents the crack mesh aperture obtained from step S21, in meters (m). Indicates the initial crack mesh aperture or the crack mesh aperture updated in the previous time step, m; This represents the initial volume of the crack mesh or the volume of the crack mesh updated in the previous time step, m. 3 ; This represents the initial carbon dioxide density in the fractured mesh or the carbon dioxide density after the previous time step update, in kg / m³. 3 .

[0046] Preferably, in step S23, the permeability of the crack mesh is calculated using the following formula:

[0047] (8)

[0048] In the formula: The permeability of the crack mesh is expressed in mD. The crack mesh aperture is represented in meters (m).

[0049] In step S24, the expression for determining whether the crack pressure has converged is:

[0050] (9)

[0051] In the formula: This represents the crack mesh pressure obtained from step S23 of the simulation, in MPa. This represents the initial value of the slot mesh pressure given in step S21, in MPa; The convergence tolerance is set by the user.

[0052] Preferably, in step S3, the permeability is calculated using the following formula:

[0053] (10)

[0054] In the formula: Permeability, expressed in mD; Indicates effective stress, in MPa; The value represents the crack mesh aperture, in meters (m).

[0055] The beneficial effects of this invention are:

[0056] This invention characterizes matrix-fracture fluid exchange by combining an embedded discrete fracture model with a component model, thus solving the problem that conventional fracturing simulation software struggles to reflect the dynamic changes in supercritical carbon dioxide density, viscosity, and other physical properties during fracturing. This invention enables more accurate integrated fracturing-production simulation and provides technical support for tracking the migration path of carbon dioxide throughout the fracturing-well shut-in-production cycle and evaluating the feasibility of post-fracturing carbon dioxide sequestration. Attached Figure Description

[0057] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0058] Figure 1 This is a schematic diagram of the integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to the present invention;

[0059] Figure 2 This is a schematic diagram of the crack mesh and the matrix mesh after refining the mesh near the crack in a specific embodiment;

[0060] Figure 3 Here are the fracture aperture and pressure distribution diagrams after the fracturing simulation in a specific embodiment; where (a) is the fracture aperture distribution diagram; and (b) is the fracture pressure distribution diagram.

[0061] Figure 4 Here are the matrix pressure and carbon dioxide composition distribution diagrams after the fracturing simulation in a specific embodiment; where (a) is the matrix pressure distribution diagram; and (b) is the matrix carbon dioxide composition distribution diagram.

[0062] Figure 5The diagram shows the distribution of matrix pressure and carbon dioxide composition after 10 days of well sealing in a specific embodiment; where (a) is the matrix pressure distribution diagram; and (b) is the matrix carbon dioxide composition distribution diagram.

[0063] Figure 6 The diagram shows the distribution of matrix pressure and carbon dioxide composition after 3 years of constant pressure production in a specific embodiment; where (a) is the matrix pressure distribution diagram and (b) is the matrix carbon dioxide composition distribution diagram. Detailed Implementation

[0064] The present invention will be further described below with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specified, the embodiments and technical features described in this application can be combined with each other. It should also be pointed out that, unless otherwise indicated, all technical and scientific terms used in this application have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains. The terms "comprising" or "including" and similar words used in this invention refer to elements or objects preceding the word that encompass the elements or objects listed following the word and their equivalents, without excluding other elements or objects.

[0065] like Figure 1 As shown, this invention provides an integrated simulation method for shale gas supercritical carbon dioxide fracturing production, comprising the following steps:

[0066] S1: Based on natural cracks and in-situ stress, set potential crack paths, generate crack meshes from the potential crack paths, and assign values ​​to the crack meshes.

[0067] It should be noted that, since the embedded discrete crack model will be used subsequently, the crack mesh generated must conform to the format required by the embedded discrete crack model. This is existing technology, and the specific format will not be elaborated here. Furthermore, the embedded discrete crack model is a crack characterization method that can realistically reflect the actual geometric morphology of cracks. It is also existing technology, and the specific model will not be elaborated here.

[0068] In one specific embodiment, when assigning values ​​to the fracture mesh, an initial fracture mesh pressure and an initial fracture aperture are given, the initial permeability of the fracture mesh is set to be equal to the matrix permeability, and all fracture meshes are initially saturated with only carbon dioxide.

[0069] In one specific embodiment, after generating the crack mesh, the step of refining the matrix mesh containing the crack mesh is further included.

[0070] S2: The planar three-dimensional displacement discontinuity method is used to simulate crack deformation and crack propagation. The embedded discrete crack model and component model are combined to characterize the matrix-crack fluid exchange. Supercritical carbon dioxide fracturing simulation is carried out through iterative coupling.

[0071] In a specific embodiment, when using the planar three-dimensional displacement discontinuity method to simulate crack deformation, the expression for crack deformation is:

[0072] (1)

[0073] In the formula: and Let m represent the tangential displacement discontinuities parallel to and orthogonal to the crack surface. This represents the discontinuity of normal displacement perpendicular to the crack surface, which is equal to the crack mesh aperture. m; superscript , Indicates the crack grid number; and These represent the in-situ stress components acting tangentially to and orthogonal to the crack surface, in MPa; This represents the in-situ stress component acting in the direction perpendicular to the crack surface, expressed in MPa. This represents the crack mesh pressure, in MPa. Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress in the same direction of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress of crack mesh i in the orthogonal direction, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress of crack mesh i in the orthogonal direction, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress in the same direction of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the tangential stress of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the tangential stress in the other direction of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Indicates the total number of open crack meshes;

[0074] When using the planar three-dimensional displacement discontinuity method to simulate crack propagation, the expression for judging crack propagation is:

[0075] (2)

[0076] In the formula: This represents the Type I stress concentration factor of the crack mesh. ; This represents the Young's modulus of the rock, in MPa. Indicates the Poisson's ratio of the rock; The crack mesh half-length is represented in meters (m). Indicates the type I fracture toughness of rock. .

[0077] In one specific embodiment, the component model can characterize the multi-component competitive adsorption mechanism and porous media diffusion mechanism in shale gas;

[0078] The multi-component competitive adsorption mechanism is as follows:

[0079] (3)

[0080] In the formula: This represents the adsorption capacity of component i, in kg / m³. 3 ; This represents the Langmuir adsorption capacity of component i, in kg / m³. 3 ; Represents the mole fraction of component i; Indicates matrix pore pressure, in MPa; The Langmuir adsorption pressure of component i is expressed in MPa. Indicates the total number of components;

[0081] The diffusion mechanism of the porous medium is as follows:

[0082] (4)

[0083] In the formula: This represents the diffusion flux of component i in phase α, in kg / (m²). 2 ·s); Indicates the porosity of the matrix rock; Indicates the saturation degree of the α phase; Indicates the tortuosity of the matrix rock; m represents the diffusion coefficient of component i in phase α. 2 / s; Represents the gradient operator; The density of the α phase is expressed in kg / m³. 3 ; This represents the mole fraction of component i in the α phase.

[0084] In a specific embodiment, the simulation of supercritical carbon dioxide fracturing through iterative coupling includes the following sub-steps:

[0085] S21: Based on the initial crack mesh pressure or the crack mesh pressure calculated in the previous time, give an initial value of crack mesh pressure, and use the planar three-dimensional displacement discontinuity method to update the opening degree of the opened crack mesh based on this initial pressure value.

[0086] S22: The crack mesh aperture updated in step S21 is used to calculate the equivalent carbon dioxide injection time Δt using the following formula. f :

[0087] (7)

[0088] In the formula: The volume of carbon dioxide injected into the ground per unit time is expressed in m. 3 / d; This represents the ground density of carbon dioxide, in kg / m³. 3 ; This represents the crack mesh aperture obtained from step S21, in meters (m). Indicates the initial crack mesh aperture or the crack mesh aperture updated in the previous time step, m; This represents the initial volume of the crack mesh or the volume of the crack mesh updated in the previous time step, m. 3 ; This represents the initial carbon dioxide density in the fractured mesh or the carbon dioxide density after the previous time step update, in kg / m³. 3 ;

[0089] S23: Import the fracture mesh aperture obtained in step S21 and the fracture mesh permeability calculated therefrom into the reservoir simulator, and use the component model to simulate and obtain dt+Δt. f Crack mesh pressure after time;

[0090] S24: Determine whether the crack pressure converges based on the initial crack mesh pressure given in step S21 and the crack mesh pressure obtained by simulation in step S23.

[0091] If convergence occurs, the crack pressure cycle is complete, and the process proceeds to step S25.

[0092] If convergence is not achieved, update the initial value of the crack mesh pressure and repeat steps S21-S24.

[0093] S25: Determine whether the crack tip is extending forward;

[0094] If so, then activate the corresponding crack mesh and repeat steps S21-S25;

[0095] If not, the crack propagation cycle is completed, and the process proceeds to step S26;

[0096] S26: Determine whether the current fracturing simulation time has reached the fracturing end time;

[0097] If so, then complete the fracturing simulation;

[0098] If not, repeat steps S21-S26.

[0099] In a specific embodiment, in step S21, the initial value of the crack mesh pressure is given by the following formula:

[0100] (5)

[0101] In the formula: This represents the initial pressure value of the given crack mesh, in MPa; This represents the initial fracture mesh pressure or the fracture mesh pressure calculated at the previous time step, in MPa; Indicates the parameters of crack pressure disturbance;

[0102] In step S24, the initial value of the crack mesh pressure is updated using the following formula:

[0103] (6)

[0104] In the formula: This indicates the updated crack mesh pressure; This represents the iteration parameter, whose value is between 0 and 1; This represents the crack mesh pressure obtained from step S23 of the simulation, in MPa.

[0105] In a specific embodiment, in step S23, the permeability of the fracture mesh is calculated using the following formula:

[0106] (8)

[0107] In the formula: The permeability of the crack mesh is expressed in mD. The crack mesh aperture is represented in meters (m).

[0108] In step S24, the expression for determining whether the crack pressure has converged is:

[0109] (9)

[0110] In the formula: The convergence tolerance is set by the user.

[0111] S3: Based on the fracture opening and the permeability calculated from it after the fracturing simulation, production simulation is carried out using the pressure field and composition distribution after the fracturing simulation as initial conditions.

[0112] In one specific embodiment, the permeability is calculated using the following formula:

[0113] (10)

[0114] In the formula: Permeability, expressed in mD; Indicates effective stress, in MPa; The value represents the crack mesh aperture, in meters (m).

[0115] In a specific embodiment, taking a target shale gas reservoir as an example, the integrated simulation method for supercritical carbon dioxide fracturing production of shale gas described in this invention is used to simulate its fracturing production process. In this embodiment, the reservoir simulator containing an embedded discrete model module is UNCONG, and the established fracture mesh and the matrix mesh containing the fracture mesh after refinement are shown below. Figure 2 As shown. The specific method for densification is as follows: Identify matrix meshes whose vertical distance from the center of the matrix mesh to the crack surface is less than the densification range; densify the identified matrix meshes according to the densification level. The densification level is the number of subdivisions of the matrix mesh in the x and y directions.

[0116] In this embodiment, the total duration of the fracturing simulation is 5 minutes, and the fracture aperture distribution obtained after the fracturing simulation is as follows: Figure 3 As shown, the pressure and composition changes in the near-fracture region caused by the filtration of carbon dioxide into the matrix are as follows: Figure 4 As shown.

[0117] The fracture geometry parameters obtained from the fracturing simulation were used for subsequent production simulation, and the pressure field and composition distribution from the post-fracturing reservoir numerical simulation were used as the initial conditions for the production simulation. After shutting down the injection well and simulating a 10-day well suffocation period, the matrix pressure and carbon dioxide composition distribution were as follows: Figure 5 As shown. With a bottomhole flowing pressure of 20 MPa, the matrix pressure and carbon dioxide composition distribution after 3 years of well production are as follows. Figure 6 As shown, after three years of production, the carbon dioxide injected into the matrix during the fracturing stage is almost entirely produced from the production well.

[0118] In summary, this invention employs a planar three-dimensional displacement discontinuity method to describe fracture deformation and propagation, and combines an embedded discrete fracture model and a component model to characterize matrix-fracture fluid exchange. This allows for the description of matrix pressure and composition changes in the near-fracture region during fracturing, enabling more accurate integrated simulation of supercritical carbon dioxide fracturing production. Compared to existing technologies, this invention represents a significant advancement.

[0119] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. An integrated simulation method for shale gas supercritical carbon dioxide fracturing production, characterized in that, Includes the following steps: S1: Based on natural cracks and in-situ stress, set potential crack paths, generate crack meshes from the potential crack paths, and assign values ​​to the crack meshes; S2: The planar three-dimensional displacement discontinuity method is used to simulate crack deformation and crack propagation. The embedded discrete crack model and component model are combined to characterize the matrix-crack fluid exchange. Supercritical carbon dioxide fracturing simulation is carried out through iterative coupling. When simulating crack deformation using the planar three-dimensional displacement discontinuity method, the expression for crack deformation is: (1) In the formula: and Let m represent the tangential displacement discontinuities parallel to and orthogonal to the crack surface. This represents the discontinuity of normal displacement perpendicular to the crack surface, which is equal to the crack mesh aperture. m; superscript , Indicates the crack grid number; and These represent the in-situ stress components acting tangentially to and orthogonal to the crack surface, in MPa; This represents the in-situ stress component acting in the direction perpendicular to the crack surface, expressed in MPa. This represents the crack mesh pressure, in MPa. Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress in the same direction of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress of crack mesh i in the orthogonal direction, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress of crack mesh i in the orthogonal direction, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the tangential stress in the same direction of crack mesh i, in MPa / m; Let be the influence coefficient of the tangential displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the tangential stress of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the tangential stress in the other direction of crack mesh i, in MPa / m; Let be the influence coefficient of the normal displacement discontinuity of crack mesh j on the normal stress of crack mesh i, in MPa / m; Indicates the total number of open crack meshes; When using the planar three-dimensional displacement discontinuity method to simulate crack propagation, the expression for judging crack propagation is: (2) In the formula: This represents the Type I stress concentration factor of the crack mesh. ; This represents the Young's modulus of the rock, in MPa. Indicates the Poisson's ratio of the rock; The crack mesh half-length is represented in meters (m). Indicates the type I fracture toughness of rock. ; The component model can characterize the multi-component competitive adsorption mechanism and porous media diffusion mechanism in shale gas. The multi-component competitive adsorption mechanism is as follows: (3) In the formula: This represents the adsorption capacity of component i, in kg / m³. 3 ; This represents the Langmuir adsorption capacity of component i, in kg / m³. 3 ; Represents the mole fraction of component i; Indicates matrix pore pressure, in MPa; The Langmuir adsorption pressure of component i is expressed in MPa. Indicates the total number of components; The diffusion mechanism of the porous medium is as follows: (4) In the formula: This represents the diffusion flux of component i in phase α, in kg / (m²). 2 ·s); Indicates the porosity of the matrix rock; Indicates the saturation degree of the α phase; Indicates the tortuosity of the matrix rock; m represents the diffusion coefficient of component i in phase α. 2 / s; Represents the gradient operator; The density of the α phase is expressed in kg / m³. 3 ; This represents the mole fraction of component i in phase α; S3: Based on the fracture opening and the permeability calculated from it after the fracturing simulation, production simulation is carried out using the pressure field and composition distribution after the fracturing simulation as initial conditions.

2. The integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to claim 1, characterized in that, In step S1, when assigning values ​​to the fracture mesh, the initial fracture mesh pressure and initial fracture aperture are given, the initial permeability of the fracture mesh is made equal to the matrix permeability, and all fracture meshes are saturated with carbon dioxide at the initial moment.

3. The integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to claim 1, characterized in that, In step S1, after generating the crack mesh, the step further includes refining the matrix mesh containing the crack mesh.

4. The integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to claim 1, characterized in that, Step S2, the simulation of supercritical carbon dioxide fracturing through iterative coupling, specifically includes the following sub-steps: S21: Based on the initial crack mesh pressure or the crack mesh pressure calculated in the previous time, give an initial value of crack mesh pressure, and use the planar three-dimensional displacement discontinuity method to update the opening degree of the opened crack mesh based on this initial pressure value. S22: The equivalent carbon dioxide injection time Δt is obtained by calculating the crack mesh aperture updated in step S21. f ; S23: Import the fracture mesh aperture obtained in step S21 and the fracture mesh permeability calculated therefrom into the reservoir simulator, and use the component model to simulate and obtain dt+Δt. f Crack mesh pressure after time; S24: Determine whether the crack pressure converges based on the initial crack mesh pressure given in step S21 and the crack mesh pressure obtained by simulation in step S23. If convergence occurs, the crack pressure cycle is complete, and the process proceeds to step S25. If convergence is not achieved, update the initial value of the crack mesh pressure and repeat steps S21-S24. S25: Determine whether the crack tip is extending forward; If so, then activate the corresponding crack mesh and repeat steps S21-S25; If not, the crack propagation cycle is completed, and the process proceeds to step S26; S26: Determine whether the current fracturing simulation time has reached the fracturing end time; If so, then complete the fracturing simulation; If not, repeat steps S21-S26.

5. The integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to claim 4, characterized in that, In step S21, the initial value of the crack mesh pressure is given by the following formula: (5) In the formula: This represents the initial pressure value of the given crack mesh, in MPa; This represents the initial fracture mesh pressure or the fracture mesh pressure calculated at the previous time step, in MPa; Indicates the parameters of crack pressure disturbance; In step S24, the initial value of the crack mesh pressure is updated using the following formula: (6) In the formula: This indicates the updated crack mesh pressure; This represents the iteration parameter, whose value is between 0 and 1; This represents the crack mesh pressure obtained from step S23 of the simulation, in MPa.

6. The integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to claim 4, characterized in that, In step S22, the equivalent carbon dioxide injection time Δt f The calculation is performed using the following formula: (7) In the formula: The volume of carbon dioxide injected into the ground per unit time is expressed in m. 3 / d; This represents the ground density of carbon dioxide, in kg / m³. 3 ; This represents the crack mesh aperture obtained from step S21, in meters (m). Indicates the initial crack mesh aperture or the crack mesh aperture updated in the previous time step, m; This represents the initial volume of the crack mesh or the volume of the crack mesh updated in the previous time step, m. 3 ; This represents the initial carbon dioxide density in the fractured mesh or the carbon dioxide density after the previous time step update, in kg / m³. 3 .

7. The integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to claim 4, characterized in that, In step S23, the permeability of the fracture mesh is calculated using the following formula: (8) In the formula: The permeability of the crack mesh is expressed in mD. The crack mesh aperture is represented in meters (m). In step S24, the expression for determining whether the crack pressure has converged is: (9) In the formula: This represents the crack mesh pressure obtained from step S23 of the simulation, in MPa. This represents the initial crack mesh pressure value given in step S21, in MPa; The convergence tolerance is set by the user.

8. The integrated simulation method for shale gas supercritical carbon dioxide fracturing production according to any one of claims 1-7, characterized in that, In step S3, the permeability is calculated using the following formula: (10) In the formula: Permeability, expressed in mD; Indicates effective stress, in MPa; The value represents the crack mesh aperture, in meters (m).

Images