Enhanced geothermal system heat-fluid-solid coupling modeling method

Abstract

The application discloses a kind of based on enhanced geothermal system heat fluidization multi-field coupling modeling method, belong to oil and gas field exploitation technical field, including generating two groups of random numbers as fracture center position, and according to fracture center position, trace length and direction determine random fracture endpoint coordinate, according to the random fracture endpoint coordinate determined, form initial smooth fracture grid;According to the initial smooth fracture grid, the roughness persistence of each fracture path is controlled by Hurst index, to generate random fracture path containing roughness;According to the random fracture path containing roughness generated, establish dry hot rock one injection two mining random fracture model;According to well logging data, set the material parameters of random fracture model, add stress field, seepage field, temperature field, chemical field to random fracture model;The initial conditions and boundary conditions of stress field, seepage field, temperature field and chemical field of random fracture model are set respectively, and multi-physical field coupling is set.

Description

TECHNICAL FIELD

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

[0002] The multi-physical field coupling numerical simulation of the enhanced geothermal system (EGS) is an important research topic. The existing enhanced geothermal system model is limited to thermal-flow (TH), flow-solid (HM) or thermal-flow-solid coupling (THM), and does not consider the dynamic change of porosity and permeability caused by mineral dissolution, ignores the time-varying effect of the chemical field (C), and causes distortion of the dynamic process: firstly, the dissolution of minerals under the temperature gradient will increase the permeability in the short term, secondly, the mineral reaction is often accompanied by an endothermic / exothermic process, and ignoring the chemical reaction heat will lead to the reconstruction of the reservoir temperature field, resulting in overestimation of the heat energy extraction efficiency, and thirdly, the rock mechanics parameters are continuously weakened due to mineral dissolution, which aggravates the risk of shear fracture.

[0003] The above content is only used to assist in understanding the technical solutions of the present application and does not represent the acknowledgement of the above content as prior art. SUMMARY

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

[0005] To achieve the above purpose, the present application provides a heat-flow-solidification multi-field coupling modeling method based on an enhanced geothermal system, which comprises:

[0006] Generate two groups of random numbers as crack center positions, and determine random crack end point coordinates according to the crack center positions, trace length and direction, and form an initial smooth crack grid according to the determined random crack end point coordinates;

[0007] According to the initial smooth crack grid, the roughness persistence of each crack path is controlled by the Hurst index to generate a random crack path containing roughness;

[0008] According to the generated random crack path containing roughness, a random crack model of hot dry rock one injection and two extraction is established;

[0009] According to the logging data, the material parameters of the random crack model are set, and the stress field, seepage field, temperature field and chemical field of the random crack model are added;

[0010] The initial conditions and boundary conditions of the stress field, seepage field, temperature field and chemical field of the random crack model are set respectively, and the multi-physical field coupling is set;

[0011] According to different regions of the random crack model, the random crack model is meshed, the meshing method of super-fine free triangle is adopted to mesh the crack region, and the meshing method of fine free quadrilateral is adopted to mesh the remaining region of the reservoir;

[0012] The evolution characteristics of stress field, seepage field, temperature field and chemical field are analyzed, the recovery degree of the reservoir is calculated through surface average post-processing, the temperature change at different positions from the injection well is calculated, and the viscosity change of the fracturing fluid at different positions is simulated through the relationship between temperature and fracturing fluid viscosity.

[0013] Preferably, in the enhanced geothermal system heat flow solidification multi-field coupling modeling method, in the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field and chemical field of the random crack model respectively, and performing multi-physical field coupling setting, the coupling setting of the stress field comprises:

[0014] A solid mechanics module is set, the model boundary condition is set as a fixed constraint, the initial displacement and initial velocity are set as 0, a linear elastic material is adopted and a thermal stress coupling term is added, the maximum horizontal principal stress and the minimum horizontal principal stress are set, the displacement of the crack is calculated according to the pore fluid pressure provided by the seepage field and the temperature provided by the temperature field, the crack opening is updated, the volume strain is output to the seepage field, the stress drives the crack deformation, and the permeability is changed;

[0015] The calculation formula is as follows:

[0016] ; (1)

[0017] Wherein, G is the shear modulus, Pa;

[0018] u1 is the displacement component, m;

[0019] Lambda is the Lame constant, Pa;

[0020] Alpha B is the Biot coefficient, dimensionless;

[0021] P is the pore fluid pressure, Pa;

[0022] E is the elastic modulus, Pa;

[0023] Nu is the Poisson's ratio, dimensionless;

[0024] Alpha T is the thermal expansion coefficient, 1 / K;

[0025] T is the local average temperature under the joint action of fluid flow and solid heat conduction, unit K;

[0026] Alpha Cis the chemical expansion coefficient, 1 / (mol / m 3 );

[0027] C is the chemical concentration field, mol / m 3 ;

[0028] F is the body force, N.

[0029] Preferably, in the enhanced geothermal system heat solidification multi-field coupling modeling method, the step of setting initial conditions and boundary conditions for the stress field, the seepage field, the temperature field, and the chemical field of the random crack model respectively, and the multi-physical field coupling setting, the coupling setting of the seepage field includes:

[0030] Setting the seepage field module, setting the boundary condition as no flow, setting the initial value as the formation pressure, setting the parameters of the fluid and the matrix in the reservoir matrix and the porous medium and adding the porous elastic storage, defining the fluid in the crack and the crack geometric attribute, forming the driving pressure drop drive by setting the pressure difference, wherein the injection well pressure is higher than the formation pressure, the production well pressure is lower than the formation pressure, calculating the Darcy flow velocity according to the volume strain ε V provided by the stress field, the temperature provided by the temperature field, and the porosity ϕ provided by the chemical field, and outputting the pore fluid pressure to the solid mechanics module;

[0031] The specific calculation formula is as follows:

[0032] The seepage equation of the reservoir is:

[0033] ; (2)

[0034] The seepage equation in the crack is:

[0035] ; (3)

[0036] wherein,

[0037] ϕ is the porosity, dimensionless;

[0038] ρ f is the fluid density, kg / m 3 ;

[0039] c t is the compression coefficient, 1 / Pa;

[0040] P is the pore fluid pressure, Pa;

[0041] t is time, s;

[0042] u is the Darcy flow velocity of the fluid in the reservoir, m / s;

[0043] α Bis the Biot coefficient, dimensionless;

[0044] ε V is the volumetric strain, ε V = ε xx + ε yy + ε zz , dimensionless, ε xx is the volumetric strain in the x direction, ε yy is the volumetric strain in the y direction, ε zz is the volumetric strain in the z direction;

[0045] Q f is the fluid source-sink term, kg / (m 3 ·s);

[0046] ϕ chem is the change in porosity due to chemical reaction, dimensionless;

[0047] R chem is the reaction rate, mol / (m 3 ·s);

[0048] w is the effective aperture of the fracture, m;

[0049] ;

[0050] d f is the initial aperture of the fracture, m;

[0051] u n is the normal displacement increment, m;

[0052] w chem is the change in fracture width induced by chemistry, m;

[0053] r f is the fluid density, kg / m 3 ;

[0054] u f is the fluid Darcy velocity in the fracture, m / s;

[0055] M is the molar mass of the substance, kg / mol.

[0056] Preferably, in the enhanced geothermal system heat-solidification multi-field coupling modeling method, in the step of setting initial conditions and boundary conditions for the stress field, the seepage field, the temperature field, and the chemical field of the random fracture model respectively, and performing multi-physical field coupling setting, the coupling setting of the temperature field comprises:

[0057] The temperature field module is set, the boundary condition reservoir and fracture is set as thermal insulation, the initial condition is set as rock mass initial temperature, the parameters of fluid and matrix in reservoir matrix and porous medium are set, the fluid and fracture geometry in fracture are defined, the temperature of water injection well is set, the temperature field is calculated according to the flow velocity provided by the seepage field, the heat source provided by the chemical field and the displacement provided by the solid mechanics, and is transmitted to the seepage field, the stress field and the chemical field respectively;

[0058] The calculation formula is as follows:

[0059] The temperature field equation in the reservoir is:

[0060] ; (4)

[0061] The temperature field equation in the fracture is:

[0062] ; (5)

[0063] Wherein, (ρC) eff is effective heat capacity, J / (m 3 ·K);

[0064] ;

[0065] φ is porosity, dimensionless;

[0066] ρ m is fluid density, kg / m 3 ;

[0067] C f and C m are specific heat capacity of fluid and reservoir respectively, J / (kg·K);

[0068] ρ f is fluid density, kg / m 3 ;

[0069] T is local average temperature under the joint action of fluid flow and solid heat conduction, unit K;

[0070] u is fluid Darcy flow velocity in the reservoir, m / s;

[0071] λ eff is rock thermal conductivity, W / (m·K);

[0072] ;

[0073] λ rock is rock thermal conductivity, W / (m·K);

[0074] λ f is fluid thermal conductivity, W / (m·K);

[0075] W f is the heat source of fluid, W / m 3 ;

[0076] Z is the reaction heat source, W / m 3 ;

[0077] w is the effective opening of the crack, m;

[0078] ;

[0079] d f is the initial opening of the crack, m;

[0080] u n is the normal displacement increment, m;

[0081] w chem is the amount of change in crack width induced by chemistry, m;

[0082] t is time, s;

[0083] u f is the fluid Darcy flow rate in the crack, m / s.

[0084] Preferably, in the enhanced geothermal system heat solidification multi-field coupling modeling method, the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field and chemical field of the random crack model respectively, and coupling setting of the chemical field includes:

[0085] Setting a chemical reaction module, adding the following reaction equation, setting the initial concentration, calculating the reaction rate according to the following formula (6) and transferring it to the dilute substance transfer field, inputting the temperature T into the heat source according to formula (7), and outputting the reaction heat source to the temperature field;

[0086] The reaction equation is:

[0087] ;

[0088] Formula (6) and formula (7) are as follows:

[0089] ; (6)

[0090] ; (7)

[0091] Where, R chem is the reaction rate, mol / (m 3 ·s);

[0092] k0 is the pre-factor, mol / (m 2 ·s);

[0093] E0 is the activation energy, J / mol;

[0094] R is the ideal gas constant, 8.314 J / (mol·K);

[0095] T is the local average temperature under the combined action of fluid flow and solid heat conduction, K;

[0096] K eq (T) is the temperature-dependent equilibrium constant, dimensionless;

[0097] Q is the ion activity product, unitless;

[0098] Z is the heat source of the reaction, W / m 3 ;

[0099] R chem，i is the reaction rate of the i-th reaction, mol / (m 3 ·s);

[0100] △H i is the reaction enthalpy change of the i-th reaction, J / mol;

[0101] i is the reaction category, dimensionless;

[0102] n is the total number of independent chemical reactions in the system, dimensionless.

[0103] Preferably, in the enhanced geothermal system heat flow solidification multi-field coupling modeling method, in the step of setting initial conditions and boundary conditions for the stress field, seepage field, temperature field, and chemical field of the random crack model respectively, and performing multi-physical field coupling setting, the chemical field includes a rare substance transfer field, and the coupling setting of the rare substance transfer field includes:

[0104] Setting a rare substance transfer module, setting the boundary condition of the reservoir to be no flux, setting the initial value of the injected fracturing fluid in the reservoir to be 0, setting the basic parameters of the fluid in the reservoir matrix and porous medium, coupling the velocity with the seepage field, setting the basic parameters of the fluid and material in the crack, adding a reaction in the crack, the reaction rate of which is coupled with the reaction rate in the chemical reaction, setting the injection well to be a concentration flux, setting the outlet to be no flux, calculating the ion concentration according to the flow rate u f provided by the seepage field and the reaction rate R chem provided by the chemical reaction field through formula (8), transferring to the chemical reaction module to drive the reaction equilibrium;

[0105] Formula (8) is as follows:

[0106] ; (8)

[0107] wherein,

[0108] w is the effective opening of the fracture, m;

[0109] C is the chemical concentration field, mol / m 3 ;

[0110] t is time, s;

[0111] u f is the Darcy velocity of fluid in the fracture, m / s;

[0112] D f is the molecular diffusion coefficient, m 2 / s;

[0113] R chem is the reaction rate, mol / (m 3 ·s)。

[0114] Preferably, in the enhanced geothermal system heat flow solidification multi-field coupling modeling method, the initial conditions and boundary conditions of the stress field, the seepage field, the temperature field and the chemical field of the random fracture model are set respectively, and the multi-physical field coupling is set, including:

[0115] In COMSOL, the multi-physical field coupling is set, the solid mechanics field and the seepage field are coupled into porous elasticity, and the temperature field and the solid mechanics field are coupled into thermal expansion; the seepage field affects the convective heat transfer process of the temperature field, the temperature field affects the seepage field by affecting the physical property parameters of the fluid, the effective stress in the stress field can change the permeability of the fracture in the reservoir and then affect the seepage field, the porous elasticity of the seepage field affects its temperature field, and the chemical reaction field affects the permeability and porosity of the rock, so as to affect the fluid flow and heat transfer, and the temperature field, the seepage field and the stress field affect the chemical reaction rate.

[0116] Preferably, in the enhanced geothermal system heat flow solidification multi-field coupling modeling method, the recovery degree of the reservoir is calculated; in the step of selecting different positions from the injection well, calculating the temperature change, and simulating the change of the viscosity of the fracturing fluid at different positions through the relationship between the temperature and the viscosity of the fracturing fluid, the recovery degree of the reservoir is calculated. The formula is:

[0117] ; (9)

[0118] Wherein, η is the recovery degree, dimensionless;

[0119] V s is the volume of the reservoir, m 3 ;

[0120] ρ s is the density of the reservoir rock, kg / m 3 ;

[0121] C sThe heat capacity of the reservoir rock is expressed in J / (kg·K).

[0122] T in The initial temperature of the reservoir is ℃;

[0123] T av (t) represents the overall average temperature of the reservoir at a certain moment, in °C;

[0124] V f Let m be the reservoir volume. 3 ;

[0125] Z chem As a heat source for chemical reactions, W / m 3 ;

[0126] T inj The temperature of the injected fluid is K.

[0127] Preferably, in the enhanced geothermal system heat flow solidification multi-field coupling modeling method, the step of generating random crack paths containing roughness by controlling the roughness persistence of each crack path through the Hurst exponent based on the initial smooth crack mesh includes:

[0128] For each crack path, recursively perform the following operations:

[0129] A random displacement is applied to the midpoint along the vertical direction to generate a new point;

[0130] The line segment formed by the new point and the endpoint is further subdivided, and the displacement is repeatedly applied. This process is repeated n times. The roughness persistence is controlled by the Hurst exponent to generate random crack paths containing roughness. The iteration method is as follows:

[0131] ;

[0132] Where Δn is the displacement magnitude in each iteration;

[0133] △0 represents the initial maximum displacement amplitude;

[0134] n is the number of iterations;

[0135] H is the Hurst exponent.

[0136] The present invention has at least the following beneficial effects:

[0137] This invention systematically constructs a four-field coupled modeling method for enhanced geothermal systems (THMC), establishes a dynamic fracture aperture model, and solves the problem of traditional methods neglecting Wchem chemical aperture changes. It reconstructs the temperature field through chemical reaction heat-driven temperature field reconstruction, with chemical feedback including R... i andKeq(T);

[0138] Further, a fractal geometry-based enhanced geothermal system fracture roughness representation technology is proposed, and a recursive displacement method is used to control roughness to generate a real geothermal fracture network, and the coupling relationship between roughness and four fields is accurately quantified.

[0139] Further, a geothermal reservoir recovery degree calculation formula is established, and through the fluid enhancement coefficient beta and the chemical heat integral term, the systematic deviation caused by the traditional THM model ignoring chemical reaction and fracture geometric details is solved, providing a high-precision method for EGS reservoir performance evaluation and risk management.

[0140] Further, by establishing a simplified fracture geometry model, the influence of fracture surface roughness on multi-physical field coupling is represented, which can simulate the spatio-temporal distribution evolution law and recovery degree of temperature field, seepage field, stress field and chemical field during the exploitation of enhanced geothermal system. The dynamic representation of fracture roughness and the THMC (thermal-fluid-solid-chemical) full coupling mechanism are combined, which breaks through the limitation of traditional EGS model ignoring chemical field and fracture geometric details. BRIEF DESCRIPTION OF DRAWINGS

[0141] Figure 1 The flowchart of the enhanced geothermal system heat flow solidification multi-field coupling modeling method based on the present application;

[0142] Figure 2 The random fracture model of the enhanced geothermal system of the present application;

[0143] Figure 3 The enhanced geothermal system multi-physical field coupling seepage field distribution diagram of the present application;

[0144] Figure 4 The enhanced geothermal system multi-physical field coupling stress field distribution diagram of the present application;

[0145] Figure 5 The enhanced geothermal system random fracture grid division diagram of the present application;

[0146] Figure 6 The enhanced geothermal system multi-physical field coupling chemical field distribution diagram of the present application;

[0147] Figure 7 The enhanced geothermal system multi-physical field coupling temperature field distribution diagram of the present application;

[0148] Figure 8 The enhanced geothermal system exploitation heat flow solidification coupling relationship diagram of the present application;

[0149] Figure 9 The enhanced geothermal system fracture diagram of the present application considering roughness after simplification;

[0150] Figure 10 The enhanced geothermal system reservoir recovery degree curve diagram of the present application;

[0151] Figure 11 Figure 1 is a temperature variation diagram of the enhanced geothermal system at different positions of the present application.

[0152] The object, functional features and advantages of the present application will be further described with reference to the embodiments and the accompanying drawings. DETAILED DESCRIPTION

[0153] In the embodiments of the present application, the term "and / or" describes the association relationship of the associated objects, which means that there can be three relationships, for example, A and / or B can represent the three cases of A existing alone, A and B existing simultaneously, and B existing alone. The character " / " generally represents an "or" relationship between the front and rear associated objects.

[0154] It should be noted that the terms "first", "second", and the like in the specification and claims of the present application and the above-mentioned drawings are used to distinguish similar objects, and do not necessarily have to describe a specific order or sequence.

[0155] In the embodiments of the present application, the term "a plurality of" means two or more, and other quantifiers are similar.

[0156] In order to make the object, technical scheme and advantages of the embodiments of the present application more clear, the embodiments of the present application will be described in detail below with reference to the drawings. However, those skilled in the art can understand that in the embodiments of the present application, many technical details are proposed in order to make the reader better understand the present application. However, even without these technical details and various changes and modifications based on the following embodiments, the technical scheme claimed by the present application can be implemented. The following embodiments are divided for the convenience of description, and should not constitute any limitation on the specific implementation of the present application, and the embodiments can be combined and referred to each other without contradiction.

[0157] At present, there are few heat-flow-solidification (THMC) coupling models for enhanced geothermal systems, and less attention is paid to the influence of roughness of cracks on multi-physical fields on the basis of the heat-flow-solidification coupling model. In the seepage field, rough cracks can cause non-Darcy flow or changes in local flow rate, in the temperature field, roughness affects the heat exchange surface area, in the stress field, the irregularity of cracks can cause stress concentration, thereby affecting the closure and expansion of cracks, and in the chemical field, rough surfaces can increase the contact area of mineral reactions, accelerate the dissolution or precipitation process. Therefore, it is very necessary to establish a heat-flow-solidification multi-field coupling modeling method for enhanced geothermal system crack network considering the influence of roughness, and then accurately obtain the temperature transfer process rule and the recovery degree.

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

[0159] Two groups of random numbers are generated as crack center positions at step S100, and random crack endpoint coordinates are determined according to the crack center positions, trace lengths and directions, and an initial smooth crack grid is formed according to the determined random crack endpoint coordinates.

[0160] Specifically, two groups of [0, 1] random numbers can be generated as the horizontal and vertical coordinates of the crack center positions in the simulation area by MATLAB, and the length of the number sequence is consistent with the number of cracks N; the average trace length of the cracks is specified, a truncated normal distribution is used, and the trace length L is constrained to be L min (L min =0.1×average trace length), N trace lengths are extracted from the truncated distribution; to simplify the cracks, the crack direction is specified as 0° and 90°, and the random crack endpoint coordinates are determined according to the crack center positions, trace lengths and directions. According to the determined random crack endpoint coordinates, the straight line segments are connected to form an initial smooth crack network.

[0161] At step S200, according to the initial smooth crack grid, the roughness persistence of each crack path is controlled by the Hurst index to generate a random crack path containing roughness.

[0162] The Hurst index (Hurst Exponent) is an important indicator for describing the long-term memory, self-similarity or trend persistence of a time series, and the value range of the Hurst index is usually 0

[0163] The following operations are recursively performed for each crack path:

[0164] A random displacement is applied to the midpoint in the vertical direction to generate a new point;

[0165] The line segment formed by the new point and the endpoint is further subdivided, and the displacement is repeatedly applied, and the iteration is performed n times, and the roughness persistence is controlled by the Hurst index to generate a random crack path containing roughness. The greater the Hurst index value, the smaller the roughness, and the iteration disturbance amplitude of the nth time needs to satisfy Δn to avoid excessive distortion of the crack, and H is usually selected to be 0.2-0.3.

[0166] The iteration method is as follows:

[0167] ;

[0168] Where Δn is the displacement amplitude of each iteration;

[0169] Δ0 is the initial maximum displacement amplitude;

[0170] n is the number of iterations;

[0171] H is the Hurst index.

[0172] In some embodiments, △0 is 0.01 m, △n is 0.0001 m, H is 0.2, and thus n is 8.

[0173] Further, in order to facilitate the simplification of the cracks, the random crack path of roughness can be discretized into multiple line segments, each line segment being composed of subdivision points generated by iteration. By calling CAD through MATLAB to convert into a parameterized two-dimensional line segment, the intersection of the crack with the boundary and other cracks is determined, the through joint isolated cracks are identified and deleted, and thus the crack is simplified (as shown in Figure 2 , finally the generated random crack file is imported into COMSOL (simulation software) in the format of DXF to realize modeling (as shown in Figure 3 ).

[0174] At step S300, a dry hot rock one-injection and two-production random crack model is established according to the generated random crack path containing roughness.

[0175] Specifically, the geometric structure of the injection well and the production well is constructed on the basis of the rough crack grid model established in step S200 (500 m x 500 m). In order to ensure the stability and convergence of numerical simulation, the injection well and the production well are created at the main seepage channel of the crack network, and the wellbore is represented by a circle with a diameter of 0.1 m. Through the Boolean difference set operation of COMSOL, the wellbore area is cut from the reservoir matrix and the crack network to form a physically continuous flow channel, and the communication between the injection-production system and the crack network is ensured (as shown in Figure 4 , 5 ).

[0176] At step S400, the material parameters of the random crack model are set according to the logging data, and the stress field, the seepage field, the temperature field and the chemical field are added to the random crack model.

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

[0178] For example, the expression of elastic modulus is: -1.1178*10 -5 *X 2 -0.00925*X+18.70626, X is the depth.

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

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

[0181] Specifically, according to the logging and experimental data, the material parameters of the random fracture model are set. Expressions of parameters such as elastic modulus and Poisson's ratio varying with temperature are added according to experiments. Based on these parameters and expressions, initial conditions and boundary conditions are set for the stress field, the seepage field, the temperature field and the chemical field of the random fracture model, and multi-physical field coupling is set.

[0182] In COMSOL, the multi-physical field coupling is set. The solid mechanics field and the seepage field are coupled into porous elasticity, and the temperature field and the solid mechanics field are coupled into thermal expansion. The seepage field affects the convective heat transfer process of the temperature field. The temperature field affects the seepage field by affecting the physical property 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 porous elasticity of the seepage field affects its temperature field. Thermal stress is induced due to different temperatures at different positions. The chemical reaction field (such as mineral dissolution) affects the permeability and porosity of the rock, so as to affect the fluid flow and heat transfer. The temperature field, the seepage field and the stress field affect the chemical reaction rate.

[0183] (1) Coupling setting of stress field

[0184] The solid mechanics module is set. The model boundary conditions are set as fixed constraints. The initial displacement and the initial velocity are set as 0. Linear elastic material is adopted and a thermal stress coupling term is added. The maximum horizontal principal stress and the minimum horizontal principal stress are set. According to the pore fluid pressure provided by the seepage field and the temperature provided by the temperature field, the displacement of the fracture is calculated. The fracture opening is updated. The volumetric strain (i.e., ∇·u in formula (1)) is output to the seepage field. The stress drives the deformation of the fracture and changes the permeability. The permeability directly affects the flow capacity of the seepage field.

[0185] The calculation formula is as follows:

[0186] (1)

[0187] Wherein, G is the shear modulus, Pa;

[0188] u1 is the displacement component, m;

[0189] λ is the Lame constant, Pa;

[0190] αB Biot coefficient, dimensionless;

[0191] P is pore fluid pressure, Pa;

[0192] E is elastic modulus, Pa;

[0193] ν is Poisson's ratio, dimensionless;

[0194] α T is thermal expansion coefficient, 1 / K;

[0195] T is local average temperature under the combined action of fluid flow and solid heat conduction, unit K;

[0196] α C is chemical expansion coefficient, 1 / (mol / m 3 );

[0197] C is chemical concentration field, mol / m 3 ;

[0198] F is body force, N.

[0199] (2) Coupling of seepage field

[0200] Set the boundary condition of the seepage field module to no flow, and the initial value to the formation pressure. Set the parameters of the fluid and matrix in the reservoir matrix and porous medium and add porous elastic storage. Define the fluid in the fracture and the fracture geometry properties. Form the driving pressure drop 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 ε V provided by the stress field, the temperature T provided by the temperature field, and the porosity ϕ provided by the chemical field, the Darcy velocity u is calculated and transferred to the temperature field and the dilute substance transfer field. The output of the pore fluid pressure is fed back to the solid mechanics module.

[0201] The specific calculation formula is as follows:

[0202] The seepage equation of the reservoir is:

[0203] ; (2)

[0204] The seepage equation in the fracture is:

[0205] ; (3)

[0206] where,

[0207] ϕ is porosity, dimensionless;

[0208] ρ f is the fluid density, kg / m 3 ;

[0209] c t is the compressibility, 1 / Pa;

[0210] P is the pore fluid pressure, Pa;

[0211] t is time, s;

[0212] u is the fluid Darcy velocity in the reservoir, m / s;

[0213] α B is the Biot’s coefficient, dimensionless;

[0214] ε V is the volumetric strain, ε V = ε xx + ε yy + ε zz , dimensionless, ε xx is the volumetric strain in the x direction, ε yy is the volumetric strain in the y direction, ε zz is the volumetric strain in the z direction;

[0215] Q f is the fluid source / sink term, kg / (m 3 ·s);

[0216] ϕ chem is the change in porosity due to chemical reaction, dimensionless;

[0217] R chem is the reaction rate, mol / (m 3 ·s);

[0218] w is the effective aperture of the fracture, m;

[0219] ;

[0220] d f is the initial aperture of the fracture, m;

[0221] u n is the normal displacement increment, m;

[0222] w chem is the change in fracture width induced by the chemical, m;

[0223] r f is the fluid density, kg / m 3 ;

[0224] u f is the fluid Darcy velocity in the fracture, m / s;

[0225] M is the molar mass of the substance, kg / mol.

[0226] (3) Coupling of temperature field

[0227] The temperature field module is set, the boundary condition reservoir and fracture is set as thermal insulation, the initial condition is set as the initial temperature of rock mass, the parameters of fluid and matrix in reservoir matrix and porous medium are set, the fluid and fracture geometry in fracture are defined, the temperature of water injection well is set, the temperature field is calculated according to the flow rate provided by the seepage field, the heat source provided by the chemical field and the displacement provided by the solid mechanics, and is transmitted to the seepage field, stress field and chemical field respectively;

[0228] The calculation formula is as follows:

[0229] The temperature field equation in the reservoir is:

[0230] ; (4)

[0231] The temperature field equation in the fracture is:

[0232] ; (5)

[0233] Wherein, (ρC) eff is the effective heat capacity, J / (m 3 ·K);

[0234] ;

[0235] φ is the porosity, dimensionless;

[0236] ρ m is the fluid density, kg / m 3 ;

[0237] C f and C m are the specific heat capacity of fluid and reservoir respectively, J / (kg·K);

[0238] ρ f is the fluid density, kg / m 3 ;

[0239] T is the local average temperature under the joint action of fluid flow and solid heat conduction, unit K;

[0240] u is the fluid Darcy flow rate in the reservoir, m / s;

[0241] λ eff is the rock thermal conductivity, W / (m·K);

[0242] ;

[0243] λ rock is the rock thermal conductivity, W / (m·K);

[0244] λf is the fluid thermal conductivity, W / (m·K);

[0245] W f As a fluid heat source, W / m 3 ;

[0246] Z is the heat source for the reaction, W / m 3 ;

[0247] w represents the effective crack aperture, in meters;

[0248] ;

[0249] d f Let m be the initial crack aperture.

[0250] u n Let m be the normal displacement increment.

[0251] w chem The value is the change in chemically induced crack width, expressed in meters (m).

[0252] t represents time, in seconds;

[0253] u f The Darcy velocity of the fluid in the crack is denoted as m / s.

[0254] (4) The coupling setup of the chemical field includes:

[0255] Set up a chemical reaction module, add the following reaction equation, set the initial concentration, calculate the reaction rate transfer to the dilute substance transfer field according to the following formula (6), input the temperature T into the heat source according to the formula (7), and output the reaction heat source to the temperature field.

[0256] The reaction equation is:

[0257] ;

[0258] Formulas (6) and (7) are as follows:

[0259] (6)

[0260] (7)

[0261] Among them, R chem The reaction rate is expressed in mol / (m²). 3 ·s);

[0262] k0 is the pre-factor, mol / (m 2 ·s);

[0263] E0 is the activation energy, J / mol;

[0264] R is the ideal gas constant, 8.314 J / (mol·K);

[0265] T is the local average temperature under the combined action of fluid flow and solid heat conduction, K;

[0266] K eq (T) is the temperature-dependent equilibrium constant, dimensionless;

[0267] Q is the ion activity product, unitless;

[0268] Z is the heat source of reaction, W / m 3 ;

[0269] R chem，i is the reaction rate of the i-th reaction, mol / (m 3 ·s);

[0270] △H i is the reaction enthalpy change of the i-th reaction, J / mol;

[0271] i is the reaction category, dimensionless;

[0272] n is the total number of independent chemical reactions in the system, dimensionless.

[0273] (5) Coupling of dilute species transfer field

[0274] Set the boundary condition of the reservoir to be no flux, the initial value of the injected fracturing fluid in the reservoir is 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, the reaction rate of which is coupled with the reaction rate in the chemical reaction, set the injection well to be a concentration flux, and the outlet to be no flux, according to the flow rate u f provided by the seepage field and the reaction rate R chem provided by the chemical reaction field, calculate the ion concentration by formula (8), transfer it to the chemical reaction module, and drive the reaction equilibrium;

[0275] Formula (8) is as follows:

[0276] ; (8)

[0277] wherein,

[0278] w is the effective opening of the fracture, m;

[0279] C is the chemical substance concentration field, mol / m 3 ;

[0280] t is time, s;

[0281] u fis the Darcy velocity of fluid in the fracture, m / s;

[0282] D f is the molecular diffusion coefficient, m 2 / s;

[0283] R chem is the reaction rate, mol / (m 3 ·s).

[0284] At step S600, the random fracture model is meshed according to different regions of the random fracture model, the meshing method of super-fine free triangle is used to mesh the fracture region, and the meshing method of fine free quadrilateral is used to mesh the remaining region of the reservoir.

[0285] The dry random fracture model is meshed, the fractures are refined and super-refined according to different regions of the random fracture model, the meshing method of super-fine free triangle is used to mesh the fracture region (the maximum unit parameter size is 8 m and the minimum unit parameter size is 0.016 m), and the meshing method of fine free quadrilateral is used to mesh the remaining region of the reservoir (the maximum unit size is 16 m and the minimum unit size is 0.06 m) (as shown in FIG. 6). Figure 5 。

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

[0287] Further, in some embodiments, in order to eliminate the interference of the vertical heterogeneity of the reservoir on a single temperature monitoring point, the temperature (T) simulation result of the surface unit and the corresponding unit area are extracted for the target reservoir surface, the integrated average temperature value (T av (t)) of the reservoir at a certain moment is obtained through the area average algorithm, the obtained integrated average temperature value (T av (t)) at a certain moment is used as a key parameter to substitute into the calculation formula 10, and the recovery degree of the reservoir under the current condition is obtained.

[0288] The fracturing fluid flows from the injection well to the deep formation, and the temperature T i(t), into formula 11, the local temperature value of each position point obtained is converted into the actual viscosity value of the corresponding position fracturing fluid, and by analyzing the change trend and actual value size of the fracturing fluid viscosity, the flow performance of the fracturing fluid in the reservoir is quantitatively evaluated.

[0289] ;(9)

[0290] wherein,

[0291] T av (t) is the comprehensive average temperature, K;

[0292] T i (t) is the temperature of the i-th surface unit at time t;

[0293] A i is the area of the i-th surface unit;

[0294] n is the total number of divided units.

[0295] Specifically, after the fracture grid is divided, the solver configuration in COMSOL is carried out; the default solver configuration is selected, the solving 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 is set according to the detailed degree of the physical process of the research, so as to ensure the accuracy of the multi-physical field coupling.

[0296] After the solver configuration is completed, the evolution characteristics of the temperature field (such as Figure 8 ), the seepage field (such as Figure 9 ), the stress field (such as Figure 10 ) and the chemical field (such as Figure 11 ) are analyzed, so as to reveal the seepage heat exchange mechanism under the coupling of multi-physical fields; in the low temperature field, the non-uniform diffusion along the fracture near the injection well is mainly from the injection well, and the rough fracture enhances the local heat exchange; the flow velocity of the seepage field mainly flows along the fracture, and the preferential channel appears, and the high roughness reduces the effective permeability; the change of the effective stress in the stress field is concentrated in the fracture area in the reservoir, and the heat exchange induced thermal stress, pore fluid pressure and mineral dissolution jointly cause the change of the effective stress; the reaction in the chemical field is mainly the dissolution effect along the fracture and the matrix.

[0297] The formula for calculating the recovery degree of the reservoir is:

[0298] ; (10)

[0299] wherein, η is the recovery degree, dimensionless;

[0300] V s is the volume of the reservoir, m 3 ;

[0301] ρs ρs is the reservoir rock density, kg / m 3 ;

[0302] C s Cp is the reservoir rock heat capacity, J / (kg·K);

[0303] T in T0 is the initial temperature of the reservoir, ℃;

[0304] T av T(t) is the integrated average temperature of the reservoir at a certain time, ℃;

[0305] V f V is the reservoir volume, m 3 ;

[0306] Z chem Z is the chemical reaction heat source, W / m 3 ;

[0307] T inj Tin is the injection fluid temperature, K.

[0308] ; (11)

[0309] μ f μ is the actual viscosity value of the fracturing fluid, Pa·s;

[0310] Obviously, the above described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, other different forms of changes or variations can be made by those skilled in the art without making creative efforts, and all should belong to the protection scope of the present application.

Claims

1. A method for multi-field coupling modeling of thermal solidification based on an enhanced geothermal system, characterized in that, The method 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 length and direction, and forming an initial smooth crack grid according to the determined random crack end point coordinates; controlling roughness persistence of each crack path through a Hurst index according to the initial smooth crack grid to generate a random crack path containing roughness; establishing a random crack model of hot dry rock one injection and two production according to the generated random crack path containing roughness; setting material parameters of the random crack model according to logging data, adding a stress field, a seepage field, a temperature field and a chemical field to the random crack model; setting initial conditions and boundary conditions of the stress field, the seepage field, the temperature field and the chemical field of the random crack model respectively, and performing multi-physical field coupling setting; dividing the random crack model into grids according to different regions of the random crack model, and adopting a super-fine free triangle grid division method to divide the crack region into grids and adopting a fine free quadrilateral grid division method to divide the remaining region of the reservoir into grids; carrying out stress field, seepage field, temperature field and chemical field evolution characteristic analysis, calculating the recovery degree of the reservoir through surface average post-processing, selecting different positions from the injection well to calculate temperature changes, and simulating the viscosity changes of the fracturing fluid at different positions through the relationship between temperature and the viscosity of the fracturing fluid; the setting initial conditions and boundary conditions of the stress field, the seepage field, the temperature field and the chemical field of the random crack model respectively, and performing multi-physical field coupling setting, comprises: in COMSOL, the multi-physical field coupling setting is performed, the solid mechanics field and the seepage field are coupled into porous elasticity, and the temperature field and the solid mechanics field are coupled into thermal expansion; the seepage field affects the convection heat transfer process of the temperature field, the temperature field affects the seepage field by affecting the physical property parameters of the fluid, the effective stress in the stress field can change the permeability of the crack in the reservoir and then affect the seepage field, the porous elasticity of the seepage field affects its temperature field, and the chemical reaction field affects the permeability and porosity of the rock, so as to affect the fluid flow and heat transfer, and the temperature field, the seepage field and the stress field affect the chemical reaction rate; in the step of setting initial conditions and boundary conditions of the stress field, the seepage field, the temperature field and the chemical field of the random crack model respectively, and performing multi-physical field coupling setting, the coupling setting of the seepage field comprises: The seepage field module is set, the boundary condition is set as no flow, the initial value is set as the formation pressure, the parameters of the fluid in the reservoir matrix and the porous medium are set and the porous elastic storage is added, the fluid in the cracks and the crack geometric properties are defined, the differential pressure is set to form the driving pressure drop drive, wherein the injection well pressure is higher than the formation pressure, the production well pressure is lower than the formation pressure, the Darcy flow velocity is calculated according to the volume strain ε V provided by the stress field, the temperature provided by the temperature field and the porosity φ provided by the chemical field, the Darcy flow velocity is transmitted to the temperature field and the rare substance transmission field, and the output pore fluid pressure is fed back to the solid mechanics module; the specific calculation formula is as follows: the seepage equation of the reservoir is: ；（2） the seepage equation in the crack is: ；（3） wherein, phi is the porosity, dimensionless; p f p is the fluid density, kg / m 3 ; c t Compressibility factor, 1 / Pa; P is the pore fluid pressure, Pa; t is the time, s; u is the fluid Darcy velocity in the reservoir, m / s; α B Biot's coefficient, dimensionless; ε V For volumetric strain, ε V =ε xx +ε yy +ε zz , dimensionless, ε xx Let ε be the volumetric strain in the x-direction. yy Let ε be the volumetric strain in the y-direction. zz The volumetric strain is in the z-direction; Q f kg / (m 3 ·s) ϕ chem Dimensionless for the change in porosity due to chemical reactions. R chem for the reaction rate, mol / (m 3 ·s); w is the effective opening of the crack, m; ； d f for the initial opening of the crack, m; u n is the normal displacement increment, m; w chem is the chemically induced crack width change, m; r f for fluid density, kg / m 3 ; u f Cfis the fracture fluid's critical velocity, m / s; M is the molar mass of the substance, kg / mol; the coupling setting of the chemical field comprises: setting a chemical reaction module, adding the following reaction equation, setting the initial concentration, calculating the reaction rate according to the formula (6) and transferring it to the dilute substance transfer field, inputting the temperature T into the heat source according to the formula (7), and outputting the reaction heat source to the temperature field; the reaction equation is: ； the formula (6) and the formula (7) are as follows: ；（6） ；（7） wherein R chem is the reaction rate, mol / (m 3 ·s); k0 is a pre-factor, mol / (m 2 ·s); E0 is the activation energy, J / mol; R is an ideal gas constant, 8.314 J / (mol·K); T is the local average temperature under the joint 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 a source of reaction heat, W / m 3 ; R chem，i Rate of reaction for the ith reaction, mol / (m 3 ·s); ΔH i ΔH is the enthalpy change for the ith reaction, J / mol; i is the reaction category, dimensionless; n is the total number of independent chemical reactions in the system, dimensionless; The chemical field includes a rare substance transfer field, and the coupling setting of the rare substance transfer field comprises: The thin substance transfer module is set, the boundary condition of the reservoir is set as no flux, the initial value of the fracturing fluid injected in the reservoir is 0, the basic parameters of the fluid in the reservoir matrix and the porous medium are set, the velocity is coupled with the seepage field, the basic parameters of the fluid and the material of the fracture are set, the reaction is added in the fracture, the reaction rate is coupled with the reaction rate in the chemical reaction, the injection well is set as the concentration flux, the outlet is set as no flux, the ion concentration is calculated according to the flow rate u f and the reaction rate R chem provided by the seepage field and the chemical reaction field through formula (8), is transferred to the chemical reaction module, and the reaction equilibrium is driven; Formula (8) is as follows: ；（8） wherein, w is the effective opening of the fracture, m; C is the chemical concentration field, mol / m 3 ; t is time, s; u f Cfis the fracture fluid's critical velocity, m / s; D f D is the diffusion coefficient of the solute in the solvent, m2 / s; 2 D is the diffusion coefficient of the solute in the solvent, m2 / s; R chem for the reaction rate, mol / (m 3 ·s); The step of generating a random fracture path containing roughness according to the initial smooth fracture grid, by controlling roughness persistence through a Hurst index for each fracture path, comprises: The following operations are recursively performed for each fracture path: A random displacement is applied in the vertical direction at the midpoint to generate a new point; The line segment formed by the new point and the end point is further subdivided, the displacement is repeatedly applied, and the iteration is performed n1 times, so as to generate a random fracture path containing roughness by controlling roughness persistence through a Hurst index, and the iteration mode is as follows: ； wherein, Δn1 is the displacement amplitude of each iteration; Δ0 is the initial maximum displacement amplitude; n1 is the number of iterations; H is the Hurst index.

2. The modeling method of claim 1, wherein, In the step of setting initial conditions and boundary conditions for the stress field, the seepage field, the temperature field and the chemical field of the random fracture model respectively, and performing multi-physical field coupling setting, the coupling setting of the stress field comprises: A solid mechanics module is set, the model boundary conditions are set as fixed constraints, the initial displacement and initial velocity are set as 0, a linear elastic material is adopted and a thermal stress coupling term is added, the maximum horizontal principal stress and the minimum horizontal principal stress are set, the displacement of the fracture is calculated according to the pore fluid pressure provided by the seepage field and the temperature provided by the temperature field, the fracture opening is updated, and the volume strain ∇·u is output to the seepage field, the stress drives the deformation of the fracture, and the permeability is changed; The calculation formula is as follows: ；（1） wherein, G is the shear modulus, Pa; u1 is the displacement component, m; λ is the Lame constant, Pa; α B Biot's coefficient, dimensionless; P is the pore fluid pressure, Pa; E is the elastic modulus, Pa; ν is the Poisson's ratio, dimensionless; α T is the thermal expansion coefficient, 1 / K; T is the local average temperature under the joint action of fluid flow and solid heat conduction, 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 modeling method of claim 1, wherein, In the step of setting initial conditions and boundary conditions for the stress field, the seepage field, the temperature field and the chemical field of the random fracture model respectively, and performing multi-physical field coupling setting, the coupling setting of the temperature field comprises: A temperature field module is set, the boundary conditions of the reservoir and the fracture are set as thermal insulation, the initial conditions are set as the initial temperature of the rock mass, the parameters of the fluid in the reservoir matrix and the porous medium are set, the fluid in the fracture and the geometric properties of the fracture are defined, the temperature of the water injection well is set, the temperature field is calculated according to the flow rate provided by the seepage field, the heat source provided by the chemical field and the displacement provided by the solid mechanics, and is respectively transmitted to the seepage field, the stress field and the chemical field; The calculation formula is as follows: The temperature field equation in the reservoir is: ；（4） The temperature field equation in the fracture is: ；（5） where (pC) eff is the effective heat capacity, J / (m 3 · K); ； wherein, p m p 3 ; C f and C m Cp and Cp are the specific heat capacities of the fluid and reservoir, respectively, J / (kg·K); p f p is the fluid density, kg / m 3 ; ϕ is the porosity, dimensionless; T is the local average temperature under the joint action of fluid flow and solid heat conduction, K; λ eff λ is the thermal conductivity of the rock, W / (m K); ； λ rock λ is the thermal conductivity of the rock, W / (m K); λ f Fluid thermal conductivity, W / (m K); W f W / m2K 3 ; Z is a source of reaction heat, W / m 3 ; u is the Darcy flow rate of the fluid in the reservoir, m / s; ； d f for the initial crack opening, m; u n is the normal displacement increment, m; w chem is the chemically induced crack width change, m; w is the effective opening of the fracture, m; t is time, s; u f C is the critical velocity, m / s.

4. The modeling method of claim 1, wherein, The method comprises the following steps: calculating the recovery degree of the reservoir; selecting different positions from the injection well; calculating the temperature change; and simulating the change of the viscosity of the fracturing fluid at different positions by the relationship between the temperature and the viscosity of the fracturing fluid. ；（9） In the step of calculating the recovery degree of the reservoir in the step of simulating the change of the viscosity of the fracturing fluid at different positions by the relationship between the temperature and the viscosity of the fracturing fluid, the formula for calculating the recovery degree of the reservoir is as follows: wherein η is the recovery degree, and dimensionless. V s reservoir volume, m3 3 ; p s R is the density of the reservoir rock, kg / m3 3 ; C s Cp for reservoir rock, J / (kg K); T in Tres = initial reservoir temperature, °C; T av (t) is the integrated average reservoir temperature at a certain time, ℃; V f reservoir volume, m3 3 ; Z chem W / m 3 ; T inj For the injected fluid temperature, K.

Images