Heat flow curing multi-field coupling modeling method based on enhanced geothermal system

Abstract

The invention discloses a heat flow curing multi-field coupling modeling method based on an enhanced geothermal system, which belongs to the technical field of oil and gas field exploitation, and comprises the following steps: generating two groups of random numbers as crack center positions, and determining random crack end point coordinates according to the crack center positions, trace lengths and directions; forming an initial smooth fracture grid according to the determined random fracture endpoint coordinates; according to the initial smooth fracture grid, roughness continuity of each fracture path is controlled through a Hurst index, and a random fracture path containing roughness is generated; according to the generated random crack path containing the roughness, a hot dry rock one-injection and two-production random crack model is established; setting material parameters of the random fracture model according to logging data, and adding a stress field, a seepage field, a temperature field and a chemical field to the random fracture model; initial conditions and boundary conditions are set for a stress field, a seepage field, a temperature field and a chemical field of the random fracture model, and multi-physics field coupling setting is carried out.

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas field exploitation, and in particular to a multi-field coupling modeling method based on heat flow solidification of an enhanced geothermal system. Background Art

[0002] Numerical simulation of multi-physics coupling in enhanced geothermal systems (EGS) is an important research topic. Existing EGS models are mostly limited to heat-fluid (TH), fluid-solid (HM), or heat-fluid-solid coupled (THM). They do not consider the dynamic changes in porosity and permeability caused by mineral dissolution, and ignore the time-varying effects of the chemical field (C), resulting in distortion of the dynamic process: first, the dissolution of minerals under temperature gradients will increase permeability in the short term. Second, mineral reactions are often accompanied by endothermic / exothermic processes. Ignoring the heat of chemical reactions will lead to reconstruction of the reservoir temperature field, resulting in an overestimation of the thermal energy extraction efficiency. Third, the rock mechanical parameters continue to weaken due to mineral dissolution, exacerbating the risk of shear fracture.

[0003] The above content is only used to assist in understanding the technical solution of the present invention and does not constitute an admission that the above content is prior art. Summary of the Invention

[0004] The main purpose of the present invention is to provide a multi-field coupling modeling method based on heat flow solidification of an enhanced geothermal system, aiming to solve or partially solve the above problems.

[0005] To achieve the above objectives, the present invention provides a multi-field coupled modeling method based on heat flow solidification of an enhanced geothermal system, the multi-field coupled modeling method based on heat flow solidification of an enhanced geothermal system comprising: Generate two sets of random numbers as the crack center position, and determine the random crack endpoint coordinates according to the crack center position, trace length and direction, and form an initial smooth crack grid based on the determined random crack endpoint coordinates; Based on the initial smooth crack grid, the roughness persistence of each crack path is controlled by the Hurst exponent to generate a random crack path containing roughness; Based on the generated random fracture paths with roughness, a random fracture model of hot dry rock one-injection-two-production is established; According to the well logging data, the material parameters of the random fracture model are set, and the stress field, seepage field, temperature field, and chemical field are added to the random fracture model; Set initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model, and perform multi-physics field coupling settings; Meshing the random fracture model according to different regions of the random fracture model, meshing the fracture region using an ultra-fine free triangle meshing method, and meshing the remaining reservoir regions using a refined free quadrilateral meshing method; Conduct analysis on the evolution characteristics of stress field, seepage field, temperature field and chemical field, calculate the recovery degree of reservoir through post-processing of surface average value; select different locations from the injection well, calculate the temperature change, and simulate the change of fracturing fluid viscosity at different locations through the relationship between temperature and fracturing fluid viscosity.

[0006] Preferably, in the multi-field coupling modeling method based on thermal flow solidification of the enhanced geothermal system, in the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physics field coupling setting, the coupling setting of the stress field includes: Set up the solid mechanics module, set the model boundary conditions to fixed constraints, set the initial displacement and initial velocity to 0, use linear elastic materials and add thermal stress coupling terms, set the maximum horizontal principal stress and minimum horizontal principal stress, calculate the displacement of the crack based on the pore fluid pressure provided by the seepage field and the temperature provided by the temperature field, update the crack aperture, output volume strain and transmit it to the seepage field, so that the stress drives the crack deformation and changes the permeability; The calculation formula is as follows: ; (1) Where G is the shear modulus, Pa; u1 is the displacement component, m; λ is the Lame constant, Pa; α B is the Biot coefficient, dimensionless; P is the pore fluid pressure, Pa; E is the elastic modulus, Pa; ν is Poisson's ratio, dimensionless; α T is the coefficient of thermal expansion, 1 / K; T is the local average temperature under the combined action of fluid flow and solid heat conduction, unit K; α C is the chemical expansion coefficient, 1 / (mol / m 3 ); C is the chemical concentration field, mol / m 3 ; F is the body force, N.

[0007] Preferably, in the multi-field coupling modeling method based on thermal flow solidification of the enhanced geothermal system, in the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physics field coupling setting, the coupling setting of the seepage field includes: Set up the seepage field module, set the boundary condition to no flow, set the initial value to formation pressure, set the parameters of the fluid and matrix in the reservoir matrix and porous media and add porous elastic storage, define the fluid in the fracture and the fracture geometry properties, and form a driving pressure drop drive by setting the pressure difference, where the injection well pressure is higher than the formation pressure and the production well pressure is lower than the formation pressure. According to the volume strain ε provided by the stress field V The temperature provided by the temperature field and the porosity ϕ provided by the chemical field are used to calculate the Darcy velocity and transfer it to the temperature field and the dilute species transfer field. The output pore fluid pressure is fed back to the solid mechanics module. The specific calculation formula is as follows: The reservoir flow equation is: ; (2) The seepage equation in the crack is: ; (3) in, ϕ is the porosity, dimensionless; ρ f is the fluid density, kg / m 3 ; c t is the compressibility coefficient, 1 / Pa; P is the pore fluid pressure, Pa; t is time, s; u is the Darcy velocity of the fluid in the reservoir, m / s; α B is the Biot coefficient, dimensionless; ε V is the volume strain, ε V =ε xx +ε yy +ε zz , dimensionless, ε xx is the volume strain in the x direction, ε yy is the volume strain in the y direction, ε zz is the volume strain in the z direction; Q f is the source and sink term of the fluid, kg / (m 3 ·s); ϕ chem is the porosity change caused by chemical reaction, dimensionless; R chem is the reaction rate, mol / (m3 ·s); w is the effective opening of the crack, m; ; d f is the initial crack opening, m; u n is the normal displacement increment, m; w chem is the change in chemically induced crack width, m; r f is the fluid density, kg / m 3 ; u f is the Darcy velocity of the fluid in the fracture, m / s; M is the molar mass of the substance, kg / mol.

[0008] Preferably, in the multi-field coupling modeling method based on thermal flow solidification of the enhanced geothermal system, in the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physics field coupling setting, the coupling setting of the temperature field includes: Set up the temperature field module, set the boundary conditions of the reservoir and fractures to thermal insulation, set the initial condition to the initial temperature of the rock mass, set the parameters of the fluid and matrix in the reservoir matrix and porous medium, define the fluid and fracture geometry properties in the fracture, set the temperature of the injection well, calculate the temperature field based on the flow velocity provided by the seepage field, the heat source provided by the chemical field, and the displacement provided by solid mechanics, and transfer it to the seepage field, stress field, and chemical field respectively; The calculation formula is as follows: The temperature field equation in the reservoir is: ; (4) The temperature field equation in the crack is: ; (5) Where, (ρC) eff is the effective heat capacity, J / (m 3 ⋅K); ; ϕ is the porosity, dimensionless; ρ m is the fluid density, kg / m 3 ; C f and C m are the specific heat capacities of the fluid and reservoir, J / (kg·K), respectively; ρ f is the fluid density, kg / m 3 ; T is the local average temperature under the combined action of fluid flow and solid heat conduction, unit K; u is the Darcy velocity of the fluid in the reservoir, m / s; λ eff is the thermal conductivity of rock, W / (m·K); ; λ rock is the thermal conductivity of rock, W / (m·K); λ f is the thermal conductivity of the fluid, W / (m·K); W f is the fluid heat source, W / m 3 ; Z is the reaction heat source, W / m 3 ; w is the effective opening of the crack, m; ; d f is the initial crack opening, m; u n is the normal displacement increment, m; w chem is the change in chemically induced crack width, m; t is time, s; u f is the Darcy velocity of the fluid in the fracture, m / s.

[0009] Preferably, in the multi-field coupling modeling method based on thermal flow solidification of the enhanced geothermal system, in the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physical field coupling setting, the coupling setting of the chemical field includes: Set up the chemical reaction module, add the following reaction equation, set the initial concentration, calculate the reaction rate according to the following formula (6) and transfer it to the dilute species transfer field. According to formula (7), input the temperature T into the heat source and output the reaction heat source to the temperature field; The reaction equation is: ; Formula (6) and formula (7) are as follows: ; (6) ; (7) Among them, R chem is the reaction rate, mol / (m 3 ·s); k0 is the prefactor, mol / (m 2 ·s); E0 is the activation energy, J / mol; R is the ideal gas constant, 8.314 J / (mol·K); T is the local average temperature under the combined action of fluid flow and solid heat conduction, K; K eq (T) is the temperature-dependent equilibrium constant, dimensionless; Q is the ion activity product, unitless; Z is the reaction heat source, W / m 3 ; R chem，i is the reaction rate of the ith reaction, mol / (m 3 ·s); △H i is the reaction enthalpy change of the i-th reaction, J / mol; i is the reaction category, dimensionless; n is the total number of independent chemical reactions in the system and is dimensionless.

[0010] Preferably, in the multi-field coupling modeling method based on thermal flow solidification of the enhanced geothermal system, in the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physical field coupling setting, the chemical field includes a dilute species transfer field, and the coupling setting of the dilute species transfer field includes: Set up the dilute species transfer module, set the boundary conditions of the reservoir to no flux, set the initial value of the fracturing fluid injected into the reservoir to 0, set the basic parameters of the fluid in the reservoir matrix and porous medium, couple the velocity with the seepage field, set the basic parameters of the fluid and material in the fracture, add a reaction in the fracture, and couple its reaction rate with the reaction rate in the chemical reaction. Set the injection well to concentration flux and the outlet to no flux. According to the flow velocity u provided by the seepage field, the flow rate u is calculated by formula (8). f and the reaction rate R provided by the chemical reaction field chem , calculate the ion concentration and pass it to the chemical reaction module to drive the reaction balance; Formula (8) is as follows: ; (8) in, w is the effective opening of the crack, m; C is the chemical concentration field, mol / m 3 ; t is time, s; u f is the Darcy velocity of the fluid in the fracture, m / s; D f is the molecular diffusion coefficient, m 2 / s; Rchem is the reaction rate, mol / (m 3 ·s).

[0011] Preferably, in the multi-field coupling modeling method based on thermal flow solidification of the enhanced geothermal system, the initial conditions and boundary conditions are set for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively, and the multi-physics field coupling setting is performed, including: In COMSOL, multi-physics field coupling is set up, solid mechanics field and seepage field are coupled into poroelasticity, and temperature field and solid mechanics field are coupled into thermal expansion; the seepage field affects the convective heat transfer process of the temperature field, and the temperature field affects the seepage field by affecting the physical parameters of the fluid. The effective stress in the stress field can change the permeability of the fractures in the reservoir and thus affect the seepage field. The poroelasticity occurring in the seepage field affects its temperature field. The chemical reaction field affects the permeability and porosity of the rock, thereby affecting fluid flow and heat transfer. The temperature field, seepage field and stress field affect the chemical reaction rate.

[0012] Preferably, in the multi-field coupled modeling method based on heat flow solidification of an enhanced geothermal system, in the step of calculating the recovery degree of the reservoir; selecting different locations from the injection well, calculating the temperature change, and simulating the change in fracturing fluid viscosity at different locations through the relationship between temperature and fracturing fluid viscosity, the formula for calculating the recovery degree of the reservoir is: ; (9) Wherein, η is the degree of recovery, dimensionless; V s is the reservoir volume, m 3 ; ρ s is the reservoir rock density, kg / m 3 ; C s is the heat capacity of reservoir rock, J / (kg·K); T in is the initial temperature of the reservoir, °C; T av (t) is the comprehensive average temperature of the reservoir at a certain moment, °C; V f is the reservoir volume, m 3 ; Z chem is the heat source of chemical reaction, W / m 3 ; T inj is the injection fluid temperature, K.

[0013] Preferably, in the multi-field coupled modeling method based on heat flow solidification of an enhanced geothermal system, the step of controlling the roughness persistence of each fracture path by using the Hurst exponent based on the initial smooth fracture grid to generate a random fracture path containing roughness includes: Recursively perform the following operations on each crack path: Apply random displacement in the vertical direction at the midpoint to generate new points; The line segment formed by the new point and the endpoint is further subdivided, and the displacement is repeatedly applied. It is iterated n times, and the Hurst exponent is used to control the persistence of the roughness to generate a random crack path with roughness. The iteration method is as follows: ; Where △n is the displacement amplitude of each iteration; △0 is the initial maximum displacement amplitude; n is the number of iterations; H is the Hurst index.

[0014] The present invention has at least the following beneficial effects: The present invention systematically constructs a four-field coupled modeling method for enhanced geothermal system heat-fluid-solidification-chemical (THMC) and establishes a dynamic fracture aperture model to solve the problem of traditional neglect of Wchem chemical aperture change. The temperature field is reconstructed by chemical reaction heat drive. Chemical feedback includes R i andKeq(T); Furthermore, an enhanced geothermal system fracture roughness characterization technology based on fractal geometry is proposed. The recursive displacement method is used to control the roughness to generate a real geothermal fracture network, and the coupling relationship between roughness and the four fields is accurately quantified.

[0015] Furthermore, a formula for calculating the recovery degree of geothermal reservoirs was established. Through the fluid enhancement coefficient β and the chemical heat integral term, the systematic deviation caused by the traditional THM model ignoring chemical reactions and fracture geometry details was solved, providing a high-precision method for EGS reservoir performance evaluation and risk management.

[0016] Furthermore, by establishing a simplified fracture geometry model, the authors focused on characterizing the impact of fracture surface roughness on multi-physics field coupling. This model can numerically simulate the spatiotemporal evolution of temperature, seepage, stress, and chemical fields, as well as the extent of recovery during the enhanced geothermal system (EGS) process. Combining the dynamic characterization of fracture roughness with the fully coupled THMC (heat-fluid-solid-chemical) mechanism overcomes the limitations of traditional EGS models, which ignore chemical fields and fracture geometry details. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 This is a flow chart of the multi-field coupling modeling method based on the heat flow solidification of the enhanced geothermal system of the present invention; Figure 2 This is a diagram of a random fracture model of the geothermal system enhanced by the present invention; Figure 3 The multi-physics coupled seepage field distribution diagram of the geothermal system enhanced by the present invention; Figure 4 The multi-physics coupled stress field distribution diagram of the geothermal system enhanced by the present invention; Figure 5 A random fracture grid division diagram for the enhanced geothermal system of the present invention; Figure 6 Enhanced geothermal system multi-physics coupled chemical field distribution map for the present invention; Figure 7 The temperature field distribution diagram of the multi-physics coupling of the geothermal system enhanced by the present invention; Figure 8 This is a diagram of the coupling relationship of heat flow solidification for the enhanced geothermal system exploitation according to the present invention; Figure 9 The fracture diagram of the enhanced geothermal system considering roughness after simplification of the present invention; Figure 10 This is a graph showing the degree of reservoir recovery in the geothermal system enhanced by the present invention; Figure 11 This is a diagram of temperature changes at different locations in the enhanced geothermal system of the present invention.

[0018] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

[0019] In embodiments of the present invention, the term "and / or" describes the association relationship between associated objects, indicating that three possible relationships exist. For example, "A and / or B" can represent three situations: A exists alone, A and B exist simultaneously, and B exists alone. The character " / " generally indicates that the associated objects are in an "or" relationship.

[0020] It should be noted that the terms "first", "second", etc. in the description and claims of the present invention and the above-mentioned drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequence.

[0021] In the embodiments of the present invention, the term "plurality" refers to two or more than two, and other quantifiers are similar.

[0022] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the embodiments of the present invention will be described in detail below with reference to the accompanying drawings. However, it will be understood by those skilled in the art that in the embodiments of the present invention, many technical details are provided to enable the reader to better understand the present invention. However, even without these technical details and the various changes and modifications based on the following embodiments, the technical solutions claimed in the present invention can be implemented. The division of the following embodiments is for convenience of description and should not constitute any limitation on the specific implementation of the present invention. The various embodiments can be combined with each other and referenced to each other under the premise that there is no contradiction.

[0023] Currently, there are few heat-fluid-solidification (THMC) coupling models for enhanced geothermal systems, and based on the heat-fluid-solidification coupling models, less attention is paid to the impact of fracture roughness on multi-physical fields. In the seepage field, rough fractures may cause non-Darcy flow or changes in local flow velocity. In the temperature field, roughness affects the surface area of ​​heat exchange. In the stress field, the irregularity of the fracture may cause stress concentration, thereby affecting the closure and expansion of the fracture. In the chemical field, the rough surface may increase the contact area of ​​mineral reactions and accelerate the dissolution or precipitation process. Therefore, it is very necessary to establish a heat-fluid-solidification-solidification coupling simulation model of the fracture network of the enhanced geothermal system considering the influence of roughness, so as to accurately obtain the law of the temperature transfer process and the degree of recovery.

[0024] The present invention provides a multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system. Figure 1 A schematic diagram of the multi-field coupling modeling method for heat flow solidification of an enhanced geothermal system provided by the present invention is shown.

[0025] In step S100 , two sets of random numbers are generated as the crack center positions, and the random crack endpoint coordinates are determined according to the crack center positions, trace lengths, and directions. An initial smooth crack grid is formed according to the determined random crack endpoint coordinates.

[0026] Specifically, MATLAB can be used to generate two sets of [0,1] random numbers in the simulation area as the horizontal and vertical coordinates of the crack center position. The length of the series is consistent with the number of cracks N. The average trace length of the crack is specified, and a truncated normal distribution is used to constrain the trace length L>L min (L min = 0.1 × average trace length), N trace lengths were sampled from the truncated distribution. To simplify the cracks, the crack directions were assigned to 0° and 90°. The coordinates of the random crack endpoints were determined based on the crack center position, trace length, and direction. The determined random crack endpoints were connected into straight line segments to form an initial smooth crack network.

[0027] In step S200 , based on the initial smooth crack grid, the Hurst exponent is used to control the persistence of the roughness of each crack path to generate a random crack path containing roughness.

[0028] Among them, the Hurst Exponent is an important indicator that describes the long-term memory, self-similarity or trend persistence of a time series. The value range of the Hurst Exponent is usually 0<H<1, and different values ​​correspond to different characteristics of the time series.

[0029] Recursively perform the following operations on each crack path: Apply random displacement in the vertical direction at the midpoint to generate new points; The line segment formed by the new point and the endpoint is further subdivided, and displacement is repeatedly applied for n iterations. The Hurst exponent is used to control the persistence of the roughness, generating a random crack path with roughness. The larger the Hurst exponent value, the smaller the roughness. The perturbation amplitude of the nth iteration must meet △n to avoid excessive crack distortion. H is usually selected as 0.2-0.3.

[0030] The iteration method is as follows: ; Where △n is the displacement amplitude of each iteration; △0 is the initial maximum displacement amplitude; n is the number of iterations; H is the Hurst index.

[0031] In some embodiments, Δ0 is 0.01 m, Δn is 0.0001 m, and H is 0.2, so n is calculated to be 8.

[0032] Furthermore, in order to simplify the cracks, the random crack path of the roughness can be discretized into polylines, each line segment consisting of iteratively generated subdivision points. By calling CAD with MATLAB to convert it into a parameterized two-dimensional line segment, the intersection of the crack with the boundary and other cracks is determined, and the through-stitched isolated cracks are identified and deleted, the cracks can be simplified (e.g. Figure 2 Finally, the generated random crack file is imported into COMSOL (simulation software) in DXF format to realize modeling (as shown in Figure 3 shown); In step S300, a random fracture model of hot dry rock one-injection-two-production is established based on the generated random fracture path containing roughness.

[0033] Specifically, based on the rough fracture grid model (500m×500m) established in step S200, the geometric structure of the injection well and the production well is constructed. To ensure the stability and convergence of the numerical simulation, the injection well and the production well are created in the main seepage channel of the fracture network. The wellbore is represented by a circle with a diameter of 0.1m. Through the Boolean difference operation of COMSOL, the wellbore area is cut off from the reservoir matrix and the fracture network to form a physically continuous flow channel, ensuring the connectivity between the injection and production system and the fracture network (such as Figure 4 、 5 shown).

[0034] In step S400 , material parameters of the random fracture model are set according to the well logging data, and stress field, seepage field, temperature field, and chemical field are added to the random fracture model.

[0035] Specifically, based on experimental and logging data, parameter expressions are predefined through the COMSOL global function interface, and key parameters such as elastic modulus, thermal conductivity, and Poisson's ratio are defined as variable functions, which supports subsequent one-click call and batch modification. In the parameter design module, the initial temperature of the rock mass (150°C), injection temperature (25°C), injection rate (0.0002m 3 / s) and other parameters.

[0036] Example: Elastic modulus expression: -1.1178*10 -5 *X 2 -0.00925*X+18.70626, where X is the depth.

[0037] The model geometry is divided into five independent modules: boundary layer, matrix domain, fracture network, injection well and production well. The modules are separated using selection list work, and stress field, seepage field, temperature field and chemical field are added to the random fracture model.

[0038] In step S500, initial conditions and boundary conditions are set for the stress field, seepage field, temperature field, and chemical field of the random fracture model, and multi-physics field coupling is set.

[0039] Specifically, according to the logging and experimental data, the material parameters of the random fracture model are set, and expressions for the changes of elastic modulus, Poisson's ratio and other parameters with temperature obtained from experiments are added. Based on these parameters and expressions, the initial conditions and boundary conditions are set for the stress field, seepage field, temperature field and chemical field of the random fracture model, and multi-physics field coupling is set.

[0040] In COMSOL, multi-physics field coupling is set up, solid mechanics field and seepage field are coupled into poroelasticity, and temperature field and solid mechanics field are coupled into thermal expansion; the seepage field affects the convective heat transfer process of the temperature field, and the temperature field affects the seepage field by affecting the physical parameters of the fluid. The effective stress in the stress field can change the permeability of the fractures in the reservoir and thus affect the seepage field. The poroelasticity generated in the seepage field affects its temperature field, and the different temperatures at different locations induce thermal stress effects. The chemical reaction field (such as the dissolution of minerals) affects the permeability and porosity of the rock, thereby affecting fluid flow and heat transfer. The temperature field, seepage field and stress field affect the chemical reaction rate.

[0041] (1) Coupling setting of stress field Set up the solid mechanics module, set the model boundary conditions to fixed constraints, set the initial displacement and initial velocity to 0, adopt a linear elastic material and add a thermal stress coupling term, set the maximum and minimum horizontal principal stresses, calculate the fracture displacement based on the pore fluid pressure provided by the seepage field and the temperature provided by the temperature field, update the fracture aperture, and output the volume strain (i.e., ∇·u in formula (1)) to transmit to the seepage field. The stress drives the fracture deformation and changes the permeability. Among them, the permeability directly affects the flow capacity of the seepage field.

[0042] The calculation formula is as follows: ; (1) Where G is the shear modulus, Pa; u1 is the displacement component, m; λ is the Lame constant, Pa; α B is the Biot coefficient, dimensionless; P is the pore fluid pressure, Pa; E is the elastic modulus, Pa; ν is Poisson's ratio, dimensionless; α T is the coefficient of thermal expansion, 1 / K; T is the local average temperature under the combined action of fluid flow and solid heat conduction, unit K; α C is the chemical expansion coefficient, 1 / (mol / m 3 ); C is the chemical concentration field, mol / m 3 ; F is the body force, N.

[0043] (2) Coupling setting of seepage field Set up the seepage field module, set the boundary condition to no flow, set the initial value to formation pressure, set the parameters of the fluid and matrix in the reservoir matrix and porous media and add porous elastic storage, define the fluid in the fracture and the fracture geometry properties, and form a driving pressure drop drive by setting the pressure difference, where the injection well pressure is higher than the formation pressure and the production well pressure is lower than the formation pressure. According to the volume strain ε provided by the stress field V The temperature T provided by the temperature field and the porosity ϕ provided by the chemical field are used to calculate the Darcy velocity u, which is transferred to the temperature field and the dilute species transfer field. The output pore fluid pressure is fed back to the solid mechanics module. The specific calculation formula is as follows: The reservoir flow equation is: ; (2) The seepage equation in the crack is: ; (3) in, ϕ is the porosity, dimensionless; ρ f is the fluid density, kg / m 3 ; c t is the compressibility coefficient, 1 / Pa; P is the pore fluid pressure, Pa; t is time, s; u is the Darcy velocity of the fluid in the reservoir, m / s; α B is the Biot coefficient, dimensionless; ε V is the volume strain, ε V =ε xx +ε yy +ε zz , dimensionless, ε xx is the volume strain in the x direction, ε yy is the volume strain in the y direction, ε zz is the volume strain in the z direction; Q f is the source and sink term of the fluid, kg / (m 3 ·s); ϕ chem is the porosity change caused by chemical reaction, dimensionless; R chem is the reaction rate, mol / (m 3 ·s); w is the effective opening of the crack, m; ; d f is the initial crack opening, m; u n is the normal displacement increment, m; w chem is the change in chemically induced crack width, m; r f is the fluid density, kg / m 3 ; u f is the Darcy velocity of the fluid in the fracture, m / s; M is the molar mass of the substance, kg / mol.

[0044] (3) Coupling setting of temperature field Set up the temperature field module, set the boundary conditions of the reservoir and fractures to thermal insulation, set the initial condition to the initial temperature of the rock mass, set the parameters of the fluid and matrix in the reservoir matrix and porous medium, define the fluid and fracture geometry properties in the fracture, set the temperature of the injection well, calculate the temperature field based on the flow velocity provided by the seepage field, the heat source provided by the chemical field, and the displacement provided by solid mechanics, and transfer it to the seepage field, stress field, and chemical field respectively; The calculation formula is as follows: The temperature field equation in the reservoir is: ; (4) The temperature field equation in the crack is: ; (5) Where, (ρC) eff is the effective heat capacity, J / (m 3 ⋅K); ; ϕ is the porosity, dimensionless; ρ m is the fluid density, kg / m 3 ; C f and C m are the specific heat capacities of the fluid and reservoir, J / (kg·K), respectively; ρ f is the fluid density, kg / m 3 ; T is the local average temperature under the combined action of fluid flow and solid heat conduction, unit K; u is the Darcy velocity of the fluid in the reservoir, m / s; λ eff is the thermal conductivity of rock, W / (m·K); ; λ rock is the thermal conductivity of rock, W / (m·K); λ fis the thermal conductivity of the fluid, W / (m·K); W f is the fluid heat source, W / m 3 ; Z is the reaction heat source, W / m 3 ; w is the effective opening of the crack, m; ; d f is the initial crack opening, m; u n is the normal displacement increment, m; w chem is the change in chemically induced crack width, m; t is time, s; u f is the Darcy velocity of the fluid in the fracture, m / s.

[0045] (4) The coupling settings of the chemical field include: Set up the chemical reaction module, add the following reaction equation, set the initial concentration, calculate the reaction rate according to the following formula (6) and transfer it to the dilute species transfer field. According to formula (7), input the temperature T into the heat source and output the reaction heat source to the temperature field; The reaction equation is: ; Formula (6) and formula (7) are as follows: ; (6) ; (7) Among them, R chem is the reaction rate, mol / (m 3 ·s); k0 is the prefactor, mol / (m 2 ·s); E0 is the activation energy, J / mol; R is the ideal gas constant, 8.314 J / (mol·K); T is the local average temperature under the combined action of fluid flow and solid heat conduction, K; K eq (T) is the temperature-dependent equilibrium constant, dimensionless; Q is the ion activity product, unitless; Z is the reaction heat source, W / m 3 ; R chem，i is the reaction rate of the ith reaction, mol / (m 3 ·s); △H iis the reaction enthalpy change of the i-th reaction, J / mol; i is the reaction category, dimensionless; n is the total number of independent chemical reactions in the system and is dimensionless.

[0046] (5) Coupling of dilute species transfer fields Set up the dilute species transfer module, set the boundary conditions of the reservoir to no flux, set the initial value of the fracturing fluid injected into the reservoir to 0, set the basic parameters of the fluid in the reservoir matrix and porous medium, couple the velocity with the seepage field, set the basic parameters of the fluid and material in the fracture, add a reaction in the fracture, and couple its reaction rate with the reaction rate in the chemical reaction. Set the injection well to concentration flux and the outlet to no flux. According to the flow velocity u provided by the seepage field, the flow rate u is calculated by formula (8). f and the reaction rate R provided by the chemical reaction field chem , calculate the ion concentration and pass it to the chemical reaction module to drive the reaction balance; Formula (8) is as follows: ; (8) in, w is the effective opening of the crack, m; C is the chemical concentration field, mol / m 3 ; t is time, s; u f is the Darcy velocity of the fluid in the fracture, m / s; D f is the molecular diffusion coefficient, m 2 / s; R chem is the reaction rate, mol / (m 3 ·s).

[0047] In step S600, the random fracture model is meshed according to different regions of the random fracture model. The fracture region is meshed using an ultra-fine free triangle meshing method, and the remaining reservoir region is meshed using a refined free quadrilateral meshing method.

[0048] The random fracture model is meshed. According to the different areas of the random fracture model, the fracture area is meshed by ultra-fine free triangle meshing (the maximum unit parameter size is 8m and the minimum unit parameter size is 0.016m), and the rest of the reservoir area is meshed by refined free quadrilateral meshing (the maximum unit parameter size is 16m and the minimum unit parameter size is 0.06m) (e.g. Figure 5 ).

[0049] In step S700, the evolution characteristics of stress field, seepage field, temperature field and chemical field are analyzed, and the recovery degree of the reservoir is calculated through post-processing of surface average value; different positions from the injection well are selected to calculate the temperature change, and the change of fracturing fluid viscosity at different positions is simulated through the relationship between temperature and fracturing fluid viscosity.

[0050] Furthermore, in some embodiments, in order to eliminate the interference of vertical heterogeneity of the reservoir on a single temperature monitoring point, the temperature (T) simulation result and the corresponding unit area of ​​the surface unit are extracted for the target reservoir surface, and the comprehensive average temperature value (T) of the reservoir at a certain moment is obtained by the area averaging algorithm as shown in Formula 9. av (t)), the comprehensive average temperature value at a certain moment (T av (t)) is substituted into the calculation formula 10 as the key parameter to obtain the recovery degree of the reservoir under the current conditions.

[0051] The fracturing fluid flows from the injection well to the deep formation, and the temperature T of a surface unit at time t is obtained at the preset key position points (different fracture nodes at 10m, 50m, and 100m away from the injection well). i (t) is substituted into Formula 11, and the local temperature value of each location is converted into the actual viscosity value of the fracturing fluid at the corresponding location. By analyzing the change trend and actual value of the fracturing fluid viscosity, the flow performance of the fracturing fluid in the reservoir is quantitatively evaluated.

[0052] ;(9) in, T av (t) is the integrated average temperature, K; T i (t) is the temperature of the i-th surface unit at time t; A i is the area of ​​the i-th surface unit; n is the total number of divided units.

[0053] Specifically, after the crack mesh is divided, the solver configuration in COMSOL is performed; the default solver configuration is selected, and the solution mode is configured as transient mode. The simulation time and unit are set according to the actual situation on site, and then the output step size is set according to the level of detail of the physical process being studied to ensure the accuracy of multi-physics field coupling.

[0054] After the solver is configured, the temperature field (such as Figure 8 ), seepage field (such as Figure 9 ), stress fields (such as Figure 10 ) and chemical fields (e.g. Figure 11) evolution characteristics analysis reveals the seepage heat transfer mechanism under the coupling of multiple physical fields; the low-temperature domain in the temperature field mainly diffuses non-uniformly along the fractures near the injection well, and rough fractures enhance local heat transfer; the flow velocity in the seepage field mainly seeps along the fractures, and preferential channels appear, and high roughness reduces the effective permeability; the change of effective stress in the stress field is concentrated in the fracture area of ​​the reservoir, and the thermal stress induced by heat exchange, pore fluid pressure and mineral dissolution jointly promote the change of effective stress; the reaction of the chemical field is mainly the dissolution effect along the fractures and the matrix.

[0055] The formula for calculating the recovery degree of the reservoir is: ; (10) Wherein, η is the degree of recovery, dimensionless; V s is the reservoir volume, m 3 ; ρ s is the reservoir rock density, kg / m 3 ; C s is the heat capacity of reservoir rock, J / (kg·K); T in is the initial temperature of the reservoir, °C; T av (t) is the comprehensive average temperature of the reservoir at a certain moment, °C; V f is the reservoir volume, m 3 ; Z chem is the heat source of chemical reaction, W / m 3 ; T inj is the injection fluid temperature, K.

[0056] ; (11) μ f is the actual viscosity of the fracturing fluid, Pa·s; Obviously, the embodiments described above are only part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, those skilled in the art may make other different forms of changes or modifications without making any creative work, and all of these should fall within the scope of protection of the present invention.

Claims

1. A multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system, characterized by: include: Generate two sets of random numbers as the crack center position, and determine the random crack endpoint coordinates according to the crack center position, trace length and direction, and form an initial smooth crack grid based on the determined random crack endpoint coordinates; Based on the initial smooth crack grid, the roughness persistence of each crack path is controlled by the Hurst exponent to generate a random crack path containing roughness; Based on the generated random fracture paths with roughness, a random fracture model of hot dry rock one-injection-two-production is established; According to the well logging data, the material parameters of the random fracture model are set, and the stress field, seepage field, temperature field, and chemical field are added to the random fracture model; Set initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model, and perform multi-physics field coupling settings; Meshing the random fracture model according to different regions of the random fracture model, meshing the fracture region using an ultra-fine free triangle meshing method, and meshing the remaining reservoir regions using a refined free quadrilateral meshing method; Conduct analysis on the evolution characteristics of stress field, seepage field, temperature field and chemical field, calculate the recovery degree of reservoir through post-processing of surface average value; select different locations from the injection well, calculate the temperature change, and simulate the change of fracturing fluid viscosity at different locations through the relationship between temperature and fracturing fluid viscosity.

2. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 1, characterized in that: In the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physics field coupling setting, the stress field coupling setting includes: Set up the solid mechanics module, set the model boundary conditions to fixed constraints, set the initial displacement and initial velocity to 0, use linear elastic materials and add thermal stress coupling terms, set the maximum horizontal principal stress and minimum horizontal principal stress, calculate the displacement of the crack based on the pore fluid pressure provided by the seepage field and the temperature provided by the temperature field, update the crack aperture, output the volume strain ∇·u and transmit it to the seepage field, and the stress drives the crack deformation, changing the permeability; The calculation formula is as follows: ；（1） Where G is the shear modulus, Pa; u1 is the displacement component, m; λ is the Lame constant, Pa; α B is the Biot coefficient, dimensionless; P is the pore fluid pressure, Pa; E is the elastic modulus, Pa; ν is Poisson's ratio, dimensionless; α T is the coefficient of thermal expansion, 1 / K; T is the local average temperature under the combined action of fluid flow and solid heat conduction, unit K; α C is the chemical expansion coefficient, 1 / (mol / m 3 ); C is the chemical concentration field, mol / m 3 ; F is the body force, N.

3. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 1, characterized in that: In the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physics field coupling setting, the coupling setting of the seepage field includes: Set up the seepage field module, set the boundary condition to no flow, set the initial value to formation pressure, set the parameters of the fluid in the reservoir matrix and porous medium and add porous elastic storage, define the fluid in the fracture and the fracture geometry properties, and form a driving pressure drop drive by setting the pressure difference, where the injection well pressure is higher than the formation pressure and the production well pressure is lower than the formation pressure. According to the volume strain ε provided by the stress field V The temperature provided by the temperature field and the porosity ϕ provided by the chemical field are used to calculate the Darcy velocity and transfer it to the temperature field and the dilute species transfer field. The output pore fluid pressure is fed back to the solid mechanics module. The specific calculation formula is as follows: The reservoir flow equation is: ；（2） The seepage equation in the crack is: ；（3） in, ϕ is the porosity, dimensionless; ρ f is the fluid density, kg / m 3 ; c t is the compressibility coefficient, 1 / Pa; P is the pore fluid pressure, Pa; t is time, s; u is the Darcy velocity of the fluid in the reservoir, m / s; α B is the Biot coefficient, dimensionless; ε V is the volume strain, ε V =ε xx +ε yy +ε zz , dimensionless, ε xx is the volume strain in the x direction, ε yy is the volume strain in the y direction, ε zz is the volume strain in the z direction; Q f is the source and sink term of the fluid, kg / (m 3 ·s); ϕ chem is the porosity change caused by chemical reaction, dimensionless; R chem is the reaction rate, mol / (m 3 ·s); w is the effective opening of the crack, m; ； d f is the initial crack opening, m; u n is the normal displacement increment, m; w chem is the change in chemically induced crack width, m; r f is the fluid density, kg / m 3 ; u f is the Darcy velocity of the fluid in the fracture, m / s; M is the molar mass of the substance, kg / mol.

4. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 1, characterized in that: In the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physics field coupling setting, the coupling setting of the temperature field includes: Set up the temperature field module, set the boundary conditions of the reservoir and fractures to thermal insulation, set the initial condition to the initial temperature of the rock mass, set the parameters of the fluid in the reservoir matrix and porous medium, define the fluid in the fracture and the fracture geometry properties, set the temperature of the injection well, calculate the temperature field based on the flow velocity provided by the seepage field, the heat source provided by the chemical field, and the displacement provided by solid mechanics, and transfer it to the seepage field, stress field, and chemical field respectively; The calculation formula is as follows: The temperature field equation in the reservoir is: ；（4） The temperature field equation in the crack is: ；（5） Where, (ρC) eff is the effective heat capacity, J / (m 3 ⋅K); ； ϕ is the porosity, dimensionless; ρ m is the fluid density, kg / m 3 ; C f and C m are the specific heat capacities of the fluid and reservoir, J / (kg·K), respectively; ρ f is the fluid density, kg / m 3 ; T is the local average temperature under the combined action of fluid flow and solid heat conduction, unit K; u is the Darcy velocity of the fluid in the reservoir, m / s; λ eff is the thermal conductivity of rock, W / (m·K); ； λ rock is the thermal conductivity of rock, W / (m·K); λ f is the thermal conductivity of the fluid, W / (m·K); W f is the fluid heat source, W / m 3 ; Z is the reaction heat source, W / m 3 ; w is the effective opening of the crack, m; ； d f is the initial crack opening, m; u n is the normal displacement increment, m; w chem is the change in chemically induced crack width, m; t is time, s; u f is the Darcy velocity of the fluid in the fracture, m / s.

5. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 1, characterized in that: In the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physical field coupling setting, the coupling setting of the chemical field includes: Set up the chemical reaction module, add the following reaction equation, set the initial concentration, calculate the reaction rate according to the following formula (6) and transfer it to the dilute species transfer field. According to formula (7), input the temperature T into the heat source and output the reaction heat source to the temperature field; The reaction equation is: ； Formula (6) and formula (7) are as follows: ；（6） ；（7） Among them, R chem is the reaction rate, mol / (m 3 ·s); k0 is the prefactor, mol / (m 2 ·s); E0 is the activation energy, J / mol; R is the ideal gas constant, 8.314 J / (mol·K); T is the local average temperature under the combined action of fluid flow and solid heat conduction, K; K eq (T) is the temperature-dependent equilibrium constant, dimensionless; Q is the ion activity product, unitless; Z is the reaction heat source, W / m 3 ; R chem，i is the reaction rate of the ith reaction, mol / (m 3 ·s); △H i is the reaction enthalpy change of the i-th reaction, J / mol; i is the reaction category, dimensionless; n is the total number of independent chemical reactions in the system and is dimensionless.

6. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 5, characterized in that: In the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively and performing multi-physical field coupling setting, the chemical field includes a dilute species transfer field, and the coupling setting of the dilute species transfer field includes: Set up the dilute species transfer module, set the boundary conditions of the reservoir to no flux, set the initial value of the fracturing fluid injected into the reservoir to 0, set the basic parameters of the fluid in the reservoir matrix and porous medium, couple the velocity with the seepage field, set the basic parameters of the fluid and material in the fracture, add a reaction in the fracture, and couple its reaction rate with the reaction rate in the chemical reaction. Set the injection well to concentration flux and the outlet to no flux. According to the flow velocity u provided by the seepage field, the flow rate u is calculated by formula (8). f and the reaction rate R provided by the chemical reaction field chem , calculate the ion concentration and pass it to the chemical reaction module to drive the reaction balance; Formula (8) is as follows: ；（8） in, w is the effective opening of the crack, m; C is the chemical concentration field, mol / m 3 ; t is time, s; u f is the Darcy velocity of the fluid in the fracture, m / s; D f is the molecular diffusion coefficient, m 2 / s; R chem is the reaction rate, mol / (m 3 ·s).

7. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 1, characterized in that: The initial conditions and boundary conditions are set for the stress field, seepage field, temperature field, and chemical field of the random fracture model respectively, and multi-physics field coupling is set, including: In COMSOL, multi-physics field coupling is set up, solid mechanics field and seepage field are coupled into poroelasticity, and temperature field and solid mechanics field are coupled into thermal expansion; the seepage field affects the convective heat transfer process of the temperature field, and the temperature field affects the seepage field by affecting the physical parameters of the fluid. The effective stress in the stress field can change the permeability of the fractures in the reservoir and thus affect the seepage field. The poroelasticity occurring in the seepage field affects its temperature field. The chemical reaction field affects the permeability and porosity of the rock, thereby affecting fluid flow and heat transfer. The temperature field, seepage field and stress field affect the chemical reaction rate.

8. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 1, characterized in that: In the step of calculating the recovery degree of the reservoir; selecting different positions from the injection well, calculating the temperature change, and simulating the change in the viscosity of the fracturing fluid at different positions through the relationship between temperature and fracturing fluid viscosity, the formula for calculating the recovery degree of the reservoir is: ；（9） Where η is the degree of recovery, dimensionless; V s is the reservoir volume, m 3 ; ρ s is the reservoir rock density, kg / m 3 ; C s is the heat capacity of reservoir rock, J / (kg·K); T in is the initial temperature of the reservoir, °C; T av (t) is the comprehensive average temperature of the reservoir at a certain moment, °C; V f is the reservoir volume, m 3 ; Z chem is the heat source of chemical reaction, W / m 3 ; T inj is the injection fluid temperature, K.

9. The multi-field coupling modeling method based on heat flow solidification of enhanced geothermal system according to claim 1, characterized in that: The step of controlling the roughness persistence of each crack path by using the Hurst exponent based on the initial smooth crack grid to generate a random crack path containing roughness includes: Recursively perform the following operations on each crack path: Apply random displacement in the vertical direction at the midpoint to generate new points; The line segment formed by the new point and the endpoint is further subdivided, and the displacement is repeatedly applied. It is iterated n times, and the Hurst exponent is used to control the persistence of the roughness to generate a random crack path with roughness. The iteration method is as follows: ； Where △n is the displacement amplitude of each iteration; △0 is the initial maximum displacement amplitude; n is the number of iterations; H is the Hurst index.

Images