Wellbore near wellbore hydraulic fracture deformation and expansion thermo-hydro-mechanical coupling simulation method and system

Abstract

The application provides a near-wellbore hydraulic fracture deformation and expansion thermo-hydro-mechanical coupling simulation method and system, comprising: generalizing reservoir rock mass as matrix and fracture, establishing a thermo-hydro-mechanical coupling mathematical model, describing the mechanical process of fluid flow, heat conduction and convection, rock mass deformation and fracture expansion; adopting three-dimensional displacement discontinuity method (DDM) and virtual stress method (FSM) to calculate fracture and wellbore induced stress respectively, calculating fluid flow, temperature change, matrix deformation and matrix heat exchange of the matrix and the fracture through finite volume and finite element method; based on linear elastic fracture mechanics, combining the maximum circumferential stress criterion and Paris scaling law to calculate the fracture expansion direction and speed; adopting dynamic grid fracture front tracking technology to update the fracture morphology, forming a three-dimensional non-planar fracture expansion model, and realizing the whole process simulation of the near-wellbore hydraulic fracture.

Description

Technical Field

[0001] This application belongs to the field of oil and gas reservoir development, specifically involving a thermal-fluid-structure interaction simulation method and system for near-wellbore hydraulic fracture deformation and propagation. Background Technology

[0002] Due to the high in-situ stress and low porosity and permeability of shale oil and gas or tight oil and gas reservoirs, generating hydraulic fractures through horizontal well staged fracturing technology is crucial for increasing reservoir permeability, improving the contact area between the reservoir and fractures, and enhancing the connectivity between hydraulic and artificial fractures. During fracturing, under the high pressure of the fracturing fluid, hydraulic fractures initiate and propagate near the wellbore wall, eventually penetrating into the reservoir. However, in the near-wellbore region, the presence of the wellbore itself induces a stress field, creating localized tensile or compressive stress zones around the wellbore that vary with angle (Feng et al., 2018). This induced stress, superimposed on the far-field stress, significantly affects the magnitude and direction of the principal stresses in the local area, thus altering the fracture propagation rate and direction. Furthermore, the methods of fracturing and perforation also create stress concentration (Behrmann and Elbel, 1991), further exacerbating the complexity of near-wellbore stress and significantly impacting the fracture initiation and propagation process. Therefore, the stress concentration generated by the wellbore and fractures to some extent complicates the distribution of geostress, leading to significant tortuosity of artificial fractures and localized tearing cracks (Li et al., 2024), resulting in more complex geometry of hydraulic fractures. The directional propagation of near-wellbore fractures significantly increases the frictional resistance of fracturing fluid, raising the propagation pressure of hydraulic fractures, which in turn brings a series of problems such as high wellhead pressure, sand plugging, and difficulty in deep fracture penetration (Wei et al., 2024). Accurately predicting the stress distribution in the near-wellbore zone and elucidating the near-wellbore hydraulic fracturing initiation and propagation mechanisms have become critical issues that urgently need to be addressed in the development of deep shale oil and gas and tight oil and gas.

[0003] In the study of hydraulic fracturing fracture propagation mechanisms and simulation methods, scholars both domestically and internationally have conducted fruitful research and obtained some insights and understandings that are of great guiding significance for fracturing production. Zhang Guangqing and Chen Mian (2006), based on the principle of stress superposition and the maximum principal stress criterion, proposed a near-wellbore fracture turning model that does not consider the influence of fracture height. This pseudo-three-dimensional model does not yet consider factors such as fracture tip stress concentration, inter-fracture stress interference, and porosity elasticity effects caused by fracturing fluid loss, and requires further development and improvement. Wei Shiming et al. (2024) used the phase-field method to study the near-wellbore fracture turning mechanism in a two-dimensional case. This study found that high ground stress difference is conducive to the turning and propagation of hydraulic fractures, while increasing the perforation angle increases the risk of sand blockage. Using the displacement discontinuity method, Abdollahipour et al. (2015) and Zhang et al. (2011) developed two-dimensional fracture propagation models and studied near-wellbore fracture turning propagation and multi-fracture competitive propagation. Their research showed that the larger the angle between the initial fracture and the maximum horizontal principal stress, the higher the corresponding fracture propagation pressure. The stress shading effect between fractures inhibits the growth of some fractures, exhibiting significant differential propagation. Similar two-dimensional fracture propagation models include Chen and Economides (1999), Cherny et al. (2009), and Cheng et al. (2022). Chen et al. (2020) constructed a three-dimensional planar multi-fracture propagation model using the three-dimensional DDM method. This model uses an explicit integral algorithm to solve the fluid-structure interaction equations and calculates the fracture propagation rate using the asymptotic analytical solution at the fracture tip. It should be noted that this model does not consider the stress discontinuity and stress concentration effects generated in the wellbore; it is only applicable to propagation prediction when the fracture is far from the wellbore. Tang et al. (2019) developed a three-dimensional fracture propagation model based on the three-dimensional DDM method. This model can simulate the turning and propagation of three-dimensional non-planar fractures, but it cannot consider the effects of wellbore stress concentration and pore elasticity on fracture propagation and width distribution. To consider the stress concentration effect of the wellbore in the three-dimensional case, Nikolskiy et al. (2016) used the boundary element method to analyze the stress distribution characteristics around the longitudinal fracture of the radial double wing near the wellbore. However, this work does not allow fracture propagation and needs further improvement. Cherny et al. (2016) further developed a three-dimensional fracture propagation model based on the two-dimensional model. This model can be used to simulate the non-planar propagation of transverse and longitudinal fractures near the wellbore, but it has not yet considered the influence of pore elasticity. Rahman et al. (2002) combined analytical methods, true triaxial fracturing experiments, and boundary element numerical simulation to deeply analyze the control effect of stress interference of the main fracture on the turning and propagation process of transverse fractures near the wellbore and the propagation of tearing secondary fractures. The study showed that closely spaced segmented fracturing will significantly increase the propagation pressure of the fracture.

[0004] The above literature review, both domestic and international, reveals a lack of three-dimensional nonplanar crack propagation calculation models that consider wellbore-fracture stress concentration and stress shadowing effects. There is an urgent need to innovate corresponding numerical simulation methods and develop multi-field coupled near-wellbore full three-dimensional crack propagation models. Currently, the main challenges lie in the following two aspects:

[0005] (1) The calculation method for stress concentration in the wellbore and stress interference between propagating fractures in three-dimensional cases is not mature;

[0006] Currently, the most commonly used methods for calculating crack propagation include the finite element method (FEM) and the boundary element method (BEM). Both methods have certain limitations in simulating wellbore-crack induced stress. Numerical models based on the extended finite element method typically treat the wellbore as the stress-loaded inner boundary and use enrichment functions to describe crack deformation and stress concentration at the crack tip. This method exhibits a significant mesh-scale effect when calculating wellbore induced stress, and its computational accuracy often falls short of capturing the wellbore stress concentration effect. Furthermore, when dealing with three-dimensional multi-crack propagation problems, this method often requires highly sophisticated mesh generation techniques, and the difficulty of code implementation and program stability pose significant challenges to its three-dimensional application. Crack propagation models based on the BEM often neglect the wellbore stress concentration effect, and models that do consider the wellbore stress concentration effect are limited to solving two-dimensional problems. Overall, the main reason for this is the immaturity of methods for calculating wellbore stress concentration and stress interference between propagating cracks in three-dimensional scenarios, necessitating the development of new numerical simulation techniques.

[0007] (2) A three-dimensional thermal-fluid-structure interaction calculation model for near-wellbore fracture propagation has not yet been developed;

[0008] The boundary element method (BEM) is widely accepted in the petroleum industry for solving hydraulic fracture propagation problems. However, existing literature often neglects the thermal-fluid-structure interaction (TFI) effect of the matrix. Due to the mass and heat exchange between the matrix and hydraulic fractures, the fracture propagation rate, fracture aperture, and fracture geometry are significantly affected. Especially for repeated fracturing, the propagation of matrix pore pressure and changes in the temperature field within the matrix influence fracture propagation and direction, and the TFI mechanism requires further investigation. Therefore, employing appropriate numerical processing techniques to address the heat and mass transfer problem between hydraulic fractures and the matrix within the BEM framework, and subsequently establishing a near-wellbore three-dimensional full-dimensional fracture propagation TFI calculation model, has become one of the key challenges urgently needing to be solved in the field of numerical simulation for fracturing stimulation of shale oil and gas and tight oil and gas reservoirs. Summary of the Invention

[0009] The purpose of this application is to provide a thermo-fluid-structure interaction simulation method and system for near-wellbore hydraulic fracture deformation and propagation, which can accurately simulate the near-wellbore hydraulic fracture deformation and propagation process.

[0010] The technical solution provided in this application is as follows:

[0011] In a first aspect, this application provides a thermal-fluid-structure interaction simulation method for near-wellbore hydraulic fracture deformation and propagation, including:

[0012] S1. The reservoir rock mass is generalized into two parts: matrix and fracture. A near-wellbore fracture propagation thermal-fluid-solid coupled mathematical model is established to describe the mechanical processes of fluid flow, heat conduction and heat convection, rock mass deformation and fracture propagation within the matrix and fracture.

[0013] S2. Based on the aforementioned thermal-fluid-structure interaction mathematical model, the three-dimensional displacement discontinuity method (DDM) and the virtual stress method (FSM) are used to calculate the fracture-induced stress and wellbore-induced stress, respectively; the finite volume method and the finite element method (FV / FEM) are used to calculate the fluid flow, temperature change, matrix deformation, and fluid mass and heat exchange between the matrix and the fracture.

[0014] S3. Within the theoretical framework of linear elastic fracture mechanics, and combined with the calculation results of S2, the maximum circumferential stress criterion and the Paris scaling law model are used to calculate the crack propagation direction and propagation speed.

[0015] S4. Employing dynamic mesh fracture propagation leading edge tracking technology, the fracture geometry is updated in real time based on the fracture propagation direction and speed, forming a three-dimensional non-planar fracture propagation calculation model under wellbore boundary constraints, realizing the full-process deformation and propagation simulation of near-wellbore hydraulic fractures from initiation and turning to deep reservoir penetration.

[0016] In one possible implementation, the thermal-fluid-solid coupling mathematical model includes: fluid flow control equations in the matrix and fractures, fluid mass exchange equations between the matrix and fractures, energy conservation equations for the matrix portion of the rock mass, heat flux equations in the matrix, energy conservation equations and heat flux equations for the fracture portion, energy exchange equations between the matrix and fractures, momentum conservation equations for the matrix portion, constitutive equations for the reservoir rock mass, calculation equations for induced stress and displacement discontinuities, fracture propagation prediction equations, and calculation equations for fracture propagation length and propagation direction.

[0017] In one possible implementation, the three-dimensional displacement discontinuity method (DDM) and the virtual stress method (FSM) are used to calculate the fracture-induced stress and wellbore-induced stress, respectively. This includes: uniformly dividing the fracture and wellbore into a series of triangular meshes; each fracture triangular calculation unit contains one normal displacement discontinuity and two shear direction displacement discontinuities; the induced stress of each fracture triangular calculation unit is obtained through the mapping relationship between the displacement discontinuities and the induced stress; each wellbore triangular calculation unit contains one normal virtual stress and two shear direction virtual stresses; the induced stress of each wellbore triangular calculation unit is obtained through the mapping relationship between the virtual stress and the induced stress.

[0018] In one possible implementation, the rock mass deformation includes matrix deformation and crack deformation, wherein the crack deformation includes shear deformation and crack width variation, which are described by the normal displacement discontinuity and shear direction displacement discontinuity of the crack triangle calculation unit, respectively.

[0019] In one possible implementation, a dynamic mesh crack propagation leading edge tracking technique is employed to update the crack geometry in real time based on the crack propagation direction and velocity, including:

[0020] (i) Based on the crack propagation length and propagation direction in the current time step, a series of initial propagation vectors are obtained;

[0021] (ii) Connect the endpoints of the current crack front boundary line to be expanded and the first and last nodes of the initial expansion vector to form a new crack front;

[0022] (iii) Find the landing point of the current crack front edge on the well wall. Take the intersection point of the plane of the triangular unit closest to the well wall and the plane of the initial crack on the well wall as the landing point of the new crack front edge on the well wall.

[0023] (iv) The Savitzky-Golay filter fitting method is used to smooth the leading edge nodes of the crack to be expanded.

[0024] In one possible implementation, the filtering fitting method can be the Savitzky-Golay filtering fitting method.

[0025] Secondly, this application provides a near-wellbore hydraulic fracture deformation and propagation thermo-fluid-structure interaction simulation system, including: a memory and a processor;

[0026] The memory is used to store computer programs;

[0027] The processor is used to invoke the computer program to execute the method described above.

[0028] Thirdly, this application provides a computer-readable storage medium storing a computer program that, when run on an electronic device, causes the electronic device to perform the method described above.

[0029] Fourthly, this application provides a computer program product, including a computer program that, when run on an electronic device, causes the electronic device to perform the method described above.

[0030] The specific implementation methods of the second to fourth aspects of this application can refer to the implementation methods of the first aspect, and will not be elaborated here.

[0031] Beneficial effects:

[0032] This application employs the three-dimensional displacement discontinuity method (DDM) and the virtual stress method (FSM) to calculate fracture-induced stress and wellbore-induced stress, respectively. The wellbore and fracture are uniformly divided into triangular mesh elements. By mapping specific unknowns to the induced stress, the stress distribution is accurately captured. A boundary stress equilibrium equation is established, superimposing the wellbore-induced stress, fracture-induced stress, and far-field stress to fully characterize the complex stress field near the wellbore. This overcomes the limitations of two-dimensional models, enabling accurate quantification of wellbore stress concentration effects and wellbore-fracture stress shadowing effects in three-dimensional scenarios, thus solving the problem of insufficient stress calculation accuracy in existing models.

[0033] This application constructs a near-wellbore fracture propagation thermo-fluid-structure interaction (TFL) mathematical model, integrating the fluid flow control equations, energy conservation equations, and momentum conservation equations of the matrix and fracture to describe the mass and heat exchange logic between the matrix and fracture. Finite volume and finite element method (FV / FEM) are used to calculate fluid flow, temperature changes, matrix deformation, and coupling exchange processes within the matrix and fracture, incorporating the TFL effect into the overall simulation framework. Within the three-dimensional boundary element method framework, TFL-S complete physics field co-calculation is achieved, accurately capturing the influence of temperature effects and non-uniform fracturing fluid loss on fracture propagation pressure, fracture width, and geometry. It can effectively simulate repeated fracturing scenarios, quantifying the regulatory role of matrix pore pressure and temperature field changes on secondary fracture propagation and reversal, filling the gap in existing models for multi-field coupling simulation.

[0034] Within the framework of online elastic fracture mechanics, combining the maximum circumferential stress criterion (Richard criterion) and the Paris scaling law model, the propagation direction and velocity of fractures are accurately calculated. A dynamic mesh fracture propagation front tracking technique is employed, dynamically updating the fracture geometry through initial propagation vector generation, wellbore wall landing point locking, Savitzky-Golay filtering smoothing, and node merging optimization. This enables full-process dynamic simulation of near-wellbore hydraulic fractures from initiation and redirection to reservoir penetration, accurately reproducing the complex morphology of three-dimensional non-planar fractures and providing a reliable basis for predicting engineering problems such as wellhead high pressure and sand plugging. The technique completely solves the model non-convergence problem caused by fracture front distortion, significantly improving program robustness and enabling stable multi-step propagation simulations of multiple fracture types, providing an efficient and accurate numerical tool for fracturing design.

[0035] This model can flexibly adapt to simulation requirements under different stress conditions (isoaxial, true triaxial), different fracture types (bi-wing, quadri-wing, hexa-wing), and different engineering parameters (initial fracture angle, pumping rate, etc.). The computational logic considers the coupled responses of stress, temperature, and flow fields, and can be transferred to various deep-earth energy development and energy, gas, and carbon storage scenarios. It is not only applicable to fracturing design in tight oil and gas and shale oil and gas reservoirs, but can also be extended to wellbore stability prediction, deep geothermal fracturing stimulation, CO2 geological sealing caprock stability calculation, and oil and gas reservoir integrity assessment, showing broad application prospects. This application can be used for numerical simulation of near-wellbore hydraulic fracture propagation during reservoir stimulation in shale oil and gas and tight oil and gas reservoirs, providing a numerical simulation tool for hydraulic fracturing design in such reservoirs, improving the technological level of shale oil and gas and tight oil and gas reservoir stimulation, and promoting their large-scale and efficient development. Attached Figure Description

[0036] Figure 1 This is a schematic diagram of the tracking and smoothing process of the leading edge of the longitudinal seam near the wellbore in one embodiment of this application;

[0037] Figure 2 This is a schematic diagram of the smoothing treatment and node merging treatment of the leading edge of the radial double-wing longitudinal slot near the wellbore in one embodiment of this application;

[0038] Figure 3 This is a comparison diagram of simulation results with two-dimensional analytical reference solutions and three-dimensional numerical reference solutions in one embodiment of this application;

[0039] Figure 4 Different dimensionless parameters in one embodiment of this application Comparison of three-dimensional non-planar fracture morphology near the wellbore under certain conditions;

[0040] Figure 5 For different A comparison of crack propagation paths predicted by the two-dimensional model and the three-dimensional model of this application under the given conditions;

[0041] Figure 6 This is a three-dimensional expansion morphology diagram of the biplane longitudinal seam under different initial crack angles in one embodiment of this application;

[0042] Figure 7 This is a diagram showing the intersection of the longitudinal fracture of the double-wing structure and the horizontal symmetry plane of the wellbore under different initial fracture angles in one embodiment of this application.

[0043] Figure 8 This is a simulation result of the propagation of multiple longitudinal cracks near the wellbore under isoaxial stress conditions in one embodiment of this application;

[0044] Figure 9 This is a simulation result of the expansion of multiple longitudinal cracks near the wellbore under true triaxial stress conditions in one embodiment of this application. Detailed Implementation

[0045] To enable those skilled in the art to better understand the present application, the technical solution of the present application will be further described in detail below with reference to the embodiments and accompanying drawings.

[0046] Current hydraulic fracturing models neglect the effects of wellbore stress concentration and thermo-fluid-structure interaction on hydraulic fracture propagation, making it difficult to accurately predict the complexity of near-wellbore stress distribution and fracture geometry. To address the simulation challenges of initiation, propagation, and directional evolution of near-wellbore full three-dimensional fractures, this application proposes a thermo-fluid-structure interaction simulation method for near-wellbore hydraulic fracture deformation and propagation. The method employs the three-dimensional displacement discontinuity method (DDM) and the virtual stress method (FSM) to calculate fracture-induced stress and wellbore-induced stress, respectively. Finite volume and finite element methods (FV / FEM) are used to calculate fluid flow, temperature changes, matrix deformation, and fluid mass and heat exchange between the matrix and fracture. Within the framework of linear elastic fracture mechanics, the maximum circumferential stress criterion and the Paris scaling law model are used to calculate the fracture propagation direction and velocity. Combined with dynamic mesh fracture propagation front tracking technology, a three-dimensional non-planar fracture propagation calculation model is formed under wellbore boundary constraints. This model can consider wellbore stress concentration effects, wellbore-fracture stress shadowing effects, temperature effects, and non-uniform filtration of fracturing fluid. It can accurately predict the arbitrary angles of hydraulic fracture propagation and the eventual generation of three-dimensional non-planar fractures during the fracturing process, from initiation near the wellbore, through the near-wellbore region, and deep into the reservoir. This provides an effective simulation tool for efficient fracturing in tight oil and gas and shale oil and gas reservoirs. The method comprehensively considers the heat and mass exchange between hydraulic fractures and the matrix, heat and mass transfer of fracturing fluid in the matrix and fractures, wellbore stress concentration and wellbore-fracture stress interference, and the mechanical processes of three-dimensional fracture propagation and orientation. It provides a reliable tool for predicting and designing fracture initiation pressure, propagation pressure, three-dimensional fracture geometry, and proppant particle size distribution during the stimulation of tight oil and gas and shale oil and gas reservoirs. This method is not only applicable to the fracturing and optimization design of low-permeability oil and gas reservoirs, but also to wellbore stability prediction, deep geothermal fracturing, CO2 geological sealing and caprock stability calculation, and oil and gas storage stability and integrity assessment. It has broad application prospects in deep geothermal energy development and energy, gas and carbon storage.

[0047] The present application will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0048] Example 1:

[0049] This application provides a method for simulating near-wellbore hydraulic fracture deformation and propagation using a thermo-fluid-structure interaction method, including:

[0050] S1. The reservoir rock mass is generalized into two parts: matrix and fracture. A near-wellbore fracture propagation thermal-fluid-solid coupled mathematical model is established to describe the mechanical processes of fluid flow, heat conduction and heat convection, rock mass deformation and fracture propagation within the matrix and fracture.

[0051] S2. Based on the aforementioned thermal-fluid-structure interaction mathematical model, the three-dimensional displacement discontinuous method (DDM) and the virtual stress method (FSM) are used to calculate the fracture-induced stress and wellbore-induced stress, respectively; the finite volume method and the finite element method are used to calculate the fluid flow, temperature change, matrix deformation, and fluid mass and heat exchange between the matrix and the fracture.

[0052] S3. Within the theoretical framework of linear elastic fracture mechanics, and combined with the calculation results of S2, the maximum circumferential stress criterion and the Paris scaling law model are used to calculate the crack propagation direction and propagation speed.

[0053] S4. Employing dynamic mesh fracture propagation leading edge tracking technology, the fracture geometry is updated in real time based on the fracture propagation direction and speed, forming a three-dimensional non-planar fracture propagation calculation model under wellbore boundary constraints, realizing the full-process deformation and propagation simulation of near-wellbore hydraulic fractures from initiation and turning to deep reservoir penetration.

[0054] The following is a detailed description of the near-wellbore fracture propagation thermal-fluid-structure interaction mathematical model (thermal-fluid-structure interaction calculation model) in this application.

[0055] In some embodiments, it is assumed that the injected fluid is compressible clean water, and the filter cake effect is ignored; the reservoir rock mass is generalized into two parts: matrix (rock skeleton without fractures + pores) and fractures (fracture walls and internal space), and described using the embedded discrete fracture method (EDFM). The fluid flow in the fractures and matrix follows the cubic law and Darcy's law, respectively; it is assumed that the reservoir rock is a slightly compressible material, and its thermal conductivity, compressibility, Young's modulus, Poisson's ratio and other parameters do not change with temperature and pressure; the fluid and reservoir rock mass form a non-adiabatic system, and it is assumed that the thermal equilibrium between the pore fluid and the rock skeleton in the characterizing unit can be achieved instantaneously, and the temperature difference between them can be ignored in a single time step; the shear dilatation / shear contraction effect caused by the protrusions on the hydraulic fracture wall (i.e., smooth fractures) is ignored; the changes in the physical parameters of the fluid and rock mass caused by the water-rock reaction are ignored; it is assumed that fracture propagation is a quasi-static process, and follows the linear elastic fracture mechanics theory; this application does not consider the tearing mechanics behavior of the fracture propagation front, that is, it does not consider the propagation of type III or type I-III hydraulic fractures. This application involves four main mechanical processes: fluid flow, heat conduction and convection, rock mass deformation, and crack propagation. Each mechanical process can be quantitatively described by the following equations / criteria, as detailed below:

[0056] a) Fluid flow process, describing the flow law of fracturing fluid in the matrix and fracture, as well as the fluid mass exchange (filtration loss) between the matrix and fracture, can be described by fluid flow control equations and mass exchange equations.

[0057] Based on the law of mass conservation, the governing equations for fluid flow in the matrix and the fracture are expressed as follows:

[0058] ;

[0059] ;

[0060] The viscous flow of fluid in the matrix and the crack can be described as follows:

[0061] ;

[0062] ;

[0063] The equation for fluid mass exchange between the matrix and fractures (the flow of fracturing fluid from the fractures into the matrix pores is called filtration loss) is:

[0064] ;

[0065] b) The heat conduction and heat convection processes describe the heat transfer (conduction + convection) within the rock matrix and fractures, as well as the heat exchange between the matrix and fractures. These processes can be described by the energy conservation equation, the heat flux equation, and the heat exchange equation.

[0066] Based on the transient thermal equilibrium assumption, the energy conservation equation for the matrix portion of the rock mass can be expressed as:

[0067] ;

[0068] ;

[0069] Heat conduction follows Fourier's law; therefore, the heat flux equation in the matrix can be expressed as:

[0070] ;

[0071] Similarly, the energy conservation equation and heat flux for the crack section can be expressed as follows:

[0072] ;

[0073] ;

[0074] The energy exchange equation between the matrix and the crack can be described as follows:

[0075] ;

[0076] The physical quantities appearing in the above formulas have the following meanings: Biot coefficient, dimensionless; It is the Biot modulus, dimensionless; , These are pore pressure and fracture flow pressure, respectively, in Pa; , , , respectively, are the thermal expansion coefficients of the rock mass, matrix, and fluid, in 1 / °C; Fluid density, kg / m³ 3 ; It is a body that adapts to change, without dimensions; Indicates temperature change, °C; , Temperature changes in the matrix and cracks, respectively, in °C; Temperature is expressed in °C. Indicates time, in seconds; Porosity, % , The average flow velocities of the matrix and fractures are respectively, in m / s; , These represent the mass flow sources and sinks within the matrix and fractures, respectively, in kg / s; , These represent energy source sinks within the matrix and fractures, respectively, in J / s; , The flow fluxes into the matrix and fractures, respectively, are m. 3 / s; This represents the number of contact units between the crack and the matrix, and is dimensionless. Represents the gradient; Indicates divergence; , These are the enthalpy fluxes flowing into the matrix and the fracture, respectively, in J / (s). m 2 ); The fluid compressibility coefficient is 1 / Pa. Indicates seam width, in meters (m). The contact area between the matrix and the crack, in meters (m). 2 ; , For matrix permeability scalar and tensor, m 2 ; Viscosity is expressed as cP. Represents pore pressure and fracture flow pressure, in Pa; Enthalpy is expressed in J / kg; , Thermal conductivity of the matrix and crack, respectively, W / (m²). K); , The flow fluxes into the matrix and fractures are respectively, in J / (s) m 2 ); Represents the volume integral domain of the rock mass matrix, m 3 , express The volume of the infinitesimal element on the surface; Denotes the volume integral domain of the crack, m 3 , express The volume of the infinitesimal element on the surface; m represents the contact area between the i-th crack and the matrix. 2 , This represents the area of ​​a micro-element on the contact surface; This represents the crack normal vector, which is dimensionless.

[0077] The equations in a) and b) above can be collectively referred to as the heat and mass transfer control equations.

[0078] c) Rock mass deformation process: describing the deformation law of rock mass in the near-wellbore area under stress, which can be described by mechanical equilibrium equations.

[0079] Among them, the mechanical equilibrium equations include the momentum conservation equation, the reservoir rock mass constitutive equation, and the equations for calculating induced stress and discontinuous displacement.

[0080] For the matrix portion, neglecting the effect of inertial forces, its momentum conservation equation can be described as:

[0081] ;

[0082] in, For surface stress, Pa; For the concentrated force, N; Density of the rock mass, kg / m³ 3 ;

[0083] Treating the reservoir rock mass as a thermally porous elastic material, its constitutive equation can be expressed as:

[0084] ;

[0085] in, This represents the Cauchy total stress tensor, Pa. The effective stress tensor is expressed in Pa. Represents the rock mass drainage elastic tensor, in Pa; For the rock mass strain tensor, which is dimensionless and corresponds to the overall deformation of the rock mass.

[0086] This application employs DDM and FSM methods to calculate the induced stress field and fracture deformation of the fracture and wellbore. In a three-dimensional context, the fracture and wellbore can be uniformly divided into a series of triangular meshes. Each fracture triangular computational unit (i.e., the smallest computational unit of the fracture, hereinafter referred to as the fracture unit) contains one normal displacement discontinuity and two shear displacement discontinuities, while each wellbore triangular computational unit (i.e., the smallest computational unit of the wellbore, hereinafter referred to as the wellbore unit) contains one normal virtual stress and two shear virtual stresses; that is, each triangular computational unit contains three unknowns. For each boundary element computational unit (triangular computational unit), its mechanical equilibrium equation can be obtained by the total stress equaling the initial stress plus the fracture and wellbore induced stresses. The case of fracture closure is not considered during fracturing; therefore, the stress equilibrium equation at the boundary can be expressed as:

[0087] ;

[0088] ;

[0089] The induced stress acting on the wellbore boundary calculation unit can be written as:

[0090] ;

[0091] The induced stress acting on the calculation element at the crack boundary can be written as:

[0092] ;

[0093] The physical quantities appearing in the above formulas have the following meanings: The pressure of the fluid inside the wellbore, in Pa; , These are the far-field stress tensor and the induced stress tensor, respectively, in Pa; This represents the normal vector of a crack or well wall surface, and has no dimension. These represent the shear-induced stress along the direction of orientation, the shear-induced stress along the dip direction, and the shear-induced stress along the normal direction, respectively, in Pa;

[0094] These are the virtual stress values ​​(in Pa) of the wellbore element along the two shear directions and the normal direction, respectively. Let m be the discontinuous displacement of the crack element along the two shear directions and the normal direction, respectively.

[0095] These represent the shear-induced stress in the two directions and the normal-direction induced stress of the boundary element, respectively. Indicates the first The shear-induced stress at the geometric center of each wellbore unit, Pa; Indicates the first The inclined shear-induced stress at the geometric center of each wellbore unit, Pa; Indicates the first Normal induced stress at the geometric center of each wellbore unit, Pa; Indicates the first The shear-induced stress at the geometric center of each crack element, Pa; Indicates the first The inclined shear-induced stress at the geometric center of each crack element, Pa; Indicates the first Normal induced stress at the geometric center of a crack element, Pa; The first The virtual stress values, in Pa, of a wellbore element along two shear directions and the normal direction; These represent the discontinuous displacement of the crack element along the two shear directions and the normal direction, respectively, in meters, corresponding to the shear deformation and crack width variation (normal deformation) of the crack. , The first The unit virtual stress of the wellbore element along the s direction in the th... The influence coefficient of the magnitude of the induced normal stress at each wellbore unit is dimensionless and is the i-th The discontinuity of the unit displacement of the crack element along the s direction in the th... The influence coefficient of the magnitude of induced normal stress at each crack element, Pa / m; , The first The unit virtual stress of the wellbore element along the t direction at the 1st... The influence coefficient of the magnitude of the induced normal stress at each wellbore unit is dimensionless and is the i-th The discontinuity of the unit displacement of the crack element along the t direction in the th... The influence coefficient of the magnitude of induced normal stress at each crack element, Pa / m; , The first The unit virtual stress of the wellbore element along the normal direction in the first... The influence coefficient of the magnitude of the induced normal stress at each wellbore unit is dimensionless and is the i-th The discontinuity of unit displacement along the normal direction of each crack element in the first... The influence coefficient of the magnitude of induced normal stress at each crack element, Pa / m; , The first The unit virtual stress of the wellbore element along the s direction in the th... The influence coefficient of the magnitude of the induced dip shear stress at each wellbore unit is dimensionless and of the th order. The discontinuity of unit displacement along the s direction in each crack element is within the crack element. The influence coefficient of the magnitude of the induced shear stress at the location, Pa / m; , The first The unit virtual stress of the wellbore element along the t direction at the 1st... The influence coefficient of the magnitude of the induced dip shear stress at each wellbore unit is dimensionless and of the th order. The discontinuity of the unit displacement of the crack element along the t direction in the th... The influence coefficient of the magnitude of the induced shear stress at each crack element, Pa / m; , The first The unit virtual stress of the wellbore element along the normal direction in the first... The influence coefficient of the magnitude of the induced dip shear stress at each wellbore unit is dimensionless and of the th order. The discontinuity of unit displacement along the normal direction of each crack element in the first... The influence coefficient of the magnitude of the induced shear stress at each crack element, Pa / m; , The first The unit virtual stress of a wellbore element along the s direction in the wellbore element The influence coefficient of the magnitude of the induced shear stress at the location is dimensionless and of the th order. The discontinuity of unit displacement along the s direction in each crack element is within the crack element. The influence coefficient of the magnitude of the induced shear stress at the location, Pa / m; , The first The unit virtual stress of the wellbore element along the t direction at the 1st... The influence coefficient of the magnitude of the induced shear stress at each well unit is dimensionless and of the th order. The discontinuity of the unit displacement of the crack element along the t direction in the th... The influence coefficient of the magnitude of the induced shear stress at each crack element, Pa / m; , The first The unit virtual stress of the wellbore element along the normal direction in the first... The influence coefficient of the magnitude of the induced shear stress at each well unit is dimensionless and of the th order. The discontinuity of unit displacement along the normal direction of each crack element in the first... The influence coefficient of the magnitude of the induced shear stress at each crack element, Pa / m. , These represent the number of virtual stress (FSM) elements associated with the wellbore and the number of displacement discontinuity (DDM) elements associated with the fracture, respectively, and are dimensionless.

[0096] d) Crack propagation process, describing the crack initiation conditions, propagation direction and propagation length, described by crack propagation prediction equation, propagation length and propagation direction calculation equation.

[0097] This application employs a three-dimensional crack propagation criterion to predict crack propagation. This criterion assumes that when the composite stress concentration factor at the crack tip... Equal to the fracture toughness of rocks When the crack propagates, the crack propagation prediction equation is:

[0098] ;

[0099] in, For rock fracture toughness, MPa ; , , These are the stress intensity factors for opening, slip, and tearing types, respectively, in MPa. ; The composite stress intensity factor is MPa. ; The crack turning angle is expressed in degrees.

[0100] Considering that fracture propagation during hydraulic fracturing generally falls within the quasi-static propagation category, and to avoid the fracture propagation front exhibiting an intermittent sawtooth morphology, it is assumed that any point on the fracture front satisfies the following condition. , The Paris scaling law is then used to give the propagation length at any point on the crack front within a certain time step:

[0101] ;

[0102] in, Represents empirical coefficients, which are dimensionless; Hereinafter, m, represent the crack element propagation length and the maximum allowable crack element propagation length, respectively.

[0103] The crack propagation direction was calculated according to the maximum circumferential stress criterion. Comparison of calculations revealed that the crack orientation angles obtained using the Schöllmann and Richard criteria were very close. For simplicity in implementation, the principal stress calculations involved in the Schöllmann criterion can be avoided, and the crack propagation direction, i.e., the crack orientation angle, can be calculated directly using the Richard criterion based on the three fundamental stress intensity factors.

[0104] ;

[0105] Among them, when When, the crack turning angle is taken as "+", when and At that time, the crack turning angle is taken as " ".

[0106] like Figure 1 As shown, in some embodiments, a dynamic mesh crack propagation leading edge tracking technology (three-dimensional non-planar crack propagation leading edge tracking algorithm) is employed to update the crack geometry in real time based on the crack propagation direction and propagation speed, including:

[0107] A three-dimensional non-planar crack propagation leading edge tracking algorithm, including:

[0108] (i) Based on the crack propagation length within the current time step ( ) and expansion direction ( ), thus obtaining a series of initial expansion vectors ,See Figure 1 (a);

[0109] (ii) Connect the endpoints of the current crack front boundary line to be expanded and the first and last nodes of the initial expansion vector to form a new crack front, such as... Figure 1 (b);

[0110] (iii) Find the point where the leading edge of the current crack to be expanded lands on the wellbore wall. This is the intersection of the plane containing the triangular element closest to the wellbore at the current crack leading edge and the plane containing the initial crack at the wellbore wall. This intersection is considered the new landing point of the crack leading edge on the wellbore wall (i.e., the intersection with the wellbore). See [link to relevant documentation]. Figure 1 (c);

[0111] (iv) The Savitzky-Golay filter fitting method is used to smooth the leading edge nodes of the crack to be expanded.

[0112] For example, it is required that the two points of impact on the upper and lower sides of the longitudinal crack after smoothing should be on the same plane as the initial crack, such as... Figure 1 (d) At the same time, the distance between the landing points on the upper and lower sides and the landing point of the previous time step, as well as the distance between the adjacent nodes at the leading edge of the crack to be expanded, should meet the set constraints. Otherwise, the landing points of the previous and next time steps should be merged (e.g., the landing points of the previous time step should be merged into the landing point of the current time step), or the adjacent nodes should be merged.

[0113] The constraints set can be... ,in The distance between adjacent nodes. This represents the maximum allowable length of the crack element.

[0114] The main purpose of smoothing the crack propagation front is to avoid intermittent serrated or other deformed crack fronts that could affect the robustness of the model. Figure 2This paper demonstrates the near-wellbore 3D nonplanar fracture morphology obtained after smoothing and node merging of the fracture propagation front. Simulation results show that the unsmoothed fracture propagation front is prone to jagged edges and inconsistent morphology, leading to non-convergence and program crashes in the fracture propagation model. After smoothing the fracture front, the stability of the fracture propagation model is significantly improved. However, due to the frequent occurrence of very close local nodes during propagation, mesh inversion and overlap can easily occur, forming locally deformed fracture elements, which can again cause non-convergence and crashes in the fracture propagation program. To further improve the stability of the program, conditions will be met... The adjacent nodes are merged, and the robustness of the crack propagation procedure is significantly enhanced after the merging process. Figure 2 During the fracturing process, the crack expands from its initial morphology through 41 steps, including initiation, propagation, and reversal, eventually forming a three-dimensional non-planar crack morphology. Throughout the simulation, no program crashes occurred due to distorted meshes, indicating that the Savitzky-Golay filtering and fitting method, along with node merging, enhances the robustness of the crack propagation model and is an effective technique for 3D crack propagation front tracking.

[0115] The following section verifies the near-wellbore three-dimensional non-planar crack propagation calculation model obtained in the embodiments of this application. The main focus is on verifying whether the model can accurately predict the impact of wellbore-induced stress concentration on the global stress field and the three-dimensional non-planar crack propagation process.

[0116] a) Verification of the calculation of the stress field induced by the wellbore-fracture:

[0117] The stress distribution near a wellbore in an infinitely large space was examined, referring to the parameters given by Rahman et al. (2002): Young's modulus E = 20 GPa, Poisson's ratio ν = 0.2, wellbore length L = 8 m, wellbore radius R = 1 m, and initial fracture length. =1m, fluid pressure inside the wellbore =2.5MPa, the double-wing fracture is divided into 202 calculation units, and the wellbore is divided into 1952 calculation units.

[0118] It should be noted that for a two-dimensional analytical reference solution, the problem it considers can be degenerated from the problem of an infinitely long wellbore and an infinitely long slot height existing in a three-dimensional infinite space. Therefore, for such... Figure 3 The stress level of the solution to the three-dimensional problem shown is clearly lower than that of the two-dimensional analytical solution. (From...) Figure 3 The comparison results show that the model proposed in this application has a good agreement with the three-dimensional numerical reference solution given by Rahman et al. (2002) (including dimensionless stress). , , However, all results are lower than the two-dimensional analytical reference solution, which is consistent with the conclusion that there is a difference between the actual physical problem corresponding to the two-dimensional analytical solution and the simulated three-dimensional problem. This simulation comparison verifies that the model can accurately simulate the wellbore-fracture induced stress field, thus providing the necessary conditions for fracture propagation calculations.

[0119] b) Validation of the near-wellbore longitudinal fracture propagation model:

[0120] Hydraulic fracture propagation is mainly influenced by factors such as in-situ stress, pumping rate, fracturing fluid viscosity, and fracture toughness, while near-wellbore fracture direction and propagation, as well as fracture morphology, are primarily determined by dimensionless parameters. Control (Mogilevskaya et al., 2000):

[0121] ;

[0122] in, , These are the maximum and minimum horizontal principal stresses, respectively, in Pa;

[0123] Mogilevskaya et al. (2000) established a two-dimensional near-wellbore fracture propagation model, and based on this model, they obtained different... Two-dimensional crack propagation morphology under certain conditions. To verify the reliability of this application, three sets of computational parameters were selected for the three-dimensional model, corresponding to different dimensionless parameters, referring to the work of Mogilevskaya et al. Value: (i) : MPa =2.0 MPa ;in, The minimum vertical principal stress is given in Pa; (ii) : MPa =2.0MPa (iii) : MPa =2.0 MPa The three coordinate axes are the principal stress directions, satisfying... , , ;like Figure 4 (a) Other common parameters include: wellbore length L = 20 m, wellbore radius R = 1.0 m, and initial fracture length. =0.5 m, initial crack angle (the angle between the crack and the horizontal plane) =81°, Young's modulus E=20 GPa, Poisson's ratio =0.2, matrix pore pressure =1.5 MPa, fluid pressure inside the wellbore =2.0 MPa, fluid injection rate =0.01 m 3 / s, viscosity =1.0 cP, maximum crack propagation distance =0.3 m.

[0124] Figure 4 (b)-(d) show different dimensionless parameters Crack morphology under certain conditions. When At that time, because the horizontal stress difference was zero, the crack did not change direction and continued to grow in the initial direction, ultimately resulting in a planar crack; when When the value increases to 0.25, the horizontal stress difference is 0.5 MPa, and the cracks begin to propagate and eventually form a complex curved surface morphology; further increases... When the value reaches 0.5, the hydraulic cracks become more pronounced due to the increase in horizontal stress difference to 2 MPa, and complex three-dimensional non-planar cracks are formed.

[0125] Figure 5 This paper presents a comparison of the results calculated using the three-dimensional model proposed in this application with the two-dimensional model provided by Mogilevskaya et al. (2000). It should be noted that, similar to the previously mentioned cases, the calculation results of the two-dimensional and three-dimensional models will inevitably differ. The fundamental reason for this is that the physical problems they solve are not equivalent. Considering the third dimension, cracks can not only change direction and propagate horizontally but also vertically; the vertical direction change completely undermines the fundamental condition for transforming a three-dimensional problem into a two-dimensional problem. The comparison results show that, although there are certain differences between the simulation results of the three-dimensional model and the two-dimensional model, the crack propagation path and final morphology exhibit a consistent pattern across the observed cross-section, reflecting the reliability of the model proposed in this application in predicting the propagation of three-dimensional non-planar cracks near the wellbore.

[0126] The near-wellbore three-dimensional non-planar fracture propagation calculation model established based on the embodiments of this application can be used to study the influence of the initial fracture angle on the fracture turning distance and the three-dimensional fracture geometry. The initial fracture angle has a significant impact on the fracture initiation and propagation pressure. The larger the angle between the fracture and the direction of the minimum principal stress, the higher the geostress the fracture needs to overcome during its turning perpendicular to the direction of the minimum principal stress. The flow friction of the fracturing fluid will also increase accordingly. Under the same conditions, the lower the net pressure inside the fracture, the smaller the fracture width. This makes fracture initiation and propagation more difficult, increases the risk of sand plugging, and reduces the proppant placement effect in hydraulic fractures. For fractured-vuggy reservoirs, it is usually necessary to scientifically design parameters such as the initial fracture angle (perforation azimuth), pumping rate, and viscosity to achieve better artificial fractures connecting randomly distributed karst bodies within the reservoir. Therefore, studying the influence of the initial fracture angle on the fracture turning distance and the three-dimensional fracture geometry has important practical application value.

[0127] This application uses the three-dimensional displacement discontinuity method and virtual stress method to solve the wellbore-fracture induced stress field. Combined with dynamic mesh fracture propagation leading-edge tracking technology, it forms a three-dimensional non-planar fracture propagation calculation model under wellbore boundary constraints. This model can accurately simulate the propagation process of near-wellbore double-wing, four-wing, and six-wing fractures and the final three-dimensional non-planar fracture geometry (see rendering). Figure 6 and Figure 7 It can provide an effective tool for fracturing design, proppant design and scheme optimization for oil and gas reservoir stimulation.

[0128] Using the method described in this application, numerical simulations were performed on the propagation process of radial fractures near the wellbore in double-wing, four-wing, and six-wing configurations. Isoaxial stress conditions were investigated. MPa) and under true triaxial stress conditions ( MPa MPa MPa), the three-dimensional morphological evolution of multiple longitudinal fractures near the wellbore under the control of wellbore-fracture induced stress, the initial fractures are uniformly distributed radially on the horizontal symmetry plane of the wellbore (see... Figure 8 The main parameters of the model are: wellbore length L = 10m, wellbore radius R = 0.5m, and initial fracture length. =0.25m, Young's modulus E=30GPa, Poisson's ratio =0.2, matrix pore pressure =2.0MPa, fluid pressure inside the wellbore =2.5MPa, fluid injection rate =0.005m 3 / s (injection rate of a single-wing crack), viscosity =1.0cP MPa MPa MPa =2.0MPa The double-wing fracture profile is divided into 76 calculation units, the wellbore profile is divided into 1197 calculation units, and the maximum fracture propagation distance is... =0.15 m.

[0129] Figure 8 The simulation results of competitive propagation of multiple fractures near the wellbore under isoaxial stress conditions (two-wing, four-wing, and six-wing) are presented, including the three-dimensional geometry of the fractures and the fracture width distribution at different locations. In the two-wing fracture case, all fractures propagate along their initial direction without any change in direction, ultimately forming symmetrical three-dimensional planar artificial fractures, such as... Figure 8 (a). See also Figure 8 (b) In the case of four longitudinal cracks, mutual interference occurred during crack propagation, causing each crack to slightly deviate. Overall, both the crack geometry and crack width distribution showed good consistency and spatial symmetry among the longitudinal cracks. See [link to relevant documentation] Figure 8 (c) For the six-wing longitudinal fracture, there is a relatively obvious competitive propagation among the fractures. The middle wing on the right side is strongly suppressed after propagating for 3 steps, and then the propagation step size of this fracture becomes smaller and undergoes obvious reversal propagation. The reversal of this fracture breaks the initial symmetrical stress condition, and the upper wing on the right and the lower wing on the left subsequently exhibit a three-dimensional non-planar propagation morphology. The lower wing on the right and the upper wing on the left maintain a good competitive advantage, generating relatively large fractures with small fracture surface curvature, and the overall fracture morphology is close to planar. As can be seen from the three comparative examples, when multiple fractures appear on the wellbore wall, even under the condition of equal stress in all directions, some fractures will be significantly suppressed and exhibit differential propagation. In fact, for the case of six-wing fractures, the main reason for differential propagation is that the induced stress in the wellbore has a sinusoidal fluctuation form around the well, which makes the stress state of each fracture somewhat different, while the first two cases just ensure the consistency of the stress state of each fracture.

[0130] Compared to isoaxial stress, the propagation of multiple longitudinal fractures near the wellbore is more complex under true triaxial stress. For example... Figure 9 As shown, under the influence of wellbore-fracture induced stress, the longitudinal fractures of the two wings gradually change from the initial direction of the fracture to a direction perpendicular to the far-field minimum principal stress. The two wing fractures exhibit good symmetry in both geometry and width distribution. Figure 9 (a); Similarly, since the initial positions of the longitudinal fractures near the four wings of the wellbore have the same perimeter stress conditions, the final fracture propagation morphology also exhibits good symmetry, such as... Figure 9(b) differs in that, due to the synchronous competitive expansion of the cracks on both sides, the two cracks repel each other under the influence of stress shadow and continuously extend outwards, resulting in a significant reduction in the turning radius and turning distance (see...). Figure 9 (Top view); For the longitudinal fractures of the six wings near the wellbore, due to the true triaxial stress condition and the sinusoidal variation of the wellbore stress, each fracture must continuously adjust its propagation direction in the changing stress field, finding the local optimal propagation direction between the wellbore-fracture induced stress and the far-field stress. From Figure 9 (c) It can be observed that there are significant differences in the morphology, width, size, and direction of fractures at different locations, reflecting the complexity of near-wellbore stress conditions and the competitive propagation patterns of multiple fractures. In fracturing design, parameters such as perforation location, perforation spacing, and phase angle should be rationally selected to ensure uniform propagation of each fracture as much as possible, thereby improving the reservoir stimulation effect of hydraulic fracturing.

[0131] Example 2:

[0132] This embodiment provides a near-wellbore hydraulic fracture deformation and propagation thermo-fluid-structure interaction simulation system, including: a memory and a processor;

[0133] The memory is used to store computer programs;

[0134] The processor is configured to invoke the computer program to execute the method as described in Embodiment 1.

[0135] Example 3:

[0136] This embodiment provides a computer-readable storage medium storing a computer program. When the computer program is run on an electronic device, it causes the electronic device to perform the method described in Embodiment 1.

[0137] Example 4:

[0138] This embodiment provides a computer program product, including a computer program that, when run on an electronic device, causes the electronic device to perform the method described in Embodiment 1.

[0139] The specific implementation of the system, electronic device, computer-readable storage medium, and computer program product provided in this application can be referred to the specific embodiments of the above methods, and will not be repeated here.

[0140] Obviously, those skilled in the art should understand that the various units or steps of this application described above can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using computer-executable program code, thereby storing them in a storage device for execution by a computing device, or fabricating them separately as individual integrated circuit modules, or fabricating multiple modules or steps into a single integrated circuit module. Thus, this application is not limited to any particular combination of hardware and software.

[0141] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A thermo-fluid-structure interaction simulation method for near-wellbore hydraulic fracture deformation and propagation, characterized in that, include: S1. The reservoir rock mass is generalized into two parts: matrix and fracture. A near-wellbore fracture propagation thermal-fluid-solid coupled mathematical model is established to describe the mechanical processes of fluid flow, heat conduction and heat convection, rock mass deformation and fracture propagation within the matrix and fracture. S2. Based on the aforementioned thermo-fluid-structure interaction mathematical model, the three-dimensional displacement discontinuity method (DDM) and the virtual stress method (FSM) are used to calculate the fracture-induced stress and wellbore-induced stress, respectively. The finite volume and finite element methods were used to calculate fluid flow, temperature changes, matrix deformation, and fluid mass and heat exchange between the matrix and the fracture. S3. Within the theoretical framework of linear elastic fracture mechanics, and combined with the calculation results of S2, the maximum circumferential stress criterion and the Paris scaling law model are used to calculate the crack propagation direction and propagation speed. S4. Employing dynamic mesh crack propagation leading edge tracking technology, the crack geometry is updated in real time based on the crack propagation direction and speed, including: (i) Based on the crack propagation length and propagation direction in the current time step, a series of initial propagation vectors are obtained; (ii) Connect the endpoints of the current crack front boundary line to be expanded and the first and last nodes of the initial expansion vector to form a new crack front; (iii) Find the landing point of the current crack front edge on the well wall. Take the intersection point of the plane of the triangular unit closest to the well wall and the plane of the initial crack on the well wall as the landing point of the new crack front edge on the well wall. (iv) The leading edge nodes of the current crack to be expanded are smoothed using the filtering fitting method; Thus, a three-dimensional non-planar fracture propagation calculation model under wellbore boundary constraints was formed, realizing the full-process deformation and propagation simulation of near-wellbore hydraulic fractures from initiation and deflection to penetration into the reservoir.

2. The method according to claim 1, characterized in that, The thermal-fluid-solid coupled mathematical model includes: fluid flow control equations in the matrix and fractures, fluid mass exchange equations between the matrix and fractures, energy conservation equations for the matrix portion of the rock mass, heat flux equations in the matrix, energy conservation equations and heat flux equations for the fracture portion, energy exchange equations between the matrix and fractures, momentum conservation equations for the matrix portion, constitutive equations for the reservoir rock mass, calculation equations for induced stress and displacement discontinuities, fracture propagation prediction equations, and calculation equations for fracture propagation length and direction.

3. The method according to claim 1, characterized in that, The three-dimensional displacement discontinuity method (DDM) and the virtual stress method (FSM) are used to calculate the fracture-induced stress and wellbore-induced stress, respectively. This includes: uniformly dividing the fracture and wellbore into a series of triangular meshes; each fracture triangular calculation element contains one normal displacement discontinuity and two shear direction displacement discontinuities; the induced stress of each fracture triangular calculation element is obtained through the mapping relationship between the displacement discontinuities and the induced stress; each wellbore triangular calculation element contains one normal virtual stress and two shear direction virtual stresses; the induced stress of each wellbore triangular calculation element is obtained through the mapping relationship between the virtual stress and the induced stress.

4. The method according to claim 3, characterized in that, The rock mass deformation includes matrix deformation and crack deformation. The crack deformation includes shear deformation and crack width variation, which are described by the discontinuity of normal displacement and shear direction displacement of the crack triangle calculation unit, respectively.

5. A near-wellbore hydraulic fracture deformation and propagation thermo-fluid-structure interaction simulation system, characterized in that, include: Memory and processor; The memory is used to store computer programs; The processor is configured to invoke the computer program to perform the method as described in any one of claims 1 to 4.

6. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed on an electronic device, causes the electronic device to perform the method as described in any one of claims 1 to 4.

7. A computer program product, comprising a computer program, characterized in that, When the computer program is run on an electronic device, it causes the electronic device to perform the method as described in any one of claims 1 to 4.

Images