Fracturing simulation method and system for coupling wellbore-fracture-cave of fracture-cave medium

Abstract

The application discloses a kind of joint hole medium wellbore-fracture-cave coupling fracturing simulation method and system, method includes: constructing simulation model and discretization obtains rock matrix, fracture, wellbore and cave unit;Based on the solution of thermal-fluid-solid coupling mathematical model obtains main variable field and stress boundary condition;Based on the solution of induced stress calculation model obtains wellbore-fracture-cave joint induced stress field and fracture physical property parameter;Through iterative convergence judgment updates thermal-fluid-solid coupling parameter;According to joint induced stress field, fracture propagation criterion is judged and fracture geometry is updated;The application realizes the collaborative simulation of thermal-fluid-solid multi-field coupling and fracture dynamic expansion, can accurately track the expansion track of three-dimensional fracture in joint hole medium, provides reliable prediction means for joint hole oil and gas reservoir fracturing, geothermal development and underground engineering stability analysis.

Description

Technical Field

[0001] This invention relates to the field of oil and gas reservoir development, and in particular to a coupled fracturing simulation method and system for fracture-cavity type wellbore-fracture-cavity. Background Technology

[0002] The initiation and propagation of cracks are ubiquitous physical phenomena in nature and engineering science, widely observed in dike-induced earthquakes, meltwater fractures in ice shelves, dam failures, fatigue damage in aerospace engineering, and hydraulic fracturing in petroleum, geothermal, and mining engineering. These crack problems typically involve strong coupling between mechanical, chemical, thermodynamic, and hydraulic processes. The multiphysics coupling nature of these problems presents significant challenges to understanding rock fracture mechanisms and developing effective engineering guidelines, making highly scalable and robust numerical models indispensable tools for addressing crack problems.

[0003] Traditional finite element methods (FEMs) struggle to handle dynamic crack propagation, requiring frequent mesh remapping. Extended finite element methods (XFEM / GFEMs) eliminate the need for remapping by enriching basis functions, but face challenges such as complex stiffness matrix assembly, difficulties in multiphysics coupling, and low efficiency in tracking multiple cracks. Furthermore, most models are limited to two dimensions or neglect fluid flow within the cavity. While phase-field methods (PFMs) can avoid explicit crack tracking, their representation of cracks as diffuse regions results in high computational costs, making it difficult to balance accuracy and efficiency.

[0004] The displacement discontinuity method (DDM) in the boundary element method has unique advantages in crack simulation because it only requires boundary discretization. However, the classical DDM formula is not accurate enough when dealing with karst caves, so the virtual stress method (FSM) needs to be introduced to accurately describe the stress transitions caused by karst caves. However, existing integrated models of DDM and FSM are mainly used to analyze stress distribution in tunnels or slopes, and there are no reports on simulating the interaction between three-dimensional non-planar cracks and natural karst caves or well shafts in thermally porous elastic media.

[0005] In summary, the current full three-dimensional thermal-fluid-structure interaction numerical simulation technology for crack propagation in fractured-vuggy composite media is still immature and has the following main shortcomings:

[0006] (1) Difficulty in characterizing non-uniform stress field: In fractured media, caverns, wells and cracks coexist. Geometric discontinuities and material abrupt changes make it difficult for conventional models to accurately describe the stress disturbances induced by their combination.

[0007] (2) The evolution law of three-dimensional non-coplanar cracks is unclear: Due to the influence of multi-level crack and cavity structure, the crack propagation path is complex, and the existing quasi-static model is difficult to capture the deflection, branching and connection criteria of cracks when they approach the karst cavity. Summary of the Invention

[0008] In view of the shortcomings of the prior art, the present invention provides a coupled fracturing simulation method and system for fractured-cavitary media wellbore-fracture-cavity, which solves the problem of difficulty in accurately characterizing the complex stress interference and the distortion of three-dimensional fracture evolution trajectory prediction under the combined action of wellbore, fracture and cavity in fractured-cavitary composite media in the prior art.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] In a first aspect, the present invention provides a coupled fracturing simulation method for fractured-cavitary media wellbore-fracture-cavitary coupling, comprising the following steps:

[0011] A. Construct a coupled fracturing simulation model of wellbore-fracture-cavity, and discretize the study area into a grid to obtain the volume elements, fracture elements, wellbore elements and cavity elements of the rock matrix;

[0012] The coupled fracturing simulation model includes a thermo-fluid-structure interaction mathematical model and an induced stress calculation model;

[0013] B. Solve based on the aforementioned thermal-fluid-structure interaction mathematical model and the current geometry of the crack to obtain the principal variable field, and extract stress boundary conditions from the principal variable field;

[0014] The main variable fields include the rock mass deformation field, pressure field, and temperature field;

[0015] The stress boundary conditions include fluid pressure inside the fracture, fluid pressure inside the cavern, and fluid pressure inside the wellbore.

[0016] C. Solve the stress-induced stress calculation model and stress boundary conditions to obtain the combined induced stress field of the wellbore-fracture-cavity and the fracture physical parameters.

[0017] The fracture physical properties include fracture aperture, fracture permeability, and fracture stiffness;

[0018] D. During the iterative solution process at the current time step, compare the main variable field and crack opening in the current iteration step with those in the previous iteration step, and determine whether the preset convergence condition is met.

[0019] If the conditions are not met, then based on the crack physical properties, update the parameters in the thermal-fluid-structure interaction mathematical model and return to step B; if the conditions are met, then proceed to step E.

[0020] E. Using the current crack geometry as the input for the next time step, return to step B until the preset simulation duration is reached, and obtain the crack propagation trajectory and the corresponding master variable field.

[0021] In an optional implementation, the method further includes the following steps after step D and before step E:

[0022] Based on the joint induced stress field, it is determined whether the crack meets the propagation criterion, and whether the current crack geometry needs to be updated.

[0023] If satisfied, determine the propagation direction and amount of the crack front edge, and update the current crack geometry.

[0024] If the conditions are not met, the current crack geometry will remain unchanged.

[0025] In an optional implementation, the thermal-fluid-structure interaction mathematical model includes the momentum conservation equation, the heat and mass transfer integral equation of the rock matrix, the heat and mass transfer integral equation of the fracture, and the mass conservation equation of the karst cave.

[0026] The integral form of the momentum conservation equation is:

[0027] ;

[0028] ;

[0029] In the formula, A is the three-dimensional integral region; A is the two-dimensional surface integral region. Let i be the component of the concentrated force along the i-th direction; Let i be the component of the surface force along the i-th direction; , and These are the areas affected by surface forces, the surface of the cave, and the surface of the cracks, respectively. This is the contact stress tensor at the crack wall; , and These are the fluid pressures within the rock matrix, fissures, and caverns, respectively. The external normal vector of the cave surface; is the outward normal vector of the crack surface; , is the component of the outward normal vector of the surface of the cave or crack along the i-direction; For the volumetric forces of the rock mass matrix; The volume density of the rock matrix; The component of gravitational acceleration along the i-direction; , and These represent the displacements along the j, k, and l directions, respectively. Second-order unit tensor; For the displacement elastic tensor; The coefficient of thermal expansion; It is the bulk modulus of the drainage. This refers to the change in temperature. Let k be the component of the rock mass displacement vector in the k direction; Let be the component of the rock mass displacement vector in the l direction.

[0030] In an optional embodiment, the integral equation for heat and mass transfer of the rock matrix is ​​expressed as follows:

[0031] ;

[0032] ;

[0033] ;

[0034] ;

[0035] In the formula, For temperature; Entropy; The coefficient of thermal expansion of the rock matrix; This represents the temperature change within the rock matrix. This represents the temperature change within the crack. The equivalent volumetric heat capacity of the rock mass; For the volumetric strain of the rock mass; Biot modulus; H is the bulk modulus of the rock mass; H is the specific enthalpy of the fluid. This is the rock matrix permeability tensor; Let j be the component of gravitational acceleration along the j-direction; For fluid viscosity; is the component of the outward normal vector of the surface of the cave or crack along the j direction; These are the source and sink terms of mass in the rock mass matrix; It represents the source and sink terms of enthalpy within the rock matrix. Let the thermal conductivity tensor be... This is for simulating time.

[0036] The integral equation for heat and mass transfer in the crack is expressed as follows:

[0037] ;

[0038] ;

[0039] ;

[0040] ;

[0041] In the formula, and These are crack stiffness and fluid stiffness, respectively. The crack opening; The coefficient of thermal expansion of the fluid; Specific heat capacity of the fluid; The temperature inside the crack; Let f be the fracture permeability tensor; It is the enthalpy of the fluid; and These are the source and sink terms for the mass and enthalpy within the crack, respectively. This is the region of crack integration; Divide the contact area between the cave and the fissure into zones; The crack-matrix contact area is divided into regions;

[0042] The mass conservation equation for the karst cave is expressed as follows:

[0043] ;

[0044] In the formula, This is the integral region of the karst cave; Let i be the modulus of the i-th karst cave; It is a source and sink of mass within the karst cave.

[0045] In an optional implementation, step B includes:

[0046] Based on the aforementioned volume element, cavern element, and crack element, the momentum conservation equation in the thermal-fluid-structure interaction mathematical model is spatially discretized using the finite element method, the mass conservation equation and energy conservation equation are spatially discretized using the finite volume method, and time discretized using the finite difference method to obtain a set of nonlinear algebraic equations.

[0047] The nonlinear algebraic equations were linearized using the Newton-Raphson method, and the Jacobian matrix was assembled to solve the linearized algebraic equations, thereby obtaining the rock mass deformation field, pressure field, and temperature field.

[0048] The rock mass deformation field includes the displacement of the rock mass matrix;

[0049] The pressure field includes rock matrix pressure, fluid pressure within fractures, and fluid pressure within karst caves.

[0050] The temperature field includes the temperature distribution within the rock matrix and the temperature distribution within the fractures;

[0051] Obtain the preset fluid pressure inside the wellbore, and use the fluid pressure inside the fracture, the fluid pressure inside the cavern, and the fluid pressure inside the wellbore as the stress boundary conditions for the induced stress calculation model.

[0052] In an optional implementation, step C includes:

[0053] Based on the induced stress calculation model, the displacement discontinuity method is used to perform surface integration on the crack element through the Kelvin fundamental solution to construct the influence coefficient relationship between the displacement discontinuity of the crack element and the induced stress generated in space.

[0054] Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the wellbore unit and the induced stress generated in space.

[0055] Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the karst unit and the induced stress generated in space.

[0056] The displacement discontinuities of all fracture elements, the virtual stresses of all wellbore elements, and the virtual stresses of all karst elements are taken as the basic unknowns.

[0057] The fluid pressure in the fracture, the fluid pressure in the cave, and the fluid pressure in the wellbore are applied to the corresponding fracture unit, cave unit, and wellbore unit to establish a joint linear equation system.

[0058] Solving the combined linear equations yields the displacement discontinuities of each fracture element, the virtual stress of each wellbore element, and the virtual stress of each karst cave element.

[0059] Based on the displacement discontinuities and virtual stresses obtained from the solution, the induced stresses generated in space by the cracks, well shafts and karst caves are calculated respectively through the influence coefficient relationship;

[0060] By linearly superimposing all induced stresses, a combined induced stress field of wellbore-fracture-cavity is obtained.

[0061] In an optional implementation, step C further includes:

[0062] Based on the induced stress field of the wellbore-fracture-cavity system, the displacement discontinuity in the normal direction of the fracture unit is obtained to determine the fracture aperture.

[0063] Based on the law of cubes, the fracture permeability is calculated according to the fracture aperture.

[0064] The crack stiffness is calculated based on the crack opening and the fluid pressure inside the crack.

[0065] In an optional implementation, the induced stress field of the karst cave is expressed as a function:

[0066] ;

[0067] In the formula, , and Each cave unit is a separate unit for all cave units. Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; For cave units Virtual stress at the direction of the cave unit The influence coefficient of shear-induced stress generated at the location; For cave units Virtual stress at the direction of the cave unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units Virtual stress at the direction of the cave unit The influence coefficient of the normal induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of shear-induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of the normal induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of shear-induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of the normal induced stress generated at the location; and Separately, it is a cave unit. Virtual stress along the dip, strike, and normal directions; This represents the total number of all cave units.

[0068] The induced stress field of the crack is expressed as a function:

[0069] ;

[0070] In the formula, , and For all crack elements, in the crack element Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of the normal induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of the normal induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of the normal induced stress generated at the location; and Do not use crack units The discontinuous displacement along the dip, strike, and normal directions; The total number of all crack elements. It is a collection of closed cracks.

[0071] Secondly, this invention provides a coupled fracturing simulation system for fractured-cavity media wellbore-fracture-cavity, comprising:

[0072] The model building module is used to construct a coupled fracturing simulation model of wellbore-fracture-cavity, and to discretize the study area into a grid to obtain the volume elements, fracture elements, wellbore elements and cavity elements of the rock matrix;

[0073] The coupled fracturing simulation model includes a thermo-fluid-structure interaction mathematical model and an induced stress calculation model;

[0074] The first solution module is used to solve the problem based on the thermal-fluid-structure interaction mathematical model and the current geometry of the crack to obtain the main variable field and extract the stress boundary conditions from the main variable field.

[0075] The main variable fields include the rock mass deformation field, pressure field, and temperature field;

[0076] The stress boundary conditions include fluid pressure inside the fracture, fluid pressure inside the cavern, and fluid pressure inside the wellbore.

[0077] The second solution module is used to solve the stress field and fracture physical parameters of the wellbore-fracture-cavity joint induced stress field based on the induced stress calculation model and the stress boundary conditions.

[0078] The fracture physical properties include fracture aperture, fracture permeability, and fracture stiffness;

[0079] The coupled iterative control module is used to compare the main variable field and crack aperture in the current iteration step and the previous iteration step during the iterative solution process in the current time step, and to determine whether the preset convergence condition is met.

[0080] If the conditions are not met, the parameters in the thermal-fluid-structure interaction mathematical model are updated based on the crack physical property parameters, and the process returns to the first solution module; if the conditions are met, the trajectory acquisition module is triggered.

[0081] The trajectory acquisition module is used to take the current crack geometry as the input for the next time step, return to the first solution module, and continue until the preset simulation duration is reached, and obtain the crack propagation trajectory and the corresponding master variable field.

[0082] In an optional implementation, after the coupling iteration control module and before the trajectory acquisition module, the following further component is included:

[0083] The crack propagation update module is used to determine whether the crack meets the propagation criterion based on the joint induced stress field.

[0084] If the propagation criterion is met, the propagation direction and propagation amount of the crack front are determined according to the joint induced stress field, and the current crack geometry is updated.

[0085] If the expansion criterion is not met, the current crack geometry remains unchanged.

[0086] The beneficial effects of the embodiments provided by the present invention include:

[0087] This invention integrates the displacement discontinuity method and the virtual stress method to construct a joint solution framework that can simultaneously handle the discontinuous characteristics of fracture displacement and the discontinuous characteristics of wellbore and karst cave stress. It can accurately characterize the joint induced stress field caused by the wellbore-fracture-karst cave system, and solve the problem of stress prediction distortion caused by neglecting the mechanical interaction of the three in the existing technology. This lays a reliable stress foundation for fracture propagation simulation.

[0088] This invention employs an iterative update mechanism for the main variable field and fracture physical parameters. Through convergence judgment and time step advancement, it realizes the coupled calculation of thermal-fluid-solid multi-field evolution and dynamic fracture propagation. It can accurately track the propagation trajectory of three-dimensional fractures in fracture-vuggy media, significantly improving the accuracy and reliability of oil and gas fracturing process simulation. Attached Figure Description

[0089] To more clearly illustrate the technical solutions in the embodiments of this application 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 this application. For those skilled in the art, other embodiments can be obtained based on these drawings.

[0090] Figure 1 This is a schematic flowchart of a coupled fracturing simulation method for a fractured-cavity media wellbore-fracture-cavity provided in an embodiment of the present invention;

[0091] Figure 2 This is a schematic diagram illustrating the principle of solving a karst cave using the virtual stress method provided in this embodiment of the invention.

[0092] Figure 3 This is a schematic diagram illustrating the principle of solving cracks using the displacement discontinuity method provided in this embodiment of the invention;

[0093] Figure 4 This is a schematic diagram of the numerical discretization of the far-field in-situ stress in a vertical wellbore provided in an embodiment of the present invention;

[0094] Figure 5 This is a comparison diagram of the numerical results and analytical solutions of the disturbance stress field of a cylindrical wellbore provided in the embodiments of the present invention; wherein, Figure 5 (a) is a schematic diagram showing the variation of effective stress with radial distance from the wellbore; Figure 5 (b) is a schematic diagram of the effective stress distribution at the well wall in the π plane;

[0095] Figure 6 This is a comparison chart of the numerical results and analytical solutions for the pressure-bearing spherical karst cave provided in the embodiments of the present invention; wherein, Figure 6 (a) is a schematic diagram of a triangular grid in a spherical karst cave; Figure 6 (b) is a historical curve of radial induced displacement of the spherical karst cave; Figure 6 (c) is a graph showing the historical variation of tangential induced stress in a spherical karst cave; Figure 6 (d) is a graph showing the historical variation of radial induced stress in a spherical karst cave;

[0096] Figure 7 This is a flowchart of a coupled fracturing simulation system for a fractured-cavity media wellbore-fracture-cavity provided in an embodiment of the present invention.

[0097] Figure 8 This is a schematic diagram of a numerical model for fracture propagation in fractured-vuggy reservoirs provided in an embodiment of the present invention; wherein, Figure 8 (a) is a schematic diagram of the full three-dimensional numerical computation grid; Figure 8 (b) is a schematic diagram of the grid subdivision results of the cave-fracture-wellbore;

[0098] Figure 9 These are illustrations showing the results of competing cracks after different propagation steps, provided in an embodiment of the present invention; wherein, Figure 9 (a) is a schematic diagram of the results after 5 steps of expansion; Figure 9 (b) is a schematic diagram of the results after 10 steps of expansion; Figure 9 (c) is a schematic diagram of the results after 15 steps of expansion; Figure 9 (d) is a schematic diagram of the results after 20 steps of expansion;

[0099] Figure 10 This is a temperature and pressure distribution cloud map of a fractured-vuggy reservoir after fracturing, provided in an embodiment of the present invention; wherein, Figure 10 (a) is a temperature field contour map; Figure 10 (b) is a pressure field contour map. Detailed Implementation

[0100] The features and exemplary embodiments of various aspects of the present invention will now be described in detail. In the following detailed description, numerous specific details are set forth in order to provide a thorough understanding of the invention. However, it will be apparent to those skilled in the art that the invention may be practiced without requiring some of these specific details. The following description of embodiments is merely intended to provide a better understanding of the invention by illustrating examples of the invention.

[0101] In some of the processes described in the specification, claims, and accompanying drawings of this invention, multiple operations appear in a specific order. However, it should be clearly understood that these operations may not be performed in the order they appear herein, or they may be performed in parallel. The operation numbers, such as S1, S2, etc., are merely used to distinguish different operations and do not themselves represent any execution order. Furthermore, these processes may include more or fewer operations, and these operations may be performed sequentially or in parallel.

[0102] Example 1

[0103] like Figure 1 As shown, this embodiment provides a coupled fracturing simulation method for fractured-cavitary media wellbore-fracture-cavity, including the following steps:

[0104] S1. Construct a coupled fracturing simulation model of wellbore-fracture-cavity, and discretize the study area into a grid to obtain the volume elements, fracture elements, wellbore elements and cavity elements of the rock matrix;

[0105] The coupled fracturing simulation model includes a thermo-fluid-structure interaction mathematical model and an induced stress calculation model. The thermo-fluid-structure interaction mathematical model and the induced stress calculation model together constitute the hybrid simulation framework of this method, realizing the collaborative solution of multi-field coupling process and three-dimensional evolution of fracture in fracture-vuggy media.

[0106] In this embodiment, the thermal-fluid-structure interaction mathematical model includes the momentum conservation equation, the heat and mass transfer integral equation of the rock matrix, the heat and mass transfer integral equation of the fracture, and the mass conservation equation of the cavern, which are used to describe the coupled interaction between the stress field, seepage field and temperature field in the fracture-cavity medium. The thermal-fluid-structure interaction mathematical model stipulates the sign rule of tensile stress as positive and compressive stress as negative, which is consistent with the general calculation standard in the field of rock mechanics.

[0107] In some embodiments, to predict the propagation path of a three-dimensional crack, the study area needs to be discretized into a mesh, specifically including:

[0108] The study area is divided into multiple hexahedral units to discretize the rock matrix, resulting in a set of volumetric units; wherein, the study area is the extent of the rock mass under study;

[0109] The surfaces of fractures, caves, and wells embedded in the rock mass are discretized into multiple triangular units to obtain sets of fracture units, sets of cave units, and sets of wells units, respectively.

[0110] Specifically, the integral form of the momentum conservation equation is:

[0111] ;

[0112] ;

[0113] In the formula, A is the three-dimensional integral region; A is the two-dimensional surface integral region. Let i be the component of the concentrated force along the i-th direction; Let i be the component of the surface force along the i-th direction; , and These are the areas affected by surface forces, the surface of the cave, and the surface of the cracks, respectively. This is the contact stress tensor at the crack wall; , and These are the fluid pressures within the rock matrix, fissures, and caverns, respectively. The external normal vector of the cave surface; is the outward normal vector of the crack surface; , is the component of the outward normal vector of the surface of the cave or crack along the i-direction; For the volumetric forces of the rock mass matrix; The volume density of the rock matrix; The component of gravitational acceleration along the i-direction; , and These represent the displacements along the j, k, and l directions, respectively. Second-order unit tensor; For the displacement elastic tensor; The coefficient of thermal expansion; It is the bulk modulus of the drainage. This refers to the change in temperature. Let k be the component of the rock mass displacement vector in the k direction; Let be the component of the rock mass displacement vector in the l direction.

[0114] The functional expression for the volume density of the rock mass matrix is:

[0115] ;

[0116] In the formula, For fluid density; ρ is the density of solid rock particles; b is the Biot coefficient. Rock porosity;

[0117] Specifically, the integral equations for heat and mass transfer in rock matrix include the matrix fluid mass conservation equation and the matrix energy conservation equation.

[0118] The expression for the matrix fluid mass conservation equation is as follows:

[0119] ;

[0120] ;

[0121] The expression for the matrix energy conservation equation is as follows:

[0122] ;

[0123] ;

[0124] In the formula, For temperature; Entropy; The coefficient of thermal expansion of the rock matrix; This represents the temperature change within the rock matrix. This represents the temperature change within the crack. The equivalent volumetric heat capacity of the rock mass; For the volumetric strain of the rock mass; Biot modulus; H is the bulk modulus of the rock mass; H is the specific enthalpy of the fluid. This is the rock matrix permeability tensor; Let j be the component of gravitational acceleration along the j-direction; For fluid viscosity; is the component of the outward normal vector of the surface of the cave or crack along the j direction; These are the source and sink terms of mass in the rock mass matrix; It represents the source and sink terms of enthalpy within the rock matrix. Let the thermal conductivity tensor be... This is for simulating time.

[0125] The thermal expansion coefficient of the rock matrix is ​​expressed as follows:

[0126] ;

[0127] In the formula, The coefficient of thermal expansion of the fluid;

[0128] The equivalent volumetric heat capacity of the rock mass is expressed as follows:

[0129] ;

[0130] In the formula, The specific heat capacity of rock solid at constant pressure; The specific heat capacity at constant pressure of the fluid;

[0131] The integral equations for heat and mass transfer in the fracture include the fracture fluid mass conservation equation and the fracture energy conservation equation, and their functional expressions are as follows:

[0132] ;

[0133] ;

[0134] ;

[0135] ;

[0136] In the formula, and These are crack stiffness and fluid stiffness, respectively. The crack opening; The coefficient of thermal expansion of the fluid; Specific heat capacity of the fluid; The temperature inside the crack; Let f be the fracture permeability tensor; It is the enthalpy of the fluid; and These are the source and sink terms for the mass and enthalpy within the crack, respectively. This is the region of crack integration; Divide the contact area between the cave and the fissure into zones; The crack-matrix contact area is divided into regions;

[0137] In some embodiments, to simplify the description of free flow within a cave, reduce computational load, and ensure engineering accuracy, the cave is approximated as an equipotential body filled with fluid. Numerical simulation results of crevice flow based on the Darcy-Stokes equations show that the pressure distribution inside the cave is uniform with almost no pressure drop, and the equipotential body assumption has sufficient engineering applicability.

[0138] This method neglects the influence of temperature on cave deformation and pressure changes, and constructs a mass conservation equation for the cave based on the equipotential body assumption. The expression for this equation is:

[0139] ;

[0140] In the formula, This is the integral region of the karst cave; The modulus of the karst cave; It is a source and sink of mass within the karst cave.

[0141] In some embodiments, the wellbore in the thermal-fluid-structure interaction mathematical model is not established with separate governing equations. Its function is achieved through the mechanical, flow, and thermal boundary conditions at the well wall, which is fundamentally different from the treatment of karst caves. The fluid pressure inside the wellbore is calculated independently based on the injection construction parameters and is not used as a variable to be solved in the model.

[0142] S2. Solve based on the aforementioned thermal-fluid-structure interaction mathematical model and the current geometry of the crack to obtain the principal variable field, and obtain the stress boundary conditions through the principal variable field;

[0143] The main variable fields include rock deformation field, pressure field and temperature field;

[0144] It should be noted that, in scenarios with sufficient on-site basic data, the initial values ​​are directly assigned from the logging, well logging, and core experiment data of the target reservoir; in scenarios with insufficient on-site basic data, the initial estimated values ​​are obtained from a decoupled rock mechanics and fluid-thermal flow model: in this decoupled model, the porous medium is under hydrostatic loading, the vertical stress is calculated from the self-weight of the overlying strata, and the initial fluid pressure in fractures and caverns is equal to the matrix pore pressure at the corresponding location; given a reference depth and reference temperature, the initial temperature distribution is assigned to the study area according to the conventional geothermal gradient;

[0145] After the initial values ​​are assigned, external loads such as lateral boundary surface forces, volume forces, fluid pressures, and thermal stresses are applied to the rock mechanics model to calculate the initial deformation state of the rock mass.

[0146] In some embodiments, step S2 is used to perform spatiotemporal discretization of the thermal-fluid-structure interaction mathematical model to construct a solvable system of nonlinear algebraic equations, specifically including:

[0147] S201. Spatial discretization is performed on the continuum domain defined by the thermal-fluid-structure interaction mathematical model. A finite-dimensional subspace is introduced to approximate the pressure, temperature, and displacement fields, and a corresponding trial function space is defined.

[0148] Among them, the volume domain is the three-dimensional continuous spatial region where the rock matrix is ​​located;

[0149] Specifically, the function expression for the trial function space is:

[0150] ;

[0151] ;

[0152] ;

[0153] ;

[0154] ;

[0155] ;

[0156] In the formula, , , These are the pressure solution space, temperature solution space, and displacement solution space, respectively. , , These are the pressure trial function space, the temperature trial function space, and the displacement trial function space, respectively. The volume domain is the three-dimensional continuous spatial region where the rock matrix is ​​located; , , These are pressure boundaries, temperature boundaries, and displacement boundaries, respectively. , , These are the pressure, temperature, and displacement values ​​given at the boundary, respectively. It is a space of square-integrable functions; Let d be a first-order Sobolev space; d is the spatial dimension.

[0157] S202. Based on the defined trial function space, construct shape functions to approximate the pressure field, temperature field, and displacement field, as follows:

[0158] ;

[0159] ;

[0160] ;

[0161] In the formula, , , These are approximate solutions for the pressure field, temperature field, and displacement field, respectively. The number of solid units, cave units, or fissure units; The number of solid elements or crack elements; The total number of nodes in a body element; , , These are the shape functions for pressure, temperature, and displacement, respectively. and These are the pressure and temperature values ​​at the center of the unit or the center of the slotted unit, respectively. This represents the nodal displacement value.

[0162] S203. The momentum conservation equation is spatially discretized using the finite element method, the mass conservation equation and the energy conservation equation are spatially discretized using the finite volume method, and the time discretization is performed using the finite difference method.

[0163] S204. Substitute the shape functions into each governing equation and linearize them using the Newton-Raphson method to obtain the residual form of the thermo-fluid-structure interaction equations, thus forming a set of residual equations.

[0164] In some embodiments, the residual equation set includes displacement residual equation, matrix pressure residual equation, matrix temperature residual equation, fracture pressure residual equation, fracture temperature residual equation, and cavern pressure residual equation.

[0165] The functional expression of the residual equation system is as follows:

[0166] ;

[0167] ;

[0168] ;

[0169] ;

[0170] ;

[0171] ;

[0172] In the formula, , , , , , These represent the displacement residual vector, matrix pressure residual, matrix temperature residual, fracture pressure residual, fracture temperature residual, and cavern pressure residual of element a, respectively; n is the previous time step. This is the current time step; Let be the nodal displacement vector of element a; These are the nodal displacement vectors of element a in the current time step and the previous time step, respectively; , These are the strain-displacement matrices for solid elements a and b, respectively. It is the elasticity matrix; Let be the displacement field shape function of element a; Let be the gradient of the displacement shape function of element a within the volume domain; Let i be the pressure field shape function of the i-th matrix element; Let i be the temperature field shape function of the i-th matrix element; Let be the pressure field shape function of the i-th crack element; Let i be the pressure field shape function of the i-th karst cave unit; Let i be the pressure field shape function of the i-th karst cave unit; Let be the rock mass density value of the i-th unit in the k-th iteration step of the current time step; The time step is the time interval between two adjacent time steps. Let be the thermal conductivity transfer coefficient between the i-th unit and the j-th unit in the matrix at the contact surface; The flow potential at the contact surface between the crack element and the matrix element at the current time step; This represents the flow potential at the contact surface between two adjacent matrix elements at the current time step. This represents the total number of fracture-cavity contact surface units. Let be the volume of the i-th cavern; Let be the equivalent bulk modulus of the i-th karst cave;

[0173] S205. Based on the residual equation system, a nonlinear algebraic equation system is constructed. After linearization using the Newton-Raphson method and assembly of the Jacobian matrix, the rock deformation field, pressure field and temperature field are obtained.

[0174] The rock deformation field includes rock matrix displacement;

[0175] The pressure field includes rock matrix pressure, fluid pressure within fractures, and fluid pressure within caverns.

[0176] The temperature field includes the temperature distribution in the matrix and the temperature distribution in the cracks;

[0177] In some embodiments, the iterative solution format for the nonlinear algebraic equation system is as follows:

[0178] ;

[0179] In the formula, k is the identifier of the Newton-Raphson iteration step; k+1 is the current iteration step, and k is the previous iteration step; The increment vector of the main variables for the current iteration step is obtained by solving the following system of linear equations:

[0180] ;

[0181] ;

[0182] In the formula, It is a Jacobian matrix; This is the increment vector of the main variables at the current time step and the k-th iteration step; This is the residual vector at the current time step and the kth iteration step.

[0183] This embodiment employs an explicit scheme to estimate the specific terms in each block of the Jacobian matrix. The coefficients of the main variables are all considered known quantities and are taken from the parameters that have converged in the previous time step, which significantly reduces the computational cost of the Jacobian matrix and improves the solution efficiency. In order to numerically solve the coupled system, the control equations are transformed into a series of linearized algebraic equations solved by the Newton-Raphson method after the aforementioned spatiotemporal discretization. Therefore, when the initial estimate is close enough to the true solution of the equations, the solution of the coupled equations usually only requires 3 to 5 Newton-Raphson iterations to meet the convergence condition.

[0184] S206. Obtain the preset fluid pressure inside the wellbore, and use the fluid pressure inside the fracture, the fluid pressure inside the karst cave, and the fluid pressure inside the wellbore as the stress boundary conditions for the induced stress calculation model.

[0185] In this embodiment, the fluid pressure inside the wellbore is determined independently based on the injection parameters of the fracturing operation and the wellbore flow model; alternatively, a preset fluid pressure inside the wellbore can be used directly as the boundary condition.

[0186] S3. Solve the stress-induced stress calculation model and stress boundary conditions to obtain the combined induced stress field of the wellbore-fracture-cavity and the fracture physical parameters.

[0187] Among them, the physical properties of the crack include crack aperture, crack permeability and crack stiffness;

[0188] This embodiment solves the industry pain point that conventional methods cannot accurately characterize the stress transition of cavity structures by using the displacement discontinuity mechanical characteristics of cracks and the virtual stress method for the stress discontinuity mechanical characteristics of wells and karst caves.

[0189] like Figures 2-3 As shown, specifically, step S3 includes the following steps:

[0190] Based on the induced stress calculation model, the displacement discontinuity method is used to perform surface integration on the crack element through the Kelvin fundamental solution to construct the influence coefficient relationship between the displacement discontinuity of the crack element and the induced stress generated in space.

[0191] Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the wellbore unit and the induced stress generated in space.

[0192] Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the karst unit and the induced stress generated in space.

[0193] For all crack elements, determine whether they are in a closed state. For crack elements in a closed state, force their normal displacement discontinuity to be 0, and correct their shear stress constraint based on the Mohr-Coulomb criterion to avoid non-physical calculation results.

[0194] The functional expression for the Mohr-Coulomb criterion is:

[0195] ;

[0196] ;

[0197] In the formula, and They act on the first Total shear stress and effective normal stress on each crack element; The cohesive strength of the crack surface. The coefficient of friction of the crack surface; and The first The shear stress along the direction and the shear stress diagonally of each crack element;

[0198] The shear stress constraint of the modified crack element is expressed as follows:

[0199] ;

[0200] ;

[0201] ;

[0202] In the formula, The direction of the total shear stress on the crack surface;

[0203] The displacement discontinuities of all fracture elements, the virtual stresses of all wellbore elements, and the virtual stresses of all karst elements are taken as the basic unknowns.

[0204] The fluid pressure in the fracture, the fluid pressure in the cave, and the fluid pressure in the wellbore are applied to the corresponding fracture unit, cave unit, and wellbore unit to establish a joint linear equation system.

[0205] Solving the combined linear equations yields the displacement discontinuities of each fracture element, the virtual stress of each wellbore element, and the virtual stress of each karst cave element.

[0206] Based on the displacement discontinuities and virtual stresses obtained from the solution, the induced stresses generated in space by the cracks, well shafts and karst caves are calculated respectively through the influence coefficient relationship;

[0207] The induced stress field of the karst cave is expressed as follows:

[0208] ;

[0209] In the formula, , and Each cave unit is a separate unit for all cave units. Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; For cave units Virtual stress at the direction of the cave unit The influence coefficient of shear-induced stress generated at the location; For cave units Virtual stress at the direction of the cave unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units Virtual stress at the direction of the cave unit The influence coefficient of the normal induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of shear-induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of the normal induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of shear-induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of the normal induced stress generated at the location; and Separately, it is a cave unit. Virtual stress along the dip, strike, and normal directions; This represents the total number of all cave units.

[0210] The induced stress field in the wellbore is expressed as follows:

[0211] ;

[0212] In the formula, , and For all wellbore units, in the wellbore unit Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; Wellbore Unit Virtual stress at the location in the wellbore unit The influence coefficient of shear-induced stress generated at the location; Wellbore Unit Virtual stress at the location in the wellbore unit The influence coefficient of the tendency shear-induced stress generated at the location; Wellbore Unit Virtual stress at the location in the wellbore unit The influence coefficient of the normal induced stress generated at the location; Wellbore Unit Virtual stress at the wellbore element The influence coefficient of shear-induced stress generated at the location; Wellbore Unit Virtual stress at the wellbore element The influence coefficient of the tendency shear-induced stress generated at the location; Wellbore Unit Virtual stress at the wellbore element The influence coefficient of the normal induced stress generated at the location; Wellbore Unit Normal virtual stress at the wellbore element The influence coefficient of shear-induced stress generated at the location; Wellbore Unit The normal virtual stress at the wellbore element The influence coefficient of the tendency shear-induced stress generated at the location; Wellbore Unit Normal virtual stress at the wellbore element The influence coefficient of the normal induced stress generated at the location; and Not a wellbore unit Virtual stress along the dip, strike, and normal directions; This represents the total number of all wellbore units.

[0213] The induced stress field of the crack is expressed as follows:

[0214] ;

[0215] In the formula, , and For all crack elements, in the crack element Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of the normal induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of the normal induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of the normal induced stress generated at the location; and Other units The discontinuous displacement along the dip, strike, and normal directions; The total number of all crack elements. It is a collection of closed cracks.

[0216] By linearly superimposing all induced stresses, a combined induced stress field of wellbore-fracture-cavity is obtained;

[0217] The functional expression for the combined induced stress field of the wellbore-fracture-cavity system is as follows:

[0218] ;

[0219] In the formula, Take the cell type , or , The fluid pressure inside the cave. The fluid pressure inside the wellbore. The fluid pressure within the crack; The direction component of the total induced shear stress acting on element i; The tendency component of the total induced shear stress acting on element i; This represents the total induced normal stress acting on element i;

[0220] Based on the induced stress field of the wellbore-fracture-cavity system, the displacement discontinuity in the normal direction of the fracture unit is obtained to determine the fracture aperture.

[0221] The formula for calculating the crack aperture is:

[0222] ;

[0223] In the formula, The crack opening;

[0224] Based on the law of cubes, the fracture permeability is calculated according to the fracture aperture.

[0225] The formula for calculating crack permeability is as follows:

[0226] ;

[0227] In the formula, The permeability of the crack; This is a correction factor related to factors such as crack surface roughness;

[0228] The crack stiffness is calculated based on the crack opening and the fluid pressure inside the crack.

[0229] The formula for calculating crack stiffness is:

[0230] ;

[0231] In the formula, For crack stiffness; This represents the change in internal fluid pressure within the fracture element; This represents the change in the aperture of the crack element.

[0232] S4. During the iterative solution process at the current time step, compare the main variable field and crack opening in the current iteration step with those in the previous iteration step, and determine whether the preset convergence condition is met.

[0233] If the conditions are not met, the parameters in the thermal-fluid-structure interaction mathematical model are updated based on the crack physical property parameters, and the process returns to step S2; if the conditions are met, step S5 is executed.

[0234] In this embodiment, the function expression for the preset convergence condition is:

[0235] ;

[0236] ;

[0237] ;

[0238] ;

[0239] In the formula, n is the previous time step; n+1 is the current time step; and k is the current iteration step. Let L2 be the norm of the vector; , , and These are the convergence thresholds for pressure difference, temperature difference, displacement difference, and crack opening difference, respectively.

[0240] In some embodiments, the parameter update of the thermal-fluid-structure interaction mathematical model specifically involves:

[0241] The calculated crack aperture, crack permeability, and crack stiffness are updated in the corresponding crack element properties in the thermo-fluid-structure interaction mathematical model as input parameters for the next iteration.

[0242] S5. Determine whether the crack meets the propagation criterion based on the joint induced stress field to determine whether the current crack geometry needs to be updated.

[0243] If satisfied, determine the propagation direction and amount of the crack front edge, and update the current crack geometry.

[0244] If the conditions are not met, the current crack geometry will remain unchanged.

[0245] In some embodiments, the propagation criterion may adopt the maximum circumferential stress criterion commonly used in rock mechanics; the direction of crack propagation is determined by the maximum circumferential stress criterion: the crack propagates along the direction of the maximum circumferential normal stress, that is, the direction in which the circumferential shear stress is zero; the amount of crack propagation is calculated by combining the ratio of stress intensity factor to rock fracture toughness and the volume of injected fluid.

[0246] In some embodiments, updating the current crack geometry specifically involves adding new crack elements at the crack leading edge that satisfies the propagation criterion, according to the calculated propagation direction and propagation amount, to complete the update of the current crack geometry and prepare for the solution of the next time step.

[0247] S6. Using the current crack geometry as the input for the next time step, return to steps S2 to S5 until the preset simulation duration is reached, and obtain the crack propagation trajectory and the corresponding master variable field.

[0248] In some embodiments, when the simulation duration reaches the preset total duration, the output results include time history data of the three-dimensional propagation trajectory of the fracture throughout the entire simulation cycle, the rock mass deformation field, pressure field, temperature field, the joint induced stress field, and fracture physical parameters at each moment, providing a complete numerical simulation basis for optimizing the fracturing construction scheme of fractured-vuggy oil and gas reservoirs.

[0249] To verify the computational accuracy and engineering applicability of the coupled fracturing simulation method for fractured-cavity media wellbore-fracture-cavity described in Example 1, this implementation constructs two sets of benchmark verification models, with the specific settings as follows:

[0250] Benchmark Validation Model 1: Calculation Example of Stress Field Validation for Cylindrical Wellbore;

[0251] like Figure 4 As shown, a numerical model of a cylindrical wellbore in infinite three-dimensional space is constructed. The wellbore length is 20m and the radius is 1.0m. The initial pore pressure of the rock matrix is ​​1.5MPa, and the far-field in-situ stress is... 6.5MPa and A constant internal pressure of 2.5 MPa was applied to the inner wall of the wellbore at 1.5 MPa. The rock mass mechanical parameters were set as follows: Young's modulus 38.8 GPa, Poisson's ratio 0.15. The wellbore wall was discretized into 4708 triangular surface elements, and the rock mass matrix was discretized into hexahedral structured elements. Numerical solutions were obtained using the method described in Example 1, with the Kirsch analytical solution in cylindrical coordinates serving as the verification benchmark.

[0252] like Figure 5 As shown, the verification results indicate that: Figure 5 As shown in (a), the wellbore-induced stress is only significant in a limited area near the borehole (i.e., within 3 times the wellbore radius); beyond this range, the stress disturbance caused by the wellbore is negligible. Figure 5 As shown in (b), the effective stress at the wellbore wall varies sinusoidally with the circumferential angle. Under the condition of fine mesh, the relative error between the numerical solution and the Kirsch analytical solution in this embodiment is less than 2%, which verifies the accuracy of the method for calculating the induced stress in the wellbore.

[0253] Benchmark Validation Model 2: Validation Example of Stress Field in a Spherical Cavern;

[0254] like Figure 6As shown in (a), a numerical model of a spherical karst cave in an infinite elastic space was constructed. The karst cave radius was 15 m, and the internal fluid pressure was 1.0 MPa. The rock mechanics parameters were consistent with those of the benchmark verification model 1. The surface of the karst cave was discretized using 118, 226, and 626 triangular surface elements, respectively. It was assumed that the surrounding medium was an elastic isotropic solid with a Young's modulus of 38.8 GPa and a Poisson's ratio of 0.15. Subsequently, numerical solutions were obtained using the method described in Example 1 at different mesh resolutions, with the analytical solution of the spherical cavity pressure problem in elasticity serving as the verification benchmark.

[0255] like Figure 6 (a) Figure 6 (b) and Figure 6 As shown in (c), the verification results show that the radial stress induced by the cave is compressive stress and the tangential stress is tensile stress. The absolute value of the tangential stress is 1 / 2 of the radial stress. Both of them decrease rapidly with the distance from the center of the cave at a rate of O(r−3). Under different grid resolutions, the numerical solution and the analytical solution of this embodiment are in high agreement, which verifies the accuracy of the method for calculating the cave-induced stress.

[0256] Example 2

[0257] like Figure 7 As shown, this embodiment provides a fracture-cavity type wellbore-fracture-cavity coupled fracturing simulation system 100, including:

[0258] Model building module 101 is used to build a coupled fracturing simulation model of wellbore-fracture-cavity and to discretize the study area into a grid to obtain the volume elements, fracture elements, wellbore elements and cavity elements of the rock matrix;

[0259] The coupled fracturing simulation model includes a thermo-fluid-structure interaction mathematical model and an induced stress calculation model.

[0260] The first solution module 102 is used to solve the problem based on the thermal-fluid-structure interaction mathematical model and the current geometry of the crack to obtain the main variable field and extract stress boundary conditions from the main variable field.

[0261] The main variable fields include rock mass deformation field, pressure field and temperature field;

[0262] Among them, the stress boundary conditions include fluid pressure inside the fracture, fluid pressure inside the cavern, and fluid pressure inside the wellbore;

[0263] The second solution module 103 is used to solve the combined induced stress field and fracture physical parameters of the wellbore-fracture-cavity based on the induced stress calculation model and the stress boundary conditions;

[0264] Among them, the physical properties of the crack include crack aperture, crack permeability and crack stiffness;

[0265] The coupled iteration control module 104 is used to compare the main variable field and crack aperture in the current iteration step and the previous iteration step during the iterative solution process in the current time step, and to determine whether the preset convergence condition is met.

[0266] If the conditions are not met, the parameters in the thermal-fluid-structure interaction mathematical model are updated based on the crack physical property parameters, and the process returns to the first solution module; if the conditions are met, the crack propagation update module is used.

[0267] The crack propagation update module 105 is used to determine whether the crack meets the propagation criterion based on the joint induced stress field, and to determine whether the current crack geometry needs to be updated.

[0268] If the propagation criterion is met, the propagation direction and propagation amount of the crack front are determined according to the joint induced stress field, and the current crack geometry is updated.

[0269] If the expansion criterion is not met, the current crack geometry remains unchanged.

[0270] The trajectory acquisition module 106 is used to take the current crack geometry as the input for the next time step, return to the first solution module, until the preset simulation duration is reached, and obtain the crack propagation trajectory and the corresponding main variable field.

[0271] Example 3

[0272] This embodiment provides an electronic device, including at least one control processor and a memory for communicatively connecting to the at least one control processor;

[0273] Memory, as a non-transitory computer-readable storage medium, can be used to store non-transitory software programs and non-transitory computer-executable programs. Furthermore, memory may include high-speed random access memory, and may also include non-transitory memory, such as at least one disk storage device, flash memory device, or other non-transitory solid-state storage device. In some embodiments, memory may optionally include memory remotely located relative to the processor, and these remote memories can be connected to the processor via a network. Examples of such networks include, but are not limited to, the Internet, intranets, local area networks, mobile communication networks, and combinations thereof.

[0274] The non-transient software program and instructions required to implement the fracture-cavity coupled fracturing simulation method for a fractured-cavity wellbore-fracture-cavity system described in the above embodiments are stored in memory. When executed by a processor, the fracture-cavity coupled fracturing simulation method described in the above embodiments is executed, for example, the method described above is executed. Figure 1 Method steps S1 to S6;

[0275] The system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0276] Example 4

[0277] This embodiment provides a computer-readable storage medium storing computer-executable instructions for causing a computer to execute a coupled fracturing simulation method for fractured-cavity media wellbore-fracture-cavity as described in Embodiment 1.

[0278] It should be noted that the computer-readable storage medium in this embodiment may be a computer-readable signal medium or a computer-readable storage medium, or any combination thereof. The computer-readable storage medium may be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof.

[0279] More specific examples of computer-readable storage media may include, but are not limited to: electrical connections having one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination of the foregoing.

[0280] In this embodiment, the computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device. In this embodiment, the computer-readable signal medium can include a data signal propagated in baseband or as part of a carrier wave, carrying a computer-readable program. Such propagated data signals can take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof.

[0281] The computer-readable signal medium may also be any computer-readable storage medium other than a computer-readable storage medium, which can send, propagate or transmit a program for use by or in connection with an instruction execution system, apparatus or device.

[0282] The computer program contained on the computer-readable storage medium can be transmitted using any suitable medium, including but not limited to: wires, optical fibers, RF (radio frequency), etc., or any suitable combination thereof. The aforementioned computer-readable storage medium can be used to write a computer program for performing this embodiment in one or more programming languages ​​or combinations thereof. These programming languages ​​include object-oriented programming languages—such as Python, Java, and C++—and conventional procedural programming languages—such as C or similar programming languages. The program can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computer (e.g., via the Internet using an Internet service provider).

[0283] The above description is merely a preferred embodiment of the present invention. It will be apparent to those skilled in the art that the present invention is not limited to the details of the above exemplary embodiments, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be considered exemplary and non-limiting in all respects. The scope of the present invention is defined by the appended claims rather than the foregoing description, and thus all variations falling within the meaning and scope of the equivalents of the claims are intended to be included within the present invention.

[0284] To further illustrate the technical solution of this invention, at least one application example is provided, as follows:

[0285] Application Example 1

[0286] This case study focuses on the mechanical interactions between propagating fractures, well shafts, and caverns, examining, for example... Figure 8 (a) and Figure 8 (b) shows a fractured-vuggy reservoir with a range of 50m. 50m The rock matrix is ​​discretized into 20×20×20=8000 Cartesian structured hexahedral elements. A horizontal wellbore extends 20m along the x-direction, with its geometric center coinciding with the initial fracture geometric center. The coordinates of the initial fracture geometric center are... The geometric center coordinates of a spherical cave with a radius of 1.0m and a diameter of (25, 25, 25)m are given. The length is (20,20,20)m; two initial cutting fractures with a spacing of 6m are prefabricated in the horizontal well section.

[0287] This example sets up three sets of cavern internal pressure conditions (25MPa, 27.5MPa, 30MPa) and three sets of cavern size conditions (1.0m, 1.75m, 2.5m), with a total simulation time of 3600s, to investigate the effects of cavern pressure and cavern size on fracture evolution in the porosity elastic reservoir region.

[0288] The main parameters were set as follows: matrix porosity 0.25; initial matrix permeability 0.5 mD; initial matrix pore water pressure 25 MPa; initial fluid pressure inside the wellbore 27 MPa; wellbore radius 0.25 m; Young's modulus of the rock mass 38.8 GPa; Poisson's ratio 0.15; rock fracture toughness 3.0 MPa. In-situ stress is MPa; fracturing fluid injection rate The fracturing fluid viscosity is 5.0 cP.

[0289] like Figure 9 As shown, Figure 9 This study demonstrates the relative positions of the wellbore, caverns, and propagating fractures, as well as the geometric evolution of three-dimensional hydraulic fractures in the absence of natural fractures. All caverns exhibit an attractive effect on nearby fractures, with spherical caverns showing a stronger attraction than ellipsoidal ones, leading to more pronounced fracture deflection. This phenomenon can be attributed to the stress concentration characteristics around caverns of different geometries. Furthermore, the stress shadow effect between adjacent fractures also contributes to their mutually deflecting growth.

[0290] in, Figure 9 (a) shows a schematic diagram of the propagation results of two competing cracks after 5 steps in the absence of natural cracks; Figure 9 (b) shows a schematic diagram of the expansion results of two competing cracks after 10 steps in the absence of natural cracks; Figure 9 (c) shows a schematic diagram of the expansion results of two competing cracks after 15 steps in the absence of natural cracks; Figure 9 (d) shows a schematic diagram of the expansion results of two competing cracks after 20 steps in the absence of natural cracks;

[0291] At the same time, in comparison Figure 9 and Figure 10 The crack propagation morphology reveals that existing natural cracks not only affect the crack propagation process but also the final crack geometry. At this point, the attraction of the karst cave to the propagating crack is weakened, which is partly attributed to the slippage of the natural crack, which induces additional stress disturbances and may thus hinder the crack from propagating towards the karst cave.

[0292] like Figure 10 As shown, Figure 10The diagram shows the temperature and pressure field contours in the presence of existing natural cracks. Figure 10 (a) is a temperature field contour map; Figure 10 (b) shows the pressure field cloud map. It can be observed that about an hour after hydraulic fracturing, only a very small area experienced cooling. Therefore, as long as a large number of secondary fractures induced by cooling are not generated near the main fracture, the effect of thermally induced stress on fracture propagation during the fracturing stage can be ignored. Most natural fractures are connected to propagation fractures, and the fluid pressure in these connected fractures can reach up to 32 MPa. In contrast, only a few isolated natural fractures and the tip region of propagation fractures still maintain a low fluid pressure.

Claims

1. A coupled fracturing simulation method for fractured-cavity media wellbore-fracture-cavity, characterized in that, Includes the following steps: A. Construct a coupled fracturing simulation model of wellbore-fracture-cavity, and discretize the study area into a grid to obtain the volume elements, fracture elements, wellbore elements and cavity elements of the rock matrix; The coupled fracturing simulation model includes a thermo-fluid-structure interaction mathematical model and an induced stress calculation model; B. Solve based on the aforementioned thermal-fluid-structure interaction mathematical model and the current geometry of the crack to obtain the principal variable field, and extract stress boundary conditions from the principal variable field; The main variable fields include the rock mass deformation field, pressure field, and temperature field; The stress boundary conditions include fluid pressure inside the fracture, fluid pressure inside the cavern, and fluid pressure inside the wellbore. C. Solve the stress-induced stress calculation model and stress boundary conditions to obtain the combined induced stress field of the wellbore-fracture-cavity and the fracture physical parameters. The fracture physical properties include fracture aperture, fracture permeability, and fracture stiffness; D. During the iterative solution process at the current time step, compare the main variable field and crack opening in the current iteration step with those in the previous iteration step, and determine whether the preset convergence condition is met. If the conditions are not met, then based on the crack physical properties, update the parameters in the thermal-fluid-structure interaction mathematical model and return to step B; if the conditions are met, then proceed to step E. E. Use the current crack geometry as the input for the next time step, return to step B, and continue until the preset simulation duration is reached, and obtain the crack propagation trajectory and the corresponding main variable field. Step C includes: Based on the induced stress calculation model, the displacement discontinuity method is used to perform surface integration on the crack element through the Kelvin fundamental solution to construct the influence coefficient relationship between the displacement discontinuity of the crack element and the induced stress generated in space. Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the wellbore unit and the induced stress generated in space. Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the karst unit and the induced stress generated in space. The displacement discontinuities of all fracture elements, the virtual stresses of all wellbore elements, and the virtual stresses of all karst elements are taken as the basic unknowns. The fluid pressure in the fracture, the fluid pressure in the cave, and the fluid pressure in the wellbore are applied to the corresponding fracture unit, cave unit, and wellbore unit to establish a joint linear equation system. Solving the combined linear equations yields the displacement discontinuities of each fracture element, the virtual stress of each wellbore element, and the virtual stress of each karst cave element. Based on the displacement discontinuities and virtual stresses obtained from the solution, the induced stresses generated in space by the cracks, well shafts and karst caves are calculated respectively through the influence coefficient relationship; By linearly superimposing all induced stresses, a combined induced stress field of wellbore-fracture-cavity is obtained.

2. The method according to claim 1, characterized in that, The steps following step D and before step E include: Based on the joint induced stress field, it is determined whether the crack meets the propagation criterion, and whether the current crack geometry needs to be updated. If satisfied, determine the propagation direction and amount of the crack front edge, and update the current crack geometry. If the conditions are not met, the current crack geometry will remain unchanged.

3. The method according to claim 1, characterized in that, The thermal-fluid-structure interaction mathematical model includes the momentum conservation equation, the heat and mass transfer integral equation of the rock matrix, the heat and mass transfer integral equation of the fracture, and the mass conservation equation of the karst cave. The integral form of the momentum conservation equation is: ； ； In the formula, A is the three-dimensional integral region; A is the two-dimensional surface integral region. Let i be the component of the concentrated force along the i-th direction; Let i be the component of the surface force along the i-th direction; , and These are the areas affected by surface forces, the surface of the cave, and the surface of the cracks, respectively. This is the contact stress tensor at the crack wall; , and These are the fluid pressures within the rock matrix, fissures, and caverns, respectively. The external normal vector of the cave surface; is the outward normal vector of the crack surface; , is the component of the outward normal vector of the surface of the cave or crack along the i-direction; For the volumetric forces of the rock mass matrix; The volume density of the rock matrix; The component of gravitational acceleration along the i-direction; , and These represent the displacements along the j, k, and l directions, respectively. Second-order unit tensor; For the displacement elastic tensor; The coefficient of thermal expansion; It is the bulk modulus of the drainage. This refers to the change in temperature. Let k be the component of the rock mass displacement vector in the k direction; Let be the component of the rock mass displacement vector in the l direction.

4. The method according to claim 3, characterized in that, The integral equation for heat and mass transfer in the rock matrix is ​​expressed as follows: ； ； ； ； In the formula, For temperature; Entropy; The coefficient of thermal expansion of the rock matrix; This represents the temperature change within the rock matrix. This represents the temperature change within the crack. The equivalent volumetric heat capacity of the rock mass; For the volumetric strain of the rock mass; Biot modulus; H is the bulk modulus of the rock mass; H is the specific enthalpy of the fluid. This is the rock matrix permeability tensor; Let j be the component of gravitational acceleration along the j-direction; For fluid viscosity; is the component of the outward normal vector of the surface of the cave or crack along the j direction; These are the source and sink terms of mass in the rock mass matrix; It represents the source and sink terms of enthalpy within the rock matrix. Let the thermal conductivity tensor be... For simulating time; The integral equation for heat and mass transfer in the crack is expressed as follows: ； ； ； ； In the formula, and These are crack stiffness and fluid stiffness, respectively. The crack opening; The coefficient of thermal expansion of the fluid; Specific heat capacity of the fluid; The temperature inside the crack; Let f be the fracture permeability tensor; It is the enthalpy of the fluid; and These are the source and sink terms for the mass and enthalpy within the crack, respectively. This is the region of crack integration; Divide the contact area between the cave and the fissure into zones; The crack-matrix contact area is divided into regions; The mass conservation equation for the karst cave is expressed as follows: ； In the formula, This is the integral region of the karst cave; Let i be the modulus of the i-th karst cave. It is a source and sink of mass within the karst cave.

5. The method according to claim 1, characterized in that, Step B includes: Based on the aforementioned volume element, cavern element, and crack element, the momentum conservation equation in the thermal-fluid-structure interaction mathematical model is spatially discretized using the finite element method, the mass conservation equation and energy conservation equation are spatially discretized using the finite volume method, and time discretized using the finite difference method to obtain a set of nonlinear algebraic equations. The nonlinear algebraic equations were linearized using the Newton-Raphson method, and the Jacobian matrix was assembled to solve the linearized algebraic equations, thereby obtaining the rock mass deformation field, pressure field, and temperature field. The rock mass deformation field includes the displacement of the rock mass matrix; The pressure field includes rock matrix pressure, fluid pressure within fractures, and fluid pressure within karst caves. The temperature field includes the temperature distribution within the rock matrix and the temperature distribution within the fractures; Obtain the preset fluid pressure inside the wellbore, and use the fluid pressure inside the fracture, the fluid pressure inside the cavern, and the fluid pressure inside the wellbore as the stress boundary conditions for the induced stress calculation model.

6. The method according to claim 1, characterized in that, Step C further includes: Based on the induced stress field of the wellbore-fracture-cavity system, the displacement discontinuity in the normal direction of the fracture unit is obtained to determine the fracture aperture. Based on the law of cubes, the fracture permeability is calculated according to the fracture aperture. The crack stiffness is calculated based on the crack opening and the fluid pressure inside the crack.

7. The method according to claim 1, characterized in that, The induced stress field of the karst cave is expressed by the following function: ； In the formula, , and Each cave unit is a separate unit for all cave units. Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; For cave units Virtual stress at the direction of the cave unit The influence coefficient of shear-induced stress generated at the location; For cave units Virtual stress at the direction of the cave unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units Virtual stress at the direction of the cave unit The influence coefficient of the normal induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of shear-induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units Virtual stress tends at the location in the cavern unit The influence coefficient of the normal induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of shear-induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of the tendency shear-induced stress generated at the location; For cave units The normal virtual stress at the location in the cavern unit The influence coefficient of the normal induced stress generated at the location; and Separately, it is a cave unit. Virtual stress along the dip, strike, and normal directions; The total number of all cave units; The induced stress field of the crack is expressed as follows: ； In the formula, , and For all crack elements, in the crack element Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous displacement at the location in the crack element The influence coefficient of the normal induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous displacement at the crack element The influence coefficient of the normal induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of shear-induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; For crack elements Discontinuous normal displacement at the crack element The influence coefficient of the normal induced stress generated at the location; and Do not use crack units The discontinuous displacement along the dip, strike, and normal directions; The total number of all crack elements. It is a collection of closed cracks.

8. A coupled fracturing simulation system for fractured-cavity media wellbore-fracture-cavity, characterized in that, include: The model building module is used to construct a coupled fracturing simulation model of wellbore-fracture-cavity, and to discretize the study area into a grid to obtain the volume elements, fracture elements, wellbore elements and cavity elements of the rock matrix; The coupled fracturing simulation model includes a thermo-fluid-structure interaction mathematical model and an induced stress calculation model; The first solution module is used to solve the problem based on the thermal-fluid-structure interaction mathematical model and the current geometry of the crack to obtain the main variable field and extract the stress boundary conditions from the main variable field. The main variable fields include the rock mass deformation field, pressure field, and temperature field; The stress boundary conditions include fluid pressure inside the fracture, fluid pressure inside the cavern, and fluid pressure inside the wellbore. The second solution module is used to solve the stress field and fracture physical parameters of the wellbore-fracture-cavity joint induced stress field based on the induced stress calculation model and the stress boundary conditions. The fracture physical properties include fracture aperture, fracture permeability, and fracture stiffness; The coupled iterative control module is used to compare the main variable field and crack aperture in the current iteration step and the previous iteration step during the iterative solution process in the current time step, and to determine whether the preset convergence condition is met. If the conditions are not met, the parameters in the thermal-fluid-structure interaction mathematical model are updated based on the crack physical property parameters, and the process returns to the first solution module; if the conditions are met, the trajectory acquisition module is triggered. The trajectory acquisition module is used to take the current crack geometry as the input for the next time step, return to the first solution module, until the preset simulation duration is reached, and obtain the crack propagation trajectory and the corresponding main variable field. The process of solving the thermal-fluid-structure interaction mathematical model and the current crack geometry to obtain the principal variable field, and extracting stress boundary conditions from the principal variable field, specifically includes the following steps: Based on the induced stress calculation model, the displacement discontinuity method is used to perform surface integration on the crack element through the Kelvin fundamental solution to construct the influence coefficient relationship between the displacement discontinuity of the crack element and the induced stress generated in space. Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the wellbore unit and the induced stress generated in space. Based on the induced stress calculation model, the virtual stress method is used to perform surface integration on the wellbore unit and the karst unit respectively through the Kelvin fundamental solution, and to construct the influence coefficient relationship between the virtual stress of the karst unit and the induced stress generated in space. The displacement discontinuities of all fracture elements, the virtual stresses of all wellbore elements, and the virtual stresses of all karst elements are taken as the basic unknowns. The fluid pressure in the fracture, the fluid pressure in the cave, and the fluid pressure in the wellbore are applied to the corresponding fracture unit, cave unit, and wellbore unit to establish a joint linear equation system. Solving the combined linear equations yields the displacement discontinuities of each fracture element, the virtual stress of each wellbore element, and the virtual stress of each karst cave element. Based on the displacement discontinuities and virtual stresses obtained from the solution, the induced stresses generated in space by the cracks, well shafts and karst caves are calculated respectively through the influence coefficient relationship; By linearly superimposing all induced stresses, a combined induced stress field of wellbore-fracture-cavity is obtained.

9. The system according to claim 8, characterized in that, Following the coupling iteration control module and preceding the trajectory acquisition module, the system further includes: The crack propagation update module is used to determine whether the crack meets the propagation criterion based on the joint induced stress field, and to determine whether the current crack geometry needs to be updated. If the propagation criterion is met, the propagation direction and propagation amount of the crack front are determined according to the joint induced stress field, and the current crack geometry is updated. If the expansion criterion is not met, the current crack geometry remains unchanged.

Images