A numerical simulation method for fracture development process of faulted fracture in a faulted combination area

Abstract

The application discloses a numerical simulation method for fracture and fissure development process in a fault-reservoir combination area, and belongs to the technical field of geothermal resource development and utilization, and comprises the following steps: a heat-flow coupling model of a stratum three-dimensional model and a fault-reservoir coupling model of the fault-reservoir combination area are constructed; the models have data exchange channels; a multi-physical field simulation is carried out on the stratum three-dimensional model; and meanwhile, the fault-reservoir coupling model is used to simulate the fracture and fissure development process in the fault-reservoir combination area. The application can accurately simulate the heat transfer process of geothermal heat in the fault-reservoir combination area, judge the geothermal heat transfer mode, and then simulate the expansion behavior of the fault and the evolution process of the fissure under different geothermal heat transfer modes in combination with relevant algorithm models. The application can accurately simulate the fracture and fissure development process, obtain the influence of complex structures such as the fault-reservoir combination area on geothermal heat transfer and geothermal fluid migration, provide a scientific basis for geothermal exploitation and structure safety evaluation, and has high practical value and popularization prospect.

Description

Technical Field

[0001] This invention relates to the field of geothermal resource development and utilization technology, specifically to a numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone. Background Technology

[0002] The characteristics of fractures and cracks significantly alter the migration and heat transfer patterns of oil, gas, and geothermal fluids. Accurately simulating the development process of fractures and cracks is an important basis for the genesis and development of geothermal resources.

[0003] Current research on geothermal geology is weak, relying mainly on electromagnetic characterization of fractures, borehole temperature measurement, and geochemical tracing to determine the characteristics of geothermal fluid migration. Among these methods, electromagnetic energy can characterize basic information such as the distribution and dip angle of fractures, geochemical parameters can determine the opening and closing of fractures, and borehole temperature measurement can analyze conductive and convective geothermal resources.

[0004] In regions with complex geological structures, fracture development is complex, and the configuration relationship between faults and reservoirs plays a crucial role in geothermal fluid seepage, heat transfer, and resource allocation. Changes in geothermal resources, in turn, affect the fracture development process. In existing technologies, Ra and Pe numbers are used to determine geothermal heat transfer modes, and then the fracture development process of different geothermal heat transfer modes is further investigated. However, there are the following shortcomings: (1) The difference between fracture and reservoir matrix is ​​not considered: the Ra formula is only applicable to homogeneous porous media and cannot reflect the fracture-reservoir matrix coupling effect; (2) The fracture dip angle is ignored: the fracture orientation has a significant impact on fluid migration, but the existing criteria do not consider this. As a result, the accuracy of geothermal heat transfer mode identification is limited, which in turn reduces the simulation accuracy of fracture development process. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone, which can effectively solve the above-mentioned problems.

[0006] The technical solution adopted in this invention is as follows:

[0007] This invention provides a numerical simulation method for the development process of fractures in a fault-reservoir combination zone, comprising the following steps:

[0008] Step S1: Construct a three-dimensional stratigraphic model of the study area; the three-dimensional stratigraphic model has several fault-reservoir combination zones; the fault-reservoir combination zone is a region with a fault in the reservoir and several fractures distributed in the vicinity of the fault.

[0009] Step S2: Construct the thermal-fluid coupling model of the three-dimensional formation model and the fault-reservoir coupling model of the fault-reservoir combination zone; the fault-reservoir coupling model is a fracture-fracture-reservoir coupling model.

[0010] Step S3, setting numerical simulation initial parameters and boundary conditions; the numerical simulation initial parameters include initial data of multi-physical fields of the thermal-flow coupling model and initial state parameters of fractured fracture reservoirs of each of the fracture-reservoir coupling models;

[0011] Step S4, establishing a data exchange channel between the thermal-flow coupling model and each of the fracture-reservoir coupling models, through which data exchange and feedback adjustment between multi-physical field data solved by the thermal-flow coupling model and state parameters of fractured fracture reservoirs solved by the fracture-reservoir coupling models are realized, and then the thermal-flow coupling model is used to simulate multi-physical fields of the stratum three-dimensional model, and the fracture-reservoir coupling model is used to simulate the development process of fractures in the fracture-reservoir combination area.

[0012] Further, step S1 is specifically:

[0013] Collecting geological data of the research area, including: geological exploration reports, drilling data, seismic exploration data and hydrogeological test data of the research area;

[0014] Analyzing and processing the geological data of the research area to determine the stratum structure parameters of the research area;

[0015] According to the stratum structure parameters of the research area, a modeling tool is used to construct the stratum three-dimensional model.

[0016] Further, the stratum structure parameters include lithology, thickness, distribution range, fracture position, fracture trend, fracture dip angle and fracture distribution data of the stratum;

[0017] The modeling tool used to construct the stratum three-dimensional model is specifically: a finite element method is used to construct the stratum three-dimensional model, and an augmented shape function combining a Heaviside discontinuous function and a crack tip term is used to describe the distribution of crack surface temperature and pressure at the crack tip, and local grid refinement processing is performed in the fracture-reservoir combination area.

[0018] Further, the thermal-flow coupling model is a model used to simulate multi-physical field data in the stratum three-dimensional model; wherein the multi-physical field data includes temperature field, pressure field and seepage velocity field.

[0019] Further, step S4 is specifically:

[0020] Step S41, according to the multi-physical field data of the stratum three-dimensional model at the previous time step and the fracture, crack and reservoir state parameters of each of the fracture-reservoir combination areas, solving the thermal-flow coupling model to obtain the multi-physical field data of the stratum three-dimensional model at the current time step;

[0021] Step S42, the fault-reservoir coupling model determines the multi-physical field data of the fault-reservoir combination zone at the current time step according to the multi-physical field data of the stratum three-dimensional model at the current time step, and determines the temperature change and the fracture pore pressure change of the fracture in the fault-reservoir combination zone according to the multi-physical field data of the fault-reservoir combination zone at the current time step; according to the temperature change and the fracture pore pressure change of the fracture, the fracture permeability and the fracture opening degree of the fracture at the current time step are obtained, and then the fracture state parameters are updated;

[0022] Step S43, based on the updated fracture state parameters, return to step S41, and obtain the multi-physical field data of the stratum three-dimensional model at the next time step by solving the heat-flow coupling model, so that the steps S41 and S42 are continuously cycled, and the fracture evolution feedback mechanism is realized in the process of simulating the multi-physical field data of the stratum three-dimensional model;

[0023] Step S44, in the process of simulating the multi-physical field data based on the fracture evolution feedback mechanism, at the selected time step, the geothermal heat transfer mode of the fault-reservoir combination zone at the current time step is identified according to the fault, fracture and reservoir state parameters;

[0024] According to the identified geothermal heat transfer mode of the fault-reservoir combination zone, the fracture and fracture development process is simulated, the fault, fracture and reservoir state parameters are simulated, and the fault, fracture and reservoir state parameters simulated are used to update the fault, fracture and reservoir state parameters in step S41 in real time.

[0025] Further, in step S42, the fracture permeability and the fracture opening degree at the current time step are obtained by using formula (1) and formula (2):

[0026] (1)

[0027] (2)

[0028] Wherein: and are the initial value of the fracture permeability and the initial value of the fracture opening degree at the current time step, respectively; , , and are the first temperature sensitivity coefficient, the second temperature sensitivity coefficient, the first pressure sensitivity coefficient and the second pressure sensitivity coefficient, respectively; and are the current temperature and the current pore pressure, respectively; and are the initial temperature and the initial pore pressure, respectively.

[0029] Further, the geothermal heat transfer mode of the fault-reservoir combination area at the current time step is identified according to the fault, fracture and reservoir state parameters, specifically:

[0030] The criterion shown in formula (3) is established by comprehensively considering the fault state parameters, the fracture state parameters and the coupling relationship between the fault and the reservoir matrix Formula:

[0031] (3)

[0032] , 0 ≤1

[0033] Wherein, g is the acceleration of gravity; β is the thermal expansion coefficient of the geothermal fluid of the fault-reservoir combination area, ΔT is the temperature difference between the upper and lower ends of the fault; H is the extension height of the fault; k f is the fracture permeability; k m is the thermal reservoir permeability; θ is the fault dip angle; μ is the kinematic viscosity of the geothermal fluid; λ is the coupling coefficient between the fault and the reservoir matrix; ɑ T is the thermal diffusion coefficient; ɑ is the empirical index; x f is the fracture heat conduction enhancement coefficient, η is the heat conduction gain factor, L f is the fracture spacing, is the fracture opening; wherein, for the case that there are multiple fractures distributed in the vicinity of the fault, the fracture spacing L f is the sum of the spacings between the fractures; the fracture opening is the sum of the openings of the fractures;

[0034] The criterion formula is used to obtain the criterion of the fault-reservoir combination area at the current time step; if the criterion ≤ an empirical value, it is identified that the geothermal heat transfer mode of the fault-reservoir combination area at the current time step is the conduction type geothermal heat transfer type; if the criterion is greater than the empirical value, it is identified that the geothermal heat transfer mode of the fault-reservoir combination area at the current time step is the convection type geothermal heat transfer type.

[0035] Further, the fracture and fissure development process is simulated according to the identified geothermal heat transfer mode of the fault-reservoir combination area, specifically:

[0036] Step A1, if the geothermal heat transfer mode is the convection type geothermal heat transfer type, the convection of the thermal fluid in the fault channel is enhanced, the local temperature gradient and the fluid pressure of the fault surface are both increased, and the thermal stress σ T and the pore pressure stress σ PThe coupling superposition effect of the pore pressure and thermal expansion is taken into account, and the effective normal stress in the fracture dip direction θ is obtained by formula (4) :

[0037] (4)

[0038] wherein: is the normal total stress in the fracture dip direction θ; is the pore pressure coupling coefficient; p is the current pore pressure; p0 is the initial pore pressure; is the thermal expansion coefficient; E is the elastic modulus; T is the current temperature; T0 is the initial temperature;

[0039] Step A2, the shear stress in the fracture dip direction θ is obtained by formula (5) :

[0040] (5)

[0041] wherein: is the in-plane shear stress component; is the vertical direction shear stress component;

[0042] Step A3, it is judged whether the sliding condition of formula (6) is established:

[0043] (6)

[0044] wherein: is the shear strength that the fracture surface still has under zero normal stress; is the internal friction angle;

[0045] Step A4, if the sliding condition of formula (6) is met, the fracture enters the shear slip stage, and if the accumulated displacement δ s exceeds the threshold value δ cr , the fracture surface expands, the adjacent cracks of the fracture extend along the fault direction to form secondary cracks, the expansion rate of the secondary cracks satisfies equation (7), the crack expansion direction is consistent with the direction of the maximum principal stress, a heat-conducting and flow-conducting channel is formed, the overall permeability and heat conductivity of the fracture zone are enhanced, and continuous deformation and thermal drive slip of the fracture morphology are triggered:

[0046] (7)

[0047] wherein: L f is the crack spacing; t is time; k v is the crack growth coefficient; m is an empirical index.

[0048] Further, the fracture crack development process is simulated according to the identified geothermal heat transfer mode of the fault-reservoir combination zone, specifically:

[0049] Step B1, if the geothermal heat transfer mode is the conduction type of geothermal heat transfer, in the fault-reservoir combination area, the geothermal energy diffuses through the conduction of the reservoir matrix, the local temperature gradient is high, which causes the thermal expansion and contraction strain in the reservoir matrix, and the reservoir matrix strain value is obtained by using formula (8) :

[0050] (8)

[0051] Wherein: is the first temperature sensitive coefficient; T is the current temperature; T0 is the initial temperature;

[0052] Step B2, with the thermal expansion and contraction strain of the reservoir matrix, the reservoir matrix expands, and tensile strain is generated around the fracture, when the tensile strain reaches the critical strain ε cr of the material, thermal induced tensile cracking or micro crack initiation is generated in the reservoir matrix or at the edge of the fracture, and the initiation probability P c is characterized by formula (9):

[0053] (9)

[0054] Wherein: n is an empirical constant;

[0055] Step B3, the micro cracks expand and connect under the long-term thermal cycle, forming a secondary heat conduction path, which gradually increases the thermal conductivity k of the reservoir matrix rock mass, but the permeability growth is limited, and the geothermal heat transfer mode is still the conduction type of geothermal heat transfer.

[0056] Further, it also includes:

[0057] In the simulation process of the development process of the fracture in the fault-reservoir combination area, a comprehensive damage criterion D is defined to characterize the fracture deformation and the damage degree of the crack :

[0058] (10)

[0059] Wherein: is the equivalent tensile strain; is the critical tensile strain; is the critical pore pressure; is the critical temperature difference; T is the current temperature; T0 is the initial temperature; p0 is the initial pore pressure; p is the current pore pressure;

[0060] If the comprehensive damage criterion D , representing that the fault-reservoir combination area has locally entered a plastic or damaged state; wherein: if it is a convective geothermal heat transfer type, the area with D≥1 is concentrated on the fault surface and the tip of the fracture; if it is a conductive geothermal heat transfer type, the area with D≥1 is distributed along the high temperature gradient area of ​​the reservoir matrix.

[0061] The numerical simulation method for the development process of fractures in a fault-reservoir combination zone provided by this invention has the following advantages:

[0062] This invention constructs a fault-reservoir coupling model, comprehensively considering the impact of the fault-reservoir coupling relationship on heat transfer and seepage. This allows for accurate simulation of geothermal heat transfer processes in fault-reservoir composite zones, identification of geothermal heat transfer modes, and simulation of fracture propagation behavior and crack evolution under different geothermal heat transfer modes using relevant algorithm models. This invention can accurately simulate the development process of fractures and cracks, thereby revealing the impact of complex structures such as fault-reservoir composite zones on geothermal heat transfer and geothermal fluid migration. It provides a scientific basis for geothermal exploitation and structural safety assessment, possessing high practical value and promising application prospects. Attached Figure Description

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

[0064] Figure 1 A flowchart of the numerical simulation method for the development process of fractures and cracks in the fault-reservoir combination zone provided by the present invention;

[0065] Figure 2 Simulation results of the numerical simulation method for the development process of fracture cracks in the fault-reservoir combination zone provided by the present invention;

[0066] Figure 3 The diagram showing the results of thermal convection during fracture when Ra* is greater than 40 is provided by this invention;

[0067] Figure 4 The diagram shows the results of thermal convection during fracture when Ra* is less than 40, as provided by this invention. Detailed Implementation

[0068] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0069] See Figure 1 This invention provides a numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone, comprising the following steps:

[0070] Step S1: Construct a three-dimensional stratigraphic model of the study area; the three-dimensional stratigraphic model has several fault-reservoir combination zones; the fault-reservoir combination zone is a region with a fault in the reservoir and several fractures distributed in the vicinity of the fault.

[0071] This step is specifically as follows:

[0072] Collect geological data for the study area, including: geological exploration reports, borehole data, seismic exploration data, and hydrogeological test data for the study area;

[0073] The geological data of the study area are analyzed and processed to accurately determine the stratigraphic structural parameters of the study area. The stratigraphic structural parameters include, but are not limited to, the lithology, thickness, distribution range, fault location, fault strike, fault dip angle, and fracture distribution data of the strata; and may further include the porosity, permeability, water conductivity, filling material composition and distribution characteristics of the fault fracture zone.

[0074] Based on the stratigraphic structure parameters of the study area, a three-dimensional model of the stratigraphy was constructed using modeling tools. Specifically, the finite element method was used to construct the three-dimensional model, and an augmented form function combining the Heaviside discontinuity function and the fracture tip term was used to describe the temperature and pressure distribution at the fracture surface at the fracture tip. At the same time, local mesh refinement was performed in the fracture-reservoir combination area to accurately capture the thermal flux gradient at the fracture-reservoir matrix interface.

[0075] Furthermore, fracture subdomains are set up in the fractured reservoir combination area and discretized with high precision. High-order finite element (p-FEM) or extended finite element (XFEM) methods are used for refinement to accurately simulate the temperature field, pressure field and seepage velocity field distribution at the corresponding locations.

[0076] Step S2: On the multiphysics numerical simulation platform, construct the thermal-fluid coupling model of the three-dimensional formation model and the fault-reservoir coupling model of the fault-reservoir combination zone; the fault-reservoir coupling model is a fracture-fracture-reservoir coupling model.

[0077] In this invention, the heat-fluid coupling model is a model used to simulate multiphysics data in the three-dimensional formation model; wherein, the multiphysics data includes temperature field, pressure field and seepage velocity field.

[0078] The present invention does not limit the specific form of the heat-fluid coupling model used. As a preferred approach, the heat-fluid coupling model includes the heat conduction equation and Darcy's law or their equivalent forms.

[0079] Specifically, the interfaces for "Heat Transfer in Porous Media" and "Darcy's Law" can be selected in a multiphysics numerical simulation platform to establish a set of coupled control equations for the heat-fluid-force three fields:

[0080] (11)

[0081] Where: ρ r c is the density of the reservoir rock. r Let be the specific heat capacity of the reservoir rock, t be time, k be the thermal conductivity, and T be the temperature. Let ρ be the temperature change, Q be the intensity of the internal heat source, and ρ be the internal heat source intensity. w For the density of the geothermal fluid, c w The specific heat capacity of geothermal fluids. For geothermal fluid seepage velocity; Darcy's Law: In the formula, k p ρ is the permeability, μ is the dynamic viscosity of the geothermal fluid, and p is the pore pressure.

[0082] Step S3: Set the initial parameters and boundary conditions for the numerical simulation; the initial parameters for the numerical simulation include the initial multiphysics data of the thermal-fluid coupling model, and the initial state parameters of the fractured reservoir of each fractured reservoir coupling model; the boundary conditions include, but are not limited to, thermal boundary and seepage boundary conditions that conform to geological conditions, and can specifically set boundary conditions such as surface, bottom boundary, and fracture surface, including heat flux boundary, isothermal boundary and seepage pressure boundary.

[0083] Step S4: Establish a data exchange channel between the thermal-fluid coupling model and each of the fault-reservoir coupling models. Through the data exchange channel, realize the data exchange and feedback adjustment between the multiphysics field data solved by the thermal-fluid coupling model and the fractured reservoir state parameters solved by the fault-reservoir coupling model. Then, perform multiphysics field simulation on the three-dimensional formation model through the thermal-fluid coupling model, and simulate the fracture development process of the fault-reservoir combination area through the fault-reservoir coupling model.

[0084] This step is specifically as follows:

[0085] Step S41: Based on the multiphysics data of the three-dimensional formation model in the previous time step and the fracture, crack and reservoir state parameters of each fault-reservoir combination zone, solve the thermal-fluid coupling model to obtain the multiphysics data of the three-dimensional formation model in the current time step.

[0086] Step S42: The fault-reservoir coupling model determines the multiphysics data of the fault-reservoir combination area at the current time step based on the multiphysics data of the formation three-dimensional model at the current time step, and determines the temperature change and fracture pore pressure change of the fractures in the fault-reservoir combination area based on the multiphysics data of the fault-reservoir combination area at the current time step; based on the temperature change and fracture pore pressure change of the fractures, the fracture permeability and fracture aperture at the current time step are obtained, and then the fracture state parameters are updated.

[0087] Specifically, the fracture permeability at the current time step is obtained using formulas (1) and (2). and crack opening :

[0088] (1)

[0089] (2)

[0090] in: and , which are the initial values ​​of fracture permeability and fracture aperture at the current time step, respectively; , , and These are the first temperature sensitivity coefficient, the second temperature sensitivity coefficient, the first pressure sensitivity coefficient, and the second pressure sensitivity coefficient, respectively. and , which are the current temperature and the current pore pressure, respectively; and , representing the initial temperature and initial pore pressure, respectively. As the temperature or pressure increases, the fracture aperture and fracture permeability increase accordingly, achieving dynamic feedback of fracture conductivity under heat-fluid coupling conditions.

[0091] Step S43: Based on the updated fracture state parameters, return to step S41. By solving the thermal-fluid coupling model, obtain the multiphysics data of the three-dimensional formation model in the next time step. Repeat steps S41 and S42 in this way to introduce a fracture evolution feedback mechanism in the process of multiphysics data simulation of the three-dimensional formation model.

[0092] Therefore, in solving the heat-fluid coupling model, this invention introduces a heat-fluid-mechanical feedback mechanism for the fracture, updating the fracture permeability k in real time based on changes in temperature and pore pressure in the fracture region. f With crack aperture b f .

[0093] Step S44: In the process of multiphysics data simulation based on fracture evolution feedback mechanism, at a selected time step, the geothermal heat transfer mode of the fracture-reservoir combination zone at the current time step is identified according to the fracture, fracture and reservoir state parameters.

[0094] Based on the identified geothermal heat transfer mode of the fault-reservoir combination area, the development process of the fracture is simulated to obtain the fracture, fracture and reservoir state parameters. The simulated fracture, fracture and reservoir state parameters are then used to update the fracture, fracture and reservoir state parameters in step S41 in real time.

[0095] In this step, the geothermal heat transfer mode of the fault-reservoir combination zone at the current time step is identified based on the parameters of fractures, cracks, and reservoir state. Specifically:

[0096] Taking into account the fracture state parameters, fracture state parameters, and fracture-reservoir matrix coupling relationship, the criterion shown in formula (3) is established. formula:

[0097] (3)

[0098] ,0< ≤1

[0099] Where: g is the gravitational acceleration; β is the thermal expansion coefficient of the geothermal fluid in the fault-reservoir combination area; ΔT is the temperature difference between the upper and lower ends of the fault; H is the fault extension height; k f k represents the crack permeability. m θ is the geothermal reservoir permeability; θ is the fracture dip angle; μ is the kinematic viscosity of the geothermal fluid; λ is the fracture-reservoir matrix coupling coefficient; α T α is the thermal diffusivity; x is the empirical exponent; f The thermal conductivity enhancement coefficient of the crack, η is the thermal conductivity gain factor, L f The crack spacing, Let L be the crack aperture; where, for the case where multiple cracks are distributed in the vicinity of the fracture, the crack spacing L is... f The sum of the distances between all cracks; crack aperture. This is the sum of the openings of all the cracks;

[0100] Criteria of the present invention The formula incorporates the fracture permeability kf and the tilt angle effect sinθ in the numerator to reflect the fracture conductivity and geometric control; the denominator uses terms... Explicitly introducing the reservoir matrix-fracture coupling effect makes it possible to achieve higher permeability (k) when the fracture is much more permeable. f ≫ k m When λ→0, the criterion result degenerates into a convection mode dominated by fracture; when the fracture permeability is close to that of the reservoir matrix (k f ≈k m ), λ≈(1 / 2)α, the criterion reflects that both the fracture and the reservoir matrix contribute to heat transfer; when the fracture has low permeability (k f ≪k m As λ→1, the contribution of convection weakens, and conduction becomes dominant.

[0101] Adopting criteria The formula yields the criterion for determining the storage failure combination region at the current time step. If the criterion If the value is less than or equal to the empirical value, then the geothermal heat transfer mode of the reservoir failure combination area at the current time step is identified as: conductive geothermal heat transfer type; if the criterion is... The empirical value identifies the geothermal heat transfer mode of the reservoir failure combination area at the current time step as: convection-type geothermal heat transfer. The empirical value is generally set to 40.

[0102] Traditional Rayleigh number criteria do not consider the differences between fractures and reservoir matrix, as well as the geometric effects of fractures, leading to significant biases in the discrimination results in fracture-developed areas. This invention proposes a new criterion as shown in formula (3). .

[0103] In this step, the development process of fractures is simulated based on the identified geothermal heat transfer pattern of the fault-reservoir combination area, specifically distinguishing the following two cases:

[0104] (a) The geothermal heat transfer mode is convection-type geothermal heat transfer:

[0105] Step A1: If the geothermal heat transfer mode is convection-type, the convection of the hot fluid in the fracture channel is enhanced, and both the local temperature gradient and fluid pressure on the fracture surface increase, generating thermal stress σ in the direction of the fracture dip angle θ. T With pore pressure stress σ P The coupling superposition effect is used to obtain the effective normal stress in the direction of the fracture dip angle θ using formula (4). :

[0106] (4)

[0107] in: The total normal stress is θ in the direction of the fracture dip angle. p is the pore pressure coupling coefficient; p is the current pore pressure; p0 is the initial pore pressure; is the coefficient of thermal expansion; E is the elastic modulus, obtained from well logging data or stress-strain curves; T is the current temperature; T0 is the initial temperature;

[0108] Step A2: Using formula (5), the shear stress in the direction of the fracture dip angle θ is obtained. :

[0109] (5)

[0110] in: These are the in-plane shear stress components; This represents the shear stress component in the vertical direction;

[0111] Step A3, determine whether the sliding condition of formula (6) is true:

[0112] (6)

[0113] in: This refers to the shear strength that the fracture surface still possesses under zero normal stress. It is the internal friction angle;

[0114] Step A4: If the slip condition of formula (6) is met, the fracture enters the shear slip stage. If the slip continues or the cumulative displacement δ s Exceeding the threshold δ cr If the fracture surface expands, the cracks adjacent to the fracture extend along the fault direction to form secondary cracks. The expansion rate of the secondary cracks satisfies equation (7). The crack expansion direction is consistent with the direction of the maximum principal stress, forming a heat conduction channel, which leads to an increase in the overall permeability and thermal conductivity of the fracture zone, and triggers continuous deformation and thermally driven slip of the fracture morphology.

[0115] (7)

[0116] Where: L f t represents the crack spacing; k represents time; k represents the crack spacing. v is the crack growth coefficient; m is the empirical index.

[0117] (ii) The geothermal heat transfer mode is the conductive type:

[0118] Step B1: If the geothermal heat transfer mode is conductive geothermal heat transfer, in the fault-reservoir combination area, geothermal energy diffuses through the reservoir matrix via conduction, resulting in a high local temperature gradient and causing thermal expansion and contraction strain within the reservoir matrix. The strain value of the reservoir matrix is ​​obtained using formula (8). :

[0119] (8)

[0120] in: T is the first temperature sensitivity coefficient; T is the current temperature; T0 is the initial temperature;

[0121] Step B2: As the reservoir matrix expands and contracts due to thermal expansion and contraction, tensile strain is generated around the fractures. When the tensile strain reaches the material's critical strain ε... cr At that time, it is usually 10 -4 -10 -3 Thermally induced tensile fracturing or microfracture initiation occurs within the reservoir matrix or at fracture edges, with an initiation probability P. c Characterized by formula (9):

[0122] (9)

[0123] Where: n is an empirical constant;

[0124] In step B3, the microfractures expand and connect under long-term thermal cycling, forming secondary heat conduction paths. This leads to a gradual increase in the thermal conductivity k of the reservoir matrix rock, but the permeability increase is limited, and the geothermal heat transfer mode remains the conduction-type geothermal heat transfer type.

[0125] In this invention, during the simulation of the fracture and crack development process in the fractured reservoir combination zone, a comprehensive damage criterion is defined to characterize the degree of fracture deformation and crack damage. :

[0126] (10)

[0127] in: For equivalent tensile strain; The critical tensile strain; The critical pore pressure; Where T is the critical temperature difference; T is the current temperature; T0 is the initial temperature; p0 is the initial pore pressure; p is the current pore pressure;

[0128] If the comprehensive damage criteria , representing that the fault-reservoir combination area has locally entered a plastic or damaged state; wherein: if it is a convective geothermal heat transfer type, the area with D≥1 is concentrated on the fault surface and the tip of the fracture; if it is a conductive geothermal heat transfer type, the area with D≥1 is distributed along the high temperature gradient area of ​​the reservoir matrix.

[0129] Therefore, this invention introduces a thermally induced crack evolution mechanism in the process of solving the coupled heat-fluid model, based on real-time data from temperature T and crack pressure. Update crack permeability k f With crack aperture b f According to the updated fracture permeability k f With crack aperture b f Calculate the improved criterion Ra * The geothermal heat transfer mode is determined accordingly; the geothermal heat transfer mode includes conductive geothermal heat transfer type and convective geothermal heat transfer type; further, based on the simulated temperature field and seepage field data, combined with the geothermal heat transfer mode, the convective fracture slip propagation behavior and conductive reservoir matrix microfracture initiation are analyzed; finally, the thermo-mechanical damage criterion D is defined to quantitatively assess the stability of the fracture zone.

[0130] The following is an example:

[0131] Addressing the shortcomings of existing technologies using Ra and Pe, the indices for determining convection and conduction in homogeneous media, this embodiment provides a numerical simulation method for the development of shallow and deep fractures in tectonically complex regions such as fault-reservoir composite areas. (Refer to...) Figures 2 to 4By comprehensively considering the fracture characteristics and the relationship between the fault and the reservoir, the geothermal heat transfer mode is identified. Then, the fracture deformation and microcrack evolution process under different geothermal heat transfer modes are accurately simulated, and the feasibility is verified by simulation.

[0132] Step 1: Collect geological data for the study area: borehole data, well logging data, core physical parameters, formation strike and dip, and fault and fracture distribution maps. Construct a 3D formation model to determine the geometric framework of the fracture-reservoir matrix, fracture dip angle θ, fracture extension height H, and the contact relationship between the fracture and the reservoir.

[0133] First, geological data for the study area is collected, including geological exploration reports, borehole data, seismic exploration data, and hydrogeological test data. Through analysis and processing of this data, parameters such as lithology, thickness, distribution range of strata, and the location, strike, dip angle, fracture zone characteristics, and permeability of complex structures like faults are determined. The influencing factors of faults mainly include the porosity, permeability, hydraulic conductivity, composition and distribution characteristics of the filling material in the fault fracture zone. Then, a corresponding three-dimensional stratigraphic model is constructed according to the determined parameters. Attention must be paid to the contact relationships between strata and the intersection relationships between faults and reservoirs to ensure the model's realism and accuracy.

[0134] Step 2: Collect or measure physical property parameters, including the basic parameters required for thermal-fluid coupling, including: 1. Reservoir matrix parameters: reservoir rock density ϼ r Specific heat capacity of reservoir rocks c r Thermal conductivity k, elastic modulus E, first temperature sensitivity coefficient 2. Crack parameters: Initial value of crack permeability k f0 Initial crack aperture value b0, crack fracture energy G f 3. Geothermal fluid parameters: Geothermal fluid density ϼ w Geothermal fluid specific heat capacity c w , where μ is the dynamic viscosity of the geothermal fluid and β is the thermal expansion coefficient of the geothermal fluid.

[0135] The thermal conduction module and the porous media flow module in COMSOL were selected to implement a heat-fluid coupling simulation. Considering that the presence of fractures alters the thermal conductivity and permeability of the formation, these fractures were treated as regions with specific thermal conductivity and permeability. Based on the characteristics of the fracture zone, lithology, and hydrogeological test data, the thermal conductivity and permeability were obtained through experiments, field tests, or literature reviews. Simultaneously, the initial temperature field of the formation was determined based on the geothermal gradient and surface temperature of the study area; the initial seepage field was determined based on the regional groundwater level distribution and hydraulic gradient. The geothermal gradient could be obtained through borehole temperature measurement data, and the surface temperature could be the multi-year average air temperature of the study area.

[0136] Step 3: Construct a 3D formation model containing faults, fractures, and reservoirs, including fault planes, fracture zones (which can be represented by surfaces or narrow domains), reservoir interfaces, and external surrounding rock boundaries.

[0137] Mesh Generation: The constructed 3D formation model is meshed, using either free tetrahedral or hexahedral meshes globally. Mesh refinement is applied in complex structural areas such as fault-reservoir assemblages, areas with large heat gradients such as stratigraphic boundaries, and areas with drastic seepage changes near fault fracture zones. For example, the minimum element size is ≤ [missing information - likely a value]. or ≤ d min (Target resolution) This is the initial value for crack aperture to improve simulation accuracy;

[0138] In regions with small heat transfer gradients and seepage variations, the mesh density can be appropriately reduced to improve computational efficiency; at crack tips / interfaces, local size reduction (at least 3 layers of fine mesh rings) should be used; mesh quality should be evaluated using parameters such as mesh distortion and aspect ratio to ensure that the mesh quality meets computational requirements.

[0139] Step 4: Discretization method determination:

[0140] If the crack path is known and relatively simple, p-FEM (higher-order FEM) is preferred; if the crack path is uncertain or there is interaction between multiple cracks, XFEM is used. The order of p-FEM is p = 2–5 (generally 2 or 3 is used for experimentation); the mesh refinement factor at the crack interface is 3–10.

[0141] Step 5: Determine the thermal-fluid coupling equations and set the solver parameters, including time step, number of iterations, convergence accuracy, etc.

[0142] The settings are determined according to the simulation timescale (e.g., short-term, medium-term, long-term) and accuracy requirements. The solution process combines the heat conduction equation and Darcy's law. In the multiphysics platform, the interfaces of "heat transfer in porous media" and "Darcy's law" are selected to establish the three-field coupled control equation set of heat-fluid-force as shown in formula (11).

[0143] For transient simulations, the time step should be chosen considering the rate of temperature and seepage change to ensure the stability of the simulation results. For steady-state simulations, a larger number of iterations and higher convergence accuracy can be set to obtain stable temperature and seepage field distributions. The COMSOL software is then used for coupled solution to obtain the distribution and changes of the temperature, seepage, and pressure fields in the formation system, clarifying the vertical migration and lateral flow characteristics of geothermal fluids.

[0144] Step 6: Introduce a thermal-fluid-mechanical feedback mechanism in the coupled solution to establish the fracture permeability k in formulas (1) and (2). fThe dynamic relationship between temperature T and pore pressure p is reflected by updating fracture parameters in real time, showing the influence of temperature rise and fluid pressure change on fracture conductivity. When the temperature rises or the fluid pressure increases, the fracture aperture and fracture permeability increase accordingly, realizing dynamic feedback of fracture conductivity under thermal-fluid coupling conditions.

[0145] Step 7: Propose an improved Rayleigh criterion formula:

[0146] The traditional Rayleigh number criterion does not consider the differences between fracture and reservoir matrix and the geometric effects of fracture, resulting in a large deviation in the discrimination results in fracture-developed areas. Therefore, a new criterion as shown in formula (3) is proposed. formula.

[0147] Based on the above comparison of fracture-reservoir matrix permeability and the combined results of numerical simulation and field verification, the following calculations were made: Value, Improved Criterion The critical value is referenced from the convection theory of porous media, and approximately 40 is taken as the criterion for the occurrence of convection. When When the value is significantly less than this, conduction dominates, indicating a conductive geothermal heat transfer type; when... When the value is significantly greater than this, convection dominates, indicating a convective geothermal transfer type. The specific applicable range of the critical value can be corrected based on numerical simulation and field observation data.

[0148] Step 8: Based on the temperature field and geothermal fluid transport results obtained from the simulation, the development process of fractures and cracks is simulated for both convective and conductive geothermal heat transfer types.

[0149] Step 9: To uniformly characterize the degree of fracture deformation and microcrack damage, a comprehensive damage criterion as shown in formula (10) is defined. .

[0150] Therefore, the improved criterion not only reveals the heat transfer mechanism under different mechanistic conditions, but also demonstrates the influence of thermal field changes on the mechanical response of the fracture system, providing a basis for geothermal extraction and structural safety assessment.

[0151] The numerical simulation method for the development process of fractures in a fault-reservoir combination zone provided by this invention can be summarized as follows:

[0152] Step 1: Based on the geological, borehole and geophysical data of the study area, construct a three-dimensional stratigraphic model containing fractures, cracks and reservoirs; select the "Heat Transfer in Porous Media" and "Darcy's Law" interfaces in the multiphysics platform to establish a heat-fluid coupling model, set heat and seepage boundary conditions consistent with the geological environment, and establish a fracture-fracture-reservoir multiphysics coupling model.

[0153] Step 2: Set up a fracture subdomain in the fracture development area of ​​the fault-reservoir combination zone, and use the high-order finite element method (p-FEM) or extended finite element method (XFEM) to characterize the temperature and velocity gradient changes inside the fracture; set the solver parameters and perform steady-state or transient solutions to obtain the temperature field, pressure field and seepage velocity field distribution.

[0154] Step 3: Introduce a thermal-fluid-mechanical feedback mechanism in the coupled solution to establish the dynamic relationship between fracture permeability kf and temperature T and pore pressure p. By updating fracture parameters in real time, the influence of temperature rise and fluid pressure change on fracture conductivity can be reflected.

[0155] Step 4: Propose an improved Rayleigh criterion formula and calculate the improved criterion Ra. ∗ According to Ra * The value can be used to distinguish between conductive and convective geothermal heat transfer modes.

[0156] Step 5: Based on the temperature field and fluid transport results obtained from the simulation, analyze the response mode of the formation structure. Under convection conditions, analyze whether the fracture may slip or propagate. Under conduction conditions, analyze whether it will induce fine cracks at the microscale.

[0157] Step 6: Define the thermo-mechanical coupling damage criterion D to comprehensively describe the risk of thermally induced deformation in fracture or reservoir matrix regions.

[0158] This invention proposes an improved criterion formula by innovatively introducing a comparison between fracture and reservoir matrix permeability, as well as factors related to fracture geometry, effectively overcoming the limitations of traditional Ra... The current criterion is only applicable to homogeneous porous media and lacks applicability in complex structural areas. This invention can more accurately identify the types of geothermal modes dominated by conduction and convection. By constructing a heat-fluid coupling model based on a multiphysics numerical simulation platform, selecting heat transfer and seepage analysis interfaces, and combining high-precision discretization methods to process fracture development zones, a heat-fluid-mechanical feedback mechanism is introduced to update fracture parameters in real time. This achieves a comprehensive reflection of fracture fluid migration and reservoir matrix heat conduction, making the model results more consistent with the characteristics of geothermal systems. All parameters can be obtained from drilling, logging, and experimental data, possessing verifiability and generalizability, ensuring the practicality and reliability of the method. Under a unified simulation framework, the fracture propagation behavior is analyzed separately for convection and conduction geothermal systems, further improving the study of the dynamic changes of fractures in geothermal systems in complex structural areas. This method not only reveals the formation mechanism of geothermal systems in fracture-controlled areas but also provides reliable technical support for the exploration, development, and structural safety assessment of geothermal resources, while having strong applicability and promotional value.

[0159] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone, characterized in that, Includes the following steps: Step S1: Construct a three-dimensional stratigraphic model of the study area; the three-dimensional stratigraphic model has several fault-reservoir combination zones; the fault-reservoir combination zone is a region with a fault in the reservoir and several fractures distributed in the vicinity of the fault. Step S2: Construct the thermal-fluid coupling model of the three-dimensional formation model and the fault-reservoir coupling model of the fault-reservoir combination zone; the fault-reservoir coupling model is a fracture-fracture-reservoir coupling model. Step S3: Set the initial parameters and boundary conditions for the numerical simulation; the initial parameters for the numerical simulation include the initial multiphysics data of the thermal-fluid coupling model, and the initial state parameters of the fractured reservoir for each fractured reservoir coupling model. Step S4: Establish a data exchange channel between the thermal-fluid coupling model and each of the fault-reservoir coupling models. Through the data exchange channel, realize the data exchange and feedback adjustment between the multiphysics field data solved by the thermal-fluid coupling model and the fractured reservoir state parameters solved by the fault-reservoir coupling model. Then, perform multiphysics field simulation on the three-dimensional formation model through the thermal-fluid coupling model, and simulate the fracture development process of the fault-reservoir combination area through the fault-reservoir coupling model. Step S4 is as follows: Step S41: Based on the multiphysics data of the three-dimensional formation model in the previous time step and the fracture, crack and reservoir state parameters of each fault-reservoir combination zone, solve the thermal-fluid coupling model to obtain the multiphysics data of the three-dimensional formation model in the current time step. Step S42: The fault-reservoir coupling model determines the multiphysics data of the fault-reservoir combination area at the current time step based on the multiphysics data of the formation three-dimensional model at the current time step, and determines the temperature change and fracture pore pressure change of the fractures in the fault-reservoir combination area based on the multiphysics data of the fault-reservoir combination area at the current time step; based on the temperature change and fracture pore pressure change of the fractures, the fracture permeability and fracture aperture at the current time step are obtained, and then the fracture state parameters are updated. Step S43: Based on the updated fracture state parameters, return to step S41. By solving the thermal-fluid coupling model, obtain the multiphysics data of the three-dimensional formation model in the next time step. Repeat steps S41 and S42 in this way to introduce a fracture evolution feedback mechanism in the process of multiphysics data simulation of the three-dimensional formation model. Step S44: In the process of multiphysics data simulation based on fracture evolution feedback mechanism, at a selected time step, the geothermal heat transfer mode of the fracture-reservoir combination zone at the current time step is identified according to the fracture, fracture and reservoir state parameters. Based on the identified geothermal heat transfer mode of the fault-reservoir combination area, the development process of the fracture is simulated to obtain the fracture, fracture and reservoir state parameters. The simulated fracture, fracture and reservoir state parameters are then used to update the fracture, fracture and reservoir state parameters in step S41 in real time.

2. The numerical simulation method for the development process of fractures in a fault-reservoir combination zone according to claim 1, characterized in that, Step S1 is as follows: Collect geological data for the study area, including: geological exploration reports, borehole data, seismic exploration data, and hydrogeological test data for the study area; Geological data of the study area are analyzed and processed to determine the stratigraphic structure parameters of the study area. Based on the stratigraphic structure parameters of the study area, a three-dimensional model of the stratigraphy was constructed using modeling tools.

3. The numerical simulation method for the development process of fractures in a fault-reservoir combination zone according to claim 2, characterized in that, The stratigraphic structural parameters include the lithology, thickness, distribution range, fracture location, fracture strike, fracture dip angle, and fracture distribution data of the strata; The specific steps of constructing the three-dimensional formation model using modeling tools are as follows: the three-dimensional formation model is constructed using the finite element method, and an augmented form function combining the Heaviside discontinuity function and the fracture tip term is used to describe the distribution of temperature and pressure on the fracture surface at the fracture tip. At the same time, local mesh refinement is performed in the fracture-reservoir combination zone.

4. The numerical simulation method for the development process of fractures in a fault-reservoir combination zone according to claim 1, characterized in that, The thermal-fluid coupling model is used to simulate the multiphysics data in the three-dimensional formation model; wherein, the multiphysics data includes temperature field, pressure field and seepage velocity field.

5. The numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone according to claim 1, characterized in that, In step S42, the fracture permeability at the current time step is obtained using formulas (1) and (2). and crack opening : （1） （2） in: and These are the initial values ​​of fracture permeability and fracture aperture at the current time step, respectively. , , and These are the first temperature sensitivity coefficient, the second temperature sensitivity coefficient, the first pressure sensitivity coefficient, and the second pressure sensitivity coefficient, respectively. and These are the current temperature and the current pore pressure, respectively. and These are the initial temperature and the initial pore pressure, respectively.

6. The numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone according to claim 1, characterized in that, The process of identifying the geothermal heat transfer mode of the fault-reservoir combination zone at the current time step based on fault, fracture, and reservoir state parameters specifically includes: Taking into account the fracture state parameters, fracture state parameters, and fracture-reservoir matrix coupling relationship, the criterion shown in formula (3) is established. formula: （3） ，0< ≤1 Where: g is the gravitational acceleration; β is the thermal expansion coefficient of the geothermal fluid in the fault-reservoir combination area; ΔT is the temperature difference between the upper and lower ends of the fault; H is the fault extension height; k f k represents the crack permeability. m θ is the geothermal reservoir permeability; θ is the fracture dip angle; μ is the kinematic viscosity of the geothermal fluid; λ is the fracture-reservoir matrix coupling coefficient; α T α is the coefficient of thermal expansion; α is the empirical exponent; x f The thermal conductivity enhancement factor of the crack, η is the thermal conductivity gain factor, L f The crack spacing, Let L be the crack aperture; where, for the case where multiple cracks are distributed in the vicinity of the fracture, the crack spacing L is... f The sum of the spacing between all cracks; crack aperture. This is the sum of the openings of all the cracks; Adopting criteria The formula yields the criterion for determining the storage failure combination region at the current time step. If the criterion If the value is less than or equal to the empirical value, then the geothermal heat transfer mode of the reservoir failure combination area at the current time step is identified as: conductive geothermal heat transfer type; if the criterion is... Based on the empirical value, the geothermal heat transfer mode of the reservoir failure combination area at the current time step is identified as: convection geothermal heat transfer type.

7. The numerical simulation method for the development process of fractures in a fault-reservoir combination zone according to claim 1, characterized in that, The process of simulating the development of fractures based on the identified geothermal heat transfer patterns in the fault-reservoir combination area is as follows: Step A1: If the geothermal heat transfer mode is convection-type, the convection of the hot fluid in the fracture channel is enhanced, and both the local temperature gradient and fluid pressure on the fracture surface increase, generating thermal stress σ in the direction of the fracture dip angle θ. T With pore pressure stress σ P The coupling superposition effect is used to obtain the effective normal stress in the direction of the fracture dip angle θ using formula (4). : （4） in: The total normal stress is θ in the direction of the fracture dip angle. p is the pore pressure coupling coefficient; p is the current pore pressure; p0 is the initial pore pressure; Where is the coefficient of thermal expansion; E is the elastic modulus; T is the current temperature; T0 is the initial temperature; Step A2: Using formula (5), the shear stress in the direction of the fracture dip angle θ is obtained. : （5） in: These are the in-plane shear stress components; This represents the shear stress component in the vertical direction; Step A3, determine whether the sliding condition of formula (6) is true: （6） in: This refers to the shear strength that the fracture surface still possesses under zero normal stress. It is the internal friction angle; Step A4: If the slip condition of formula (6) is met, the fracture enters the shear slip stage. If the slip continues or the cumulative displacement δ s Exceeding the threshold δ cr If the fracture surface expands, the cracks adjacent to the fracture extend along the fault direction to form secondary cracks. The expansion rate of the secondary cracks satisfies equation (7). The crack expansion direction is consistent with the direction of the maximum principal stress, forming a heat conduction channel, which leads to an increase in the overall permeability and thermal conductivity of the fracture zone, and triggers continuous deformation and thermally driven slip of the fracture morphology. （7） Where: L f t represents the crack spacing; k represents time; k represents the crack spacing. v is the crack growth coefficient; m is the empirical index.

8. The numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone according to claim 1, characterized in that, The process of simulating the development of fractures based on the identified geothermal heat transfer patterns in the fault-reservoir combination area is as follows: Step B1: If the geothermal heat transfer mode is conductive geothermal heat transfer, in the fault-reservoir combination area, geothermal energy diffuses through the reservoir matrix via conduction, resulting in a high local temperature gradient and causing thermal expansion and contraction strain within the reservoir matrix. The strain value of the reservoir matrix is ​​obtained using formula (8). : （8） in: T is the first temperature sensitivity coefficient; T is the current temperature; T0 is the initial temperature; Step B2: As the reservoir matrix expands and contracts due to thermal expansion and contraction, tensile strain is generated around the fractures. When the tensile strain reaches the material's critical strain ε... cr At this time, thermally induced tensile fracturing or microfracture initiation occurs within the reservoir matrix or at the fracture edge, with an initiation probability P. c Characterized by formula (9): （9） Where: n is an empirical constant; In step B3, the microfractures expand and connect under long-term thermal cycling, forming secondary heat conduction paths. This leads to a gradual increase in the thermal conductivity k of the reservoir matrix rock, but the permeability increase is limited, and the geothermal heat transfer mode remains the conduction-type geothermal heat transfer type.

9. The numerical simulation method for the development process of fractures and cracks in a fault-reservoir combination zone according to claim 1, characterized in that, Also includes: In simulating the fracture and crack development process in the aforementioned fault-reservoir combination zone, a comprehensive damage criterion is defined to characterize the degree of fracture deformation and crack damage. : （10） in: For equivalent tensile strain; The critical tensile strain; The critical pore pressure; Where T is the critical temperature difference; T is the current temperature; T0 is the initial temperature; p0 is the initial pore pressure; p is the current pore pressure; If the comprehensive damage criteria , representing that the fault-reservoir combination area has locally entered a plastic or damaged state; wherein: if it is a convective geothermal heat transfer type, the area with D≥1 is concentrated on the fault surface and the tip of the fracture; if it is a conductive geothermal heat transfer type, the area with D≥1 is distributed along the high temperature gradient area of ​​the reservoir matrix.

Images