Numerical simulation method for hydraulic fracture of hot dry rock based on continuous-discontinuous unit method

Abstract

The application discloses a kind of dry hot rock hydraulic fracture numerical simulation methods based on continuous-discontinuous unit method, first, solid stress field calculation model, fluid field calculation model, temperature field calculation model are established, and the fracture criterion of dry hot rock meets maximum tensile stress criterion and Mohr-Coulomb criterion;Then according to the actual situation of dry hot rock, natural fracture development dry hot rock geothermal reservoir hydraulic fracturing model is established;Finally, the established fracturing model is combined with each calculation model described above to solve iteratively, and in the solving calculation process, whether new crack is generated in hydraulic fracturing model is judged using fracture criterion and model crack updating is carried out, until the dry hot rock crack propagation numerical simulation process is completed;The application carries out numerical simulation on deep dry hot rock hydraulic fracture propagation under the premise of considering thermal stress and the influence of natural crack, and the simulation accuracy is higher through test verification, so as to provide data support for subsequent dry hot rock geothermal reservoir fracturing construction parameter optimization.

Description

Technical Field

[0001] This invention belongs to the field of geothermal development technology of hot dry rock, and specifically relates to a numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method. Background Technology

[0002] Medium-deep and deep dry hot rock reservoirs are usually composed of granite, characterized by high temperature, poor matrix permeability, hard texture, and abundant natural fractures. Constructing an enhanced geothermal system (EGS) through hydraulic fracturing is a key technology for the efficient development of geothermal energy.

[0003] Establishing a complex, interconnected network of fractures underground using hydraulic fracturing technology is a major challenge in the development of geothermal reservoir fracturing (EGS) in hot dry rock. Due to significant geological influences on actual experiments, current methods primarily involve building models and conducting numerical simulations of fracture propagation under different geothermal reservoir conditions to provide data support for subsequent construction. However, existing numerical simulation models for fracture propagation in geothermal reservoirs are mostly two-dimensional or small-scale three-dimensional models. Some models fail to couple temperature field calculations, cannot analyze the hydraulic fracture propagation patterns under the influence of thermal stress and natural fractures, and lack accurate and efficient field-scale thermo-fluid-solid coupled three-dimensional fracture propagation numerical models. Artificial fracture propagation is difficult to predict and control, ultimately failing to effectively guide the optimization design of fracturing construction parameters in hot dry rock geothermal reservoirs.

[0004] Therefore, the research direction of this invention is to provide a new method that can numerically simulate the propagation of hydraulic fractures in deep dry hot rock under the premise of considering the influence of thermal stress and natural fractures, so as to provide data support for the optimization of subsequent hydraulic fracturing construction parameters of dry hot rock geothermal reservoirs. Summary of the Invention

[0005] To address the problems existing in the prior art, this invention provides a numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method. This method can numerically simulate the propagation of hydraulic fractures in deep hot dry rock under the premise of considering the influence of thermal stress and natural fractures, thereby providing data support for the optimization of subsequent hydraulic fracturing construction parameters in hot dry rock geothermal reservoirs.

[0006] To achieve the above objectives, the technical solution adopted by this invention is: a numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method, comprising the following steps:

[0007] Step 1: Determine the solid stress field calculation model: Based on d'Alembert's principle, establish the motion control equations of the solid stress field. Then, considering the fluid pressure and thermal stress conditions, combine Hooke's law to determine the solid stress field calculation model, which is used to calculate the solid stress data of dry hot rock.

[0008] Step 2: Determine the fluid field calculation model: Based on the principle of mass conservation, establish the fluid field governing equations. Then, combining Darcy's law and Gaussian divergence theorem, determine the fluid field calculation model, which is used to calculate the porosity data of hot dry rock considering the influence of solid stress field on pore seepage.

[0009] Step 3: Determine the temperature field calculation model: First, establish the temperature field control equations. Then, based on Fourier's law of thermal conductivity and Gaussian divergence theorem, use the finite volume method to solve the heat conduction process of the block element, thereby determining the temperature field calculation model, which is used to calculate the temperature data at different locations of the dry hot rock.

[0010] Step 4: Determine the fracture criteria: Set the fracture of hot dry rock to satisfy the maximum tensile stress criterion and the Mohr-Coulomb criterion, which are used to determine tensile fracture and shear fracture.

[0011] Step 5: Establish a hydraulic fracturing model for dry hot rock geothermal reservoirs with naturally developed fractures: First, establish a hydraulic fracturing fracture propagation model for dry hot rock reservoirs without naturally developed fractures, and then establish a hydraulic fracturing model for dry hot rock geothermal reservoirs with naturally developed fractures based on this model.

[0012] Step Six: Numerical Simulation of Fracture Propagation in Hot Dry Rocks: The hydraulic fracturing model of the hot dry rock geothermal reservoir with natural fractures established in Step Five is used in conjunction with the solid stress field calculation model, fluid field calculation model, and temperature field calculation model determined in Steps One to Three. The model is solved explicitly iteratively using dynamic relaxation technology. During the solution calculation process, the model is updated to determine whether new fractures are generated in the hydraulic fracturing model of the hot dry rock geothermal reservoir based on the fracture criterion determined in Step Four. At the same time, the aperture of the new fractures is calculated until the numerical simulation process of fracture propagation in hot dry rocks is completed.

[0013] Furthermore, step one specifically includes:

[0014] According to d'Alembert's principle, when an object is in a state of dynamic "equilibrium," the governing equations for the motion of the solid stress field are:

[0015]

[0016] In the formula, σ′ ij Let x be the total effective stress tensor, Pa; j Let m be the coordinate; b be the coordinates. i For volume force, N·m -3 ;ρ s Density of rock, kg·m -3 ;u i t is displacement, m; t is time, s; c is damping coefficient, N·s·m -4 ;

[0017] Due to the effects of pore pressure and temperature changes, the constitutive relation of solids needs to consider the effects of fluid pressure and thermal stress. The effective stress of porous media is:

[0018] σ′ ij =σ ij -α B p m δ ij

[0019] Thermal stress caused by temperature changes is calculated by the following formula:

[0020]

[0021] For linear elastic materials, the stress-strain relationship satisfies Hooke's Law:

[0022]

[0023] In the formula, σ ij Let α be the total effective stress tensor, Pa; B p is the Biot coefficient, dimensionless; m δ represents pore pressure, Pa; ij Kronecker notation; σ T α is thermal stress, Pa; E is elastic modulus, Pa; α T The coefficient of thermal expansion is given in °C. -1 T represents the current temperature, in °C. ref Reference temperature (°C); ν is Poisson's ratio; G is shear modulus (Pa), which can be derived from... Calculate; ε ij Let be the strain tensor, which is dimensionless; λ be the Lamé coefficient, Pa, which can be derived from... calculate;

[0024] The constitutive relation of hot dry rock material is:

[0025]

[0026] Strain and displacement satisfy the geometric equation:

[0027]

[0028] In the formula, u i,j and u j,i All are first-order partial derivatives of displacement with respect to coordinates, and have no dimension.

[0029] Furthermore, step two specifically involves:

[0030] Based on the principle of mass conservation, the governing equations of the fluid field are established as follows:

[0031]

[0032] In the formula, S f Pa is the fluid compressibility coefficient. -1 v is the fluid velocity, in m·s -1 ;q ap For source (sink) items, s -1 ;

[0033] For pore fluid flow, Darcy's law applies:

[0034]

[0035] In the formula, v i Let be the node velocity in the i-direction, m·s -1 ;p E The total pressure at the pore nodes is expressed in Pa and k. i Permeability in the i-direction, μm 2 μ is the fluid viscosity, Pa·s; k s The permeability is relative and dimensionless, calculated based on the average saturation of the pore units.

[0036]

[0037] In the formula, N e s represents the total number of nodes in the porous element, dimensionless; Ek Let be the saturation of the k-th node of the porosity element; it is dimensionless.

[0038] According to Gaussian divergence theorem, the nodal seepage velocity can be expressed as:

[0039]

[0040] In the formula, V is the volume of the pore element, m 3 N m S represents the total number of surfaces contained in a pore element, dimensionless; j Let m be the area of ​​the j-th face of the pore element. 2 ; p is the total pressure p at all nodes on the j-th face of the pore element. E The average value, Pa; n ij Let i be the component of the unit normal vector of the j-th face of the pore element in the i-th direction, which is dimensionless;

[0041] The flow rate of the pore element was calculated as follows:

[0042]

[0043] In the formula, N d is the total number of faces associated with this node in the bulk element, dimensionless; v is the seepage velocity vector of the pore element, m / s; n jN is the unit outward normal vector of the j-th face of the pore element, dimensionless; j is the total number of nodes on the j-th face of the porosity element, which is dimensionless;

[0044] Therefore, the total pore flow rate at a node is:

[0045]

[0046] In the formula, N c q represents the total number of elements connected to the node. j For the flow rate of the j-th pore element connected to the node, m 3 ·s -1 The change in node saturation is:

[0047]

[0048] The node saturation is then:

[0049] s E (t1)=s E (t0)+Δs E

[0050] In the formula, Q app For external flow boundary conditions, m 3 ·s -1 φ represents nodal porosity, which is dimensionless; V n Let m be the node volume. 3 ;

[0051] When the node saturation is less than 1, the node fluid pressure is considered to be 0; when the node saturation reaches 1, the change in node fluid pressure is...

[0052]

[0053] The total pressure of the fluid at the node is:

[0054] p E (t1)=p E (t0)+Δp E

[0055] The calculation method for fracture seepage is similar to that for pore seepage. The pore elements in the above calculation process are replaced with fracture elements. When performing fracture seepage calculations, the volumes and surfaces involved in the pore seepage calculation process need to be reduced to surfaces and line segments, respectively. For fracture seepage, its flow law follows the cubic law, so the nodal velocity of the fracture element is:

[0056]

[0057] In the formula, p FThe total fluid pressure at the fracture element node is Pa; w e Let m be the opening of the crack element node, which can be calculated based on the displacement difference between the nodes of the block elements on both sides of the crack element.

[0058] The average pore pressure of a pore element is obtained by averaging the total fluid pressure at all nodes of the pore element. Similarly, the average pressure of the fracture element is obtained by averaging the total fluid pressure at all nodes on the fracture element. The fluid flow rate between the pore and fracture elements can be calculated using Darcy's law:

[0059]

[0060] In the formula, k E Permeability of pore units, in μm 2 ; d is the perpendicular distance between the centroid of the pore element and the surface of the fracture element, in meters; A is the contact area between the pore element and the fracture element, in square meters. 2 ;

[0061] By Q EF These are used as flow boundary conditions for pore flow and fracture flow calculations, respectively, to achieve coupling between pores and fractures. Considering the influence of the solid stress field on pore flow, the porosity of hot dry rock is calculated by the following formula:

[0062] φ=φ0+α B ε v

[0063] In the formula, φ0 is the initial porosity, which is dimensionless; ε v The volumetric strain of the block element is dimensionless.

[0064] Furthermore, step three specifically includes:

[0065] Temperature field calculations are influenced by heat conduction, heat convection, and heat source terms. Based on the assumption of local thermal equilibrium, the governing equations for the temperature field of hot dry rock are obtained as follows:

[0066]

[0067] In the formula, (ρC) eff =ρ f C f (1-φ)+ρ s C s φ, λ eff =λ f (1-φ)+λ s φ;T m Q represents the reservoir rock temperature, in °C. fr As a heat source, W·m -3 ;ρ f For fluid density, kg·m-3 C f Specific heat capacity of the fluid, J·kg -1 ·℃ -1 ;ρ s Density of the rock skeleton, kg·m -3 C s Specific heat capacity of the rock skeleton, J·kg -1 ·℃ -1 v is the fluid velocity vector, in m·s -1 ;λ f The thermal conductivity of the fluid is W·m. -1 ·℃ -1 ;λ s The thermal conductivity of the rock skeleton is W·m. -1 ·℃ -1 ;

[0068] Heat conduction follows Fourier's law of thermal conductivity. Based on Gaussian divergence theorem, the finite volume method is used to solve for the heat conduction process of a bulk element. First, the element heat flow rate is calculated based on the nodal temperatures.

[0069]

[0070] In the formula, q Ti Let be the heat flux velocity of the bulk element in the i-direction, W·m. -2 V represents the volume of a block unit, in meters. 3 N m S is the total number of faces contained in the block unit, dimensionless; j Let m be the area of ​​the j-th face of the block element. 2 ; Let T be the temperature of all nodes on the j-th face of the block element. m The average value, ℃; n ij Let i be the component of the unit normal vector of the j-th face of the block element in the i-th direction, which is dimensionless;

[0071] Calculate the heat flow rate through each block element node based on the element heat flow rate:

[0072]

[0073] In the formula, N d q represents the total number of faces associated with this node in the block element, dimensionless; T Let W·m be the heat flux vector of the bulk element. -2 ;n j N is the unit outward normal vector of the j-th face of the block element, dimensionless; j The total number of nodes on the j-th face of the block element is dimensionless.

[0074] When there are multiple units, the heat flux at the common node needs to be superimposed to obtain the total heat flux generated by heat conduction at a single node:

[0075]

[0076] In the formula, N c Q represents the total number of block elements connected to this node, dimensionless; Tj Let W be the heat flux at the j-th block unit connected to this node;

[0077] The heat source term is the heat flux boundary at each block element node:

[0078] Q app =V n Q fr

[0079] In the formula, V n Let m be the volume of a block element node. 3 ;

[0080] The fluid velocity in the thermal convection term is taken from the fluid seepage velocity calculated in the previous section, and then the heat flow rate generated by thermal convection is derived based on the Gaussian divergence theorem:

[0081]

[0082] In the formula, v x Let be the component of the pore seepage velocity in the x-direction, in m·s. -1 ;v y Let be the component of the pore flow velocity in the y-direction, in m·s. -1 ;v z Let be the component of the pore flow velocity in the z-direction, in m·s. -1 ;n xj Let n be the x-component of the unit normal vector of the j-th face of the block element, dimensionless; yj Let n be the dimensionless component of the unit normal vector of the j-th face of the block element in the y-direction; zj Let z be the component of the unit normal vector of the j-th face of the block element in the z-direction, which is dimensionless;

[0083] Similarly, for multiple units, the convective heat flow at the common node is superimposed to obtain the total heat flow generated by thermal convection at a single node:

[0084]

[0085] In the formula, N n This represents the total number of nodes in the block element, which is dimensionless.

[0086] The temperature change at a node is:

[0087]

[0088] The node temperature after time Δt is: T m (t+Δt)=T m (t)+ΔT m .

[0089] Furthermore, step four specifically includes:

[0090] In the numerical simulation, the fracture behavior of the rock is described by the fracture of the interface contact spring. The relative displacement of the two ends of the contact spring at the interface and the spring force of the contact spring satisfy Hooke's law. Therefore, the normal and tangential test contact forces in the next time step are expressed by the incremental method:

[0091]

[0092] In the formula, F n F s These are the normal and tangential forces on the contact spring, respectively, in N; Δu n , Δu s Let m and K be the normal and tangential relative displacements of the two ends of the contact spring, respectively. n K s These are the normal and tangential stiffnesses, respectively, in N·m. -1 ;

[0093] Calculate the normal and tangential contact forces of the contact spring based on the nodal displacement. When the normal contact force satisfies the maximum tensile stress criterion:

[0094] -F n (t1)≥σ t (t0)A c

[0095] At this point, tensile fracture occurs at the interface unit, and the normal contact force and tensile strength are corrected as follows:

[0096]

[0097] If the tangential contact force satisfies the Mohr-Coulomb criterion:

[0098]

[0099] At this point, shear failure occurs at the interface unit, and the tangential contact force and cohesion are corrected as follows:

[0100]

[0101] In the formula, σ t The tensile strength of the material is expressed in Pa and A. c m is the equivalent area of ​​the node. 2 ; is the internal friction angle of the material, °; c is the cohesive force of the material, Pa; t0 and t1 represent the current moment and the next moment, respectively;

[0102] After fracture occurs, the crack aperture is obtained by the displacement difference between the nodes at both ends of the interface contact spring:

[0103] w = |(u1-u2)·n 12 |

[0104] In the formula, u1 and u2 are the displacements of the connecting nodes at both ends of the spring, m and n, respectively. 12 Let be the unit normal vector of the contact surface, which is dimensionless.

[0105] Compared with existing technologies, this invention first establishes solid stress field calculation models, fluid field calculation models, and temperature field calculation models based on the condition of hot dry rocks, and determines that the fracture criteria of hot dry rocks satisfy the maximum tensile stress criterion and the Mohr-Coulomb criterion; then, based on the actual situation of the hot dry rocks to be simulated, a hydraulic fracturing model of hot dry rock geothermal reservoirs with natural fractures is established; finally, the established fracturing model is combined with the above-mentioned solid stress field calculation models, fluid field calculation models, and temperature field calculation models, and a dynamic relaxation technique is used for explicit iterative solution, and fracture criterion is used to determine the hydraulic fracturing of hot dry rock geothermal reservoirs during the solution calculation process. The model is updated to reflect new cracks and the opening of new cracks is calculated until the numerical simulation of crack propagation in hot dry rock is completed. This invention establishes a three-dimensional thermal-fluid-structure interaction crack propagation mathematical model based on the continuous-discontinuous element method. Under the premise of considering the influence of thermal stress and natural cracks, it performs numerical simulation of hydraulic crack propagation in deep hot dry rock. The simulation accuracy has been verified by experiments to be high. It can analyze the influence of thermal stress and natural cracks on hydraulic crack propagation during the fracturing process, thereby providing data support for the optimization of subsequent fracturing construction parameters in hot dry rock geothermal reservoirs. Attached Figure Description

[0106] Figure 1 This is a flowchart of a numerical simulation of crack propagation in hot dry rock according to an embodiment of the present invention;

[0107] Figure 2 This is a hydraulic fracturing fracture propagation model for a dry hot rock reservoir with underdeveloped natural fractures, established in the embodiments of the present invention.

[0108] Figure 3 This is a hydraulic fracturing fracture propagation model for a naturally fractured hot dry rock reservoir established in this embodiment of the invention;

[0109] Figure 4 It is the Penny three-dimensional hydraulic fracturing numerical model established in the verification experiment;

[0110] Figure 5This is a numerical simulation diagram of the crack width and pressure distribution in the Penny crack during the verification experiment;

[0111] Figure 6 This is a comparison chart of the numerical simulation and analytical solution results of hydraulic fracturing fracture propagation in Case A of the verification experiment;

[0112] Figure 7 This is a comparison chart of the numerical simulation and analytical solution results of hydraulic fracturing fracture propagation in Case B of the verification experiment;

[0113] Figure 8 This is a cloud map of crack width distribution under different grid sizes in the verification experiment;

[0114] Figure 9 It is a graph showing the evolution of slit length and slit width for different grid sizes in the verification experiment;

[0115] Figure 10 It is a numerical model for the intersection of hydraulic cracks and natural cracks in the verification experiment;

[0116] Figure 11 This is a comparison chart of numerical simulation and indoor experimental results of the intersection behavior of hydraulic cracks and natural cracks in the verification experiment. Detailed Implementation

[0117] The present invention will be further described below.

[0118] like Figure 1 As shown, the present invention includes the following steps:

[0119] Step 1: Determine the solid stress field calculation model: Based on d'Alembert's principle, establish the governing equations of motion for the solid stress field. Then, considering the conditions of fluid pressure and thermal stress, and combining Hooke's law, determine the solid stress field calculation model for calculating the solid stress data of hot dry rock. Specifically:

[0120] According to d'Alembert's principle, when an object is in a state of dynamic "equilibrium," the governing equations for the motion of the solid stress field are:

[0121]

[0122] In the formula, σ′ ij Let x be the total effective stress tensor, Pa; j Let m be the coordinate; b be the coordinates. i For volume force, N·m -3 ;ρ s Density of rock, kg·m -3 ;u i t is displacement, m; t is time, s; c is damping coefficient, N·s·m -4 ;

[0123] Due to the effects of pore pressure and temperature changes, the constitutive relation of solids needs to consider the effects of fluid pressure and thermal stress. The effective stress of porous media is:

[0124] σ′ ij =σ ij -α B p m δ ij

[0125] Thermal stress caused by temperature changes is calculated by the following formula:

[0126]

[0127] For linear elastic materials, the stress-strain relationship satisfies Hooke's Law:

[0128]

[0129] In the formula, σ ij Let α be the total effective stress tensor, Pa; B p is the Biot coefficient, dimensionless; m δ represents pore pressure, Pa; ij Kronecker notation; σ T α is thermal stress, Pa; E is elastic modulus, Pa; α T The coefficient of thermal expansion is given in °C. -1 T represents the current temperature, in °C. ref Reference temperature (°C); ν is Poisson's ratio; G is shear modulus (Pa), which can be derived from... Calculate; ε ij Let be the strain tensor, which is dimensionless; λ be the Lamé coefficient, Pa, which can be derived from... calculate;

[0130] The constitutive relation of hot dry rock material is:

[0131]

[0132] Strain and displacement satisfy the geometric equation:

[0133]

[0134] In the formula, u i,j and u j,i All are first-order partial derivatives of displacement with respect to coordinates, and have no dimension.

[0135] Step 2: Determine the fluid field calculation model: Based on the principle of mass conservation, establish the fluid field governing equations. Then, combining Darcy's law and Gaussian divergence theorem, determine the fluid field calculation model. This model is used to calculate the porosity data of hot dry rock considering the influence of solid stress field on pore seepage. Specifically:

[0136] Based on the principle of mass conservation, the governing equations of the fluid field are established as follows:

[0137]

[0138] In the formula, S f Pa is the fluid compressibility coefficient. -1 v is the fluid velocity, in m·s -1 ;q ap For source (sink) items, s -1 ;

[0139] For pore fluid flow, Darcy's law applies:

[0140]

[0141] In the formula, v i Let be the node velocity in the i-direction, m·s -1 ;p E The total pressure at the pore nodes is expressed in Pa and k. i Permeability in the i-direction, μm 2 μ is the fluid viscosity, Pa·s; k s The permeability is relative and dimensionless, calculated based on the average saturation of the pore units.

[0142]

[0143] In the formula, N e s represents the total number of nodes in the porous element, dimensionless; Ek Let be the saturation of the k-th node of the porosity element; it is dimensionless.

[0144] According to Gaussian divergence theorem, the nodal seepage velocity can be expressed as:

[0145]

[0146] In the formula, V is the volume of the pore element, m 3 N m S represents the total number of surfaces contained in a pore element, dimensionless; j Let m be the area of ​​the j-th face of the pore element. 2 ; p is the total pressure p at all nodes on the j-th face of the pore element. E The average value, Pa; n ij Let i be the component of the unit normal vector of the j-th face of the pore element in the i-th direction, which is dimensionless;

[0147] The flow rate of the pore element was calculated as follows:

[0148]

[0149] In the formula, N d is the total number of faces associated with this node in the bulk element, dimensionless; v is the seepage velocity vector of the pore element, m / s; n j N is the unit outward normal vector of the j-th face of the pore element, dimensionless; j is the total number of nodes on the j-th face of the porosity element, which is dimensionless;

[0150] Therefore, the total pore flow rate at a node is:

[0151]

[0152] In the formula, N c q represents the total number of elements connected to the node. j For the flow rate of the j-th pore element connected to the node, m 3 ·s -1 ;

[0153] The change in node saturation is:

[0154]

[0155] The node saturation is then:

[0156] s E (t1)=s E (t0)+Δs E

[0157] In the formula, Q app For external flow boundary conditions, m 3 ·s -1 φ represents nodal porosity, which is dimensionless; V n Let m be the node volume. 3 ;

[0158] When the node saturation is less than 1, the node fluid pressure is considered to be 0; when the node saturation reaches 1, the change in node fluid pressure is...

[0159]

[0160] The total pressure of the fluid at the node is:

[0161] p E (t1)=p E (t0)+Δp E

[0162] The calculation method for fracture seepage is similar to that for pore seepage. The pore elements in the above calculation process are replaced with fracture elements. When performing fracture seepage calculations, the volumes and surfaces involved in the pore seepage calculation process need to be reduced to surfaces and line segments, respectively. For fracture seepage, its flow law follows the cubic law, so the nodal velocity of the fracture element is:

[0163]

[0164] In the formula, p F The total fluid pressure at the fracture element node is Pa; w e Let m be the opening of the crack element node, which can be calculated based on the displacement difference between the nodes of the block elements on both sides of the crack element.

[0165] The average pore pressure of a pore element is obtained by averaging the total fluid pressure at all nodes of the pore element. Similarly, the average pressure of the fracture element is obtained by averaging the total fluid pressure at all nodes on the fracture element. The fluid flow rate between the pore and fracture elements can be calculated using Darcy's law:

[0166]

[0167] In the formula, k E Permeability of pore units, in μm 2 ; d is the perpendicular distance between the centroid of the pore element and the surface of the fracture element, in meters; A is the contact area between the pore element and the fracture element, in square meters. 2 ;

[0168] By Q EF These are used as flow boundary conditions for pore flow and fracture flow calculations, respectively, to achieve coupling between pores and fractures. Considering the influence of the solid stress field on pore flow, the porosity of hot dry rock is calculated by the following formula:

[0169] φ=φ0+α B ε v

[0170] In the formula, φ0 is the initial porosity, which is dimensionless; ε v The volumetric strain of the block element is dimensionless.

[0171] Step 3: Determine the temperature field calculation model: First, establish the temperature field governing equations. Then, based on Fourier's law of thermal conductivity and Gaussian divergence theorem, use the finite volume method to solve the heat conduction process of the bulk element, thereby determining the temperature field calculation model used to calculate temperature data at different locations in the hot dry rock. Specifically:

[0172] Temperature field calculations are influenced by heat conduction, heat convection, and heat source terms. Based on the assumption of local thermal equilibrium, the governing equations for the temperature field of hot dry rock are obtained as follows:

[0173]

[0174] In the formula, (ρC) eff =ρ f C f (1-φ)+ρ s C s φ, λ eff =λ f (1-φ)+λ s φ;T m Q represents the reservoir rock temperature, in °C. fr As a heat source, W·m -3 ;ρ f For fluid density, kg·m -3 C f Specific heat capacity of the fluid, J·kg -1 ·℃ -1 ;ρ s Density of the rock skeleton, kg·m -3 C s Specific heat capacity of the rock skeleton, J·kg -1 ·℃ -1 v is the fluid velocity vector, in m·s -1 ;λ f The thermal conductivity of the fluid is W·m. -1 ·℃ -1 ;λ s The thermal conductivity of the rock skeleton is W·m. -1 ·℃ -1 ;

[0175] Heat conduction follows Fourier's law of thermal conductivity. Based on Gaussian divergence theorem, the finite volume method is used to solve for the heat conduction process in a bulk element. First, the element heat flow rate is calculated based on the nodal temperatures.

[0176]

[0177] In the formula, q Ti Let be the heat flux velocity of the bulk element in the i-direction, W·m. -2 V represents the volume of a block unit, in meters. 3 N m S is the total number of faces contained in the block unit, dimensionless; j Let m be the area of ​​the j-th face of the block element. 2 ; Let T be the temperature of all nodes on the j-th face of the block element. m The average value, ℃; n ij Let i be the component of the unit normal vector of the j-th face of the block element in the i-th direction, which is dimensionless;

[0178] Calculate the heat flow rate through each block element node based on the element heat flow rate:

[0179]

[0180] In the formula, N d q represents the total number of faces associated with this node in the block element, dimensionless; T Let W·m be the heat flux vector of the bulk element. -2 ;n j N is the unit outward normal vector of the j-th face of the block element, dimensionless; j The total number of nodes on the j-th face of the block element is dimensionless.

[0181] When there are multiple units, the heat flux at the common node needs to be superimposed to obtain the total heat flux generated by heat conduction at a single node:

[0182]

[0183] In the formula, N c Q represents the total number of block elements connected to this node, dimensionless; Tj Let W be the heat flux at the j-th block unit connected to this node;

[0184] The heat source term is the heat flux boundary at each block element node:

[0185] Q app =V n Q fr

[0186] In the formula, V n Let m be the volume of a block element node. 3 ;

[0187] The fluid velocity in the thermal convection term is taken from the fluid seepage velocity calculated in the previous section, and then the heat flow rate generated by thermal convection is derived based on the Gaussian divergence theorem:

[0188]

[0189] In the formula, v x Let be the component of the pore seepage velocity in the x-direction, in m·s. -1 ;v y Let be the component of the pore flow velocity in the y-direction, in m·s. -1 ;v z Let be the component of the pore flow velocity in the z-direction, in m·s. -1 ;n xj Let n be the x-component of the unit normal vector of the j-th face of the block element, dimensionless; yj Let n be the dimensionless component of the unit normal vector of the j-th face of the block element in the y-direction;zj Let z be the component of the unit normal vector of the j-th face of the block element in the z-direction, which is dimensionless;

[0190] Similarly, for multiple units, the convective heat flow at the common node is superimposed to obtain the total heat flow generated by thermal convection at a single node:

[0191]

[0192] In the formula, N n This represents the total number of nodes in the block element, which is dimensionless.

[0193] The temperature change at a node is:

[0194]

[0195] The node temperature after time Δt is: T m (t+Δt)=T m (t)+ΔT m .

[0196] Step 4: Determine the fracture criteria: The fracture of hot, dry rocks is assumed to satisfy the maximum tensile stress criterion and the Mohr-Coulomb criterion, used to determine tensile and shear fractures, specifically:

[0197] In the numerical simulation process of CDEM-THM3D software, the fracture behavior of the rock is described by the fracture of the interface contact spring. The relative displacement of the two ends of the contact spring at the interface and the spring force of the contact spring satisfy Hooke's law. Therefore, the normal and tangential test contact forces in the next time step are expressed by the incremental method:

[0198]

[0199] In the formula, F n F s These are the normal and tangential forces on the contact spring, respectively, in N; Δu n , Δu s Let m and K be the normal and tangential relative displacements of the two ends of the contact spring, respectively. n K s These are the normal and tangential stiffnesses, respectively, in N·m. -1 ;

[0200] Calculate the normal and tangential contact forces of the contact spring based on the nodal displacement. When the normal contact force satisfies the maximum tensile stress criterion:

[0201] -F n (t1)≥σ t (t0)A c

[0202] At this point, tensile fracture occurs at the interface unit, and the normal contact force and tensile strength are corrected as follows:

[0203]

[0204] If the tangential contact force satisfies the Mohr-Coulomb criterion:

[0205]

[0206] At this point, shear failure occurs at the interface unit, and the tangential contact force and cohesion are corrected as follows:

[0207]

[0208] In the formula, σ t The tensile strength of the material is expressed in Pa and A. c m is the equivalent area of ​​the node. 2 ; is the internal friction angle of the material, °; c is the cohesive force of the material, Pa; t0 and t1 represent the current moment and the next moment, respectively;

[0209] After fracture occurs, the crack aperture is obtained by the displacement difference between the nodes at both ends of the interface contact spring:

[0210] w = |(u1-u2)·n 12 |

[0211] In the formula, u1 and u2 are the displacements of the connecting nodes at both ends of the spring, m and n, respectively. 12 Let be the unit normal vector of the contact surface, which is dimensionless.

[0212] Step 5: Establish a hydraulic fracturing model for a geothermal reservoir with naturally developed dry fractures: First, establish a hydraulic fracturing fracture propagation model for a dry geothermal reservoir without naturally developed fractures, and then establish a hydraulic fracturing model for a geothermal reservoir with naturally developed fractures based on this model. In this embodiment, the basic parameters of the numerical model are set by testing the physical and mechanical parameters of granite and conventional numerical settings for medium-deep (approximately 2000m) granite dry geothermal reservoirs, as shown in Table 1. The model achieves the setting of triaxial geostress by applying surface forces to the upper, left, and front surfaces, and setting normal displacement constraints on the lower, right, and rear surfaces, and sets a high-stress barrier. The model is divided into hexahedral meshes. To maintain computational accuracy and reduce the number of elements to improve computational efficiency, 1m meshes are set in the y and z directions (fracture length and fracture height directions), and 1m meshes are set in the middle and 8m meshes are set on both sides in the x direction (fracture width direction), resulting in 43,200 block elements and 133,980 interface elements (fracture elements). Figure 2 As shown;

[0213] Table 1

[0214]

[0215]

[0216] Under the influence of tectonic stress and local faults, deep hot dry rock geothermal reservoirs can develop abundant natural fractures and joints. In numerical simulations, the strength of natural fractures is usually weakened while their permeability is enhanced to distinguish between hot dry rock matrix units and natural fracture units. When modeling with CDEM-THM3D, interface units (fracture units) need to be divided along the natural fractures, and the contact stiffness, tensile strength, cohesion, and internal friction angle of the interface units at the natural fractures are reduced, while the initial aperture (i.e., permeability) of the fracture units at the natural fractures is increased. Based on the above-mentioned fracturing model of hot dry rock reservoirs without well-developed natural fractures, randomly distributed natural fractures are added to form a hydraulic fracturing model of hot dry rock geothermal reservoirs with well-developed natural fractures, such as... Figure 3 As shown in Table 2, the basic parameters of the model are set. Two hundred natural fractures are randomly generated, with their orientations uniformly distributed from 0 to 180°, and their horizontal cross-sectional lengths following a normal distribution with a mean of 20 m and a variance of 4. The tensile strength, cohesion, and internal friction coefficient of the natural fractures are set to be one-tenth that of the matrix, and the permeability of the natural fractures is set to be 1000 times that of the matrix. Maximum horizontal stress is applied in the x-direction, minimum horizontal stress in the y-direction, and vertical stress in the z-direction. The injection point is located at the center of the model, and fracturing fluid is continuously injected at a constant flow rate. This model mainly studies the process of hydraulic fracture communication activating natural fractures to form a complex fracture network in hot dry rock reservoirs, focusing on the simulation results of hydraulic fractures in the fracture length and width directions. Therefore, in order to reduce the amount of computation and improve the computation speed, the model did not set up a partition and was divided into a coarser grid (10m) in the vertical direction. At the same time, in order to ensure the simulation accuracy of crack length and crack width, a finer grid (1m) was divided in the horizontal direction, resulting in 48,775 triangular prism elements and 111,480 interface elements (crack elements).

[0217] Table 2

[0218]

[0219] Step Six: Numerical Simulation of Fracture Propagation in Hot Dry Rocks: The hydraulic fracturing model of the naturally fractured hot dry rock geothermal reservoir established in Step Five is used in conjunction with the solid stress field calculation model, fluid field calculation model, and temperature field calculation model determined in Steps One to Three. The model is solved explicitly iteratively using dynamic relaxation techniques. Figure 1As shown, by introducing a damping term in the dynamic calculation, the initially unbalanced vibration system gradually decays to the equilibrium position. The calculations of the solid stress field and fluid field are implemented using proprietary modules of the GDEM platform (a series of mechanical analysis software jointly developed by the Joint Laboratory of Discontinuous Media Mechanics and Engineering Disasters of the Chinese Academy of Sciences and Beijing Jidao Chengran Technology Co., Ltd.). The temperature field calculation considering convective heat transfer is implemented through C++ secondary development based on the reserved interface of the GDEM platform. Furthermore, CPU parallel technology is used in the solution process of each physical field to improve the calculation speed. And in the solution calculation process, the calculation is determined according to step four. The fracture criterion is used to determine whether a new fracture is generated in the hydraulic fracturing model of a hot dry rock geothermal reservoir and to update the model fracture. At the same time, the aperture of the new fracture is calculated. The nodal displacement and element strain obtained from the solid stress field calculation will change the element porosity and fracture element width, thus affecting the fluid flow calculation in the hot dry rock reservoir. The thermal convection term in the temperature field calculation depends on the fluid flow velocity in the fluid field calculation model. Fluid pressure and thermal stress caused by temperature changes are both used as external forces to affect the calculation of the solid stress field, thereby realizing the coupling of heat-fluid-solid stress and finally completing the numerical simulation process of fracture propagation in hot dry rock.

[0220] Verification experiment:

[0221] Verification Experiment 1:

[0222] Theoretically, the propagation process of three-dimensional hydraulic fracturing fractures conforms to a circular fracture model (i.e., the penny fracture model). Some scholars in the industry have derived theoretical solutions for circular fractures and further provided analytical solutions that approximate complex functions. Based on the dimensionless time τ and the filtration parameter φ, the analytical solution of the circular fracture model is divided into several regions, and analytical solution formulas for the ductile region, viscous region, filtration-free ductile region, and filtration-free viscous region are given. The calculation formulas for the dimensionless time τ and the filtration parameter φ are as follows:

[0223]

[0224] in,

[0225] μ′=12μ, C′=2C L

[0226] In the formula, E is the elastic modulus, Pa; ν is Poisson's ratio; μ is the fluid viscosity, Pa·s; K Ic For fracture toughness, Pa·m 0.5 There is a relationship between fracture toughness and tensile strength. Ic =σ t d 1 / 2 / a, where d is the average size of the rock block unit (m), and a is a numerical parameter (taken as 1); C L Let be the Caterpillar filtration coefficient, m·s-0.5 t is the fracturing time, in seconds; Q t For injection displacement, m 3 ·s -1 .

[0227] The analytical solution formula for a circular crack in the ductile zone is:

[0228]

[0229] In the formula, ρ is the normalized distance, which is the ratio of distance to half-length of the crack; w is the crack width, m; R is the half-length of the crack, m; and p is the pressure inside the crack, Pa.

[0230] The analytical solution formula for a circular crack in the viscous region is:

[0231]

[0232] in,

[0233]

[0234] In the formula, KE1(*) and KE2(*) are the first and second elliptic integrals, respectively.

[0235] Two typical cases of circular fracture propagation in ductile and viscous zones were set up. The numerical simulation results were compared with the analytical solutions to verify the accuracy of the hydraulic fracturing model of naturally fractured dry hot rock geothermal reservoirs established in this invention and the simulation of three-dimensional fracture propagation using the method of this invention. The basic parameter settings for the two cases are shown in Table 3.

[0236] Table 3

[0237]

[0238] like Figure 4 As shown, a three-dimensional hydraulic fracturing numerical model was established. The model size is 30m × 30m × 20m, consisting of two blocks and one contact surface. The blocks are divided into tetrahedral elements, and the contact surface is divided into triangular elements. The mesh size at the contact surface is 0.5m, and the mesh size at the top and bottom faces is 5m, resulting in 52,111 tetrahedral elements and 9,324 contact surface elements (fracture elements). An injection point was set at the center of the model to achieve a constant displacement of 1.5m³ / h. 3 ·min -1 Numerical simulations of crack propagation dominated by toughness and viscosity were performed with continuous injection.

[0239] Figure 5The numerical simulations of crack width and pressure distribution in Case A and Case B show that the cracks are circular, with the maximum width at the injection point and gradually decreasing outwards. In Case A (ductile zone), the pressure within the crack is higher at the injection point and almost the same at other locations, and the pressure spread is greater than the crack length, consistent with the pressure distribution pattern of a circular crack dominated by ductility. In Case B (viscous zone), the pressure within the crack shows a clear decreasing trend from the injection point outwards, and the pressure spread is less than the crack length, consistent with the pressure distribution pattern of a circular crack dominated by viscosity.

[0240] Depend on Figure 6 It can be seen that the numerical simulation results of the ductile zone agree well with the analytical solution. The average error between the numerical and analytical solutions for the fracture length is 2.1%, and the average error for the fracture width versus time curve is 0.5%. The fracture width versus radial distance curve agrees well with the theoretical solution, and the pressure distribution within the fracture conforms to the theoretical law. Figure 7 It can be seen that the combined error between the numerical simulation of the viscous region and the analytical solution is 2.8%, and the combined error of the curve of the slit width changing with time is 0.3%. Moreover, the slit width and the pressure distribution inside the slit are in good agreement with the analytical solution.

[0241] The accuracy of the proposed method in simulating three-dimensional hydraulic fracture propagation problems can be verified using ductility-dominated and viscosity-dominated circular fracture propagation models. To clarify the impact of CDEM-THM3D mesh generation on the numerical simulation results of fracture propagation, comparative simulation cases with contact surface mesh sizes of 0.2m, 0.5m, 0.8m, and 1.0m were set up. The results are as follows: Figure 8 and Figure 9 As shown.

[0242] like Figure 8 As shown in the diagram, the crack width distribution cloud map reveals that the simulated crack morphology under different mesh sizes largely conforms to the characteristics of circular crack propagation. When the mesh size is large, the crack width distribution cloud map appears coarser, while when the mesh size is small, the crack width distribution cloud map is finer. Figure 9As shown in the curves illustrating the changes in crack length and width, it can be seen that as the mesh size increases, the crack length obtained from the numerical simulation decreases while the crack width increases, consistent with previous numerical simulation results. Furthermore, the numerical and analytical solutions are quite consistent across different mesh sizes; the smaller the mesh size, the smaller the error between the numerical and analytical solutions. With a mesh size of 0.2m, the average errors between the numerical simulation and analytical solutions for crack length and width are 4.7% and 2.8%, respectively; with a mesh size of 0.5m, the average errors are 6.1% and 3.7%, respectively; with a mesh size of 0.8m, the average errors are 8.7% and 5.3%, respectively; and with a mesh size of 1.0m, the average errors are 10.8% and 7.5%, respectively. Even with a relatively coarse mesh, computational accuracy meeting engineering requirements can be achieved.

[0243] Verification Experiment 2:

[0244] To ensure the accuracy of this invention's simulation of the propagation behavior of hydraulic fractures encountering natural fractures, verification was performed using hydraulic fracturing model experiments on rock samples with naturally developed fractures, as conducted by previous researchers. Table 4 shows the conditions and results of hydraulic fracture-natural fracture intersection model experiments conducted by Zhou and Gu, respectively, where the stress difference coefficient is equal to (σ... H -σ h ) / σ h This invention establishes a numerical model of the intersection of hydraulic fractures and natural fractures based on CDEM-THM3D, such as... Figure 10 As shown. Zhou and Gu did not apply temperature conditions during the fracturing model experiments; therefore, temperature field calculations were disabled when using CDEM-THM3D to simulate fracture propagation. Two natural fractures of identical size and orientation were placed at equidistant points on either side of the injection point of the model. The angle between the orientation of the natural fractures and the x-direction was used as the approximation angle. Maintaining the approximation angle and the magnitude of the geostress consistent with the experimental conditions, numerical simulations of the intersection of hydraulic fractures and natural fractures were conducted, and the numerical simulation results were compared with the experimental results.

[0245] Table 4

[0246]

[0247]

[0248] In the table above, the intersection behavior of hydraulic fractures and natural fractures can be categorized into three types: penetration, opening, and termination. Penetration refers to the hydraulic fracture directly penetrating the natural fracture without activating it; opening refers to the natural fracture being activated regardless of whether the hydraulic fracture penetrated it; termination refers to the hydraulic fracture ceasing to propagate in front of the natural fracture without activating it.

[0249] Comparison of numerical simulation and laboratory experiment results, for example Figure 11 As shown, the numerical simulation results agree well with the experimental results. In the given case, only when the approximation angle is 30° and the stress difference coefficient is 3.3 does the numerical simulation result show "on" while the experimental result shows "off"; otherwise, the results are identical. According to the simulation results, as the approximation angle and stress difference coefficient increase, hydraulic fractures are more likely to penetrate natural fractures, which is consistent with relevant theoretical understanding.

[0250] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method, characterized in that, Includes the following steps: Step 1: Determine the solid stress field calculation model: Based on d'Alembert's principle, establish the motion control equations of the solid stress field. Then, considering the fluid pressure and thermal stress conditions, combine Hooke's law to determine the solid stress field calculation model, which is used to calculate the solid stress data of dry hot rock. Step 2: Determine the fluid field calculation model: Based on the principle of mass conservation, establish the fluid field governing equations. Then, combining Darcy's law and Gaussian divergence theorem, determine the fluid field calculation model, which is used to calculate the porosity data of hot dry rock considering the influence of solid stress field on pore seepage. Step 3: Determine the temperature field calculation model: First, establish the temperature field control equations. Then, based on Fourier's law of thermal conductivity and Gaussian divergence theorem, use the finite volume method to solve the heat conduction process of the block element, thereby determining the temperature field calculation model, which is used to calculate the temperature data at different locations of the dry hot rock. Step 4: Determine the fracture criteria: Set the fracture of hot dry rock to satisfy the maximum tensile stress criterion and the Mohr-Coulomb criterion, which are used to determine tensile fracture and shear fracture. Step 5: Establish a hydraulic fracturing model for dry hot rock geothermal reservoirs with naturally developed fractures: First, establish a hydraulic fracturing fracture propagation model for dry hot rock reservoirs without naturally developed fractures, and then establish a hydraulic fracturing model for dry hot rock geothermal reservoirs with naturally developed fractures based on this model. Step Six: Numerical Simulation of Fracture Propagation in Hot Dry Rocks: The hydraulic fracturing model of the hot dry rock geothermal reservoir with natural fractures established in Step Five is used in conjunction with the solid stress field calculation model, fluid field calculation model, and temperature field calculation model determined in Steps One to Three. The model is solved explicitly iteratively using dynamic relaxation technology. During the solution calculation process, the model is updated to determine whether new fractures are generated in the hydraulic fracturing model of the hot dry rock geothermal reservoir based on the fracture criterion determined in Step Four. At the same time, the aperture of the new fractures is calculated until the numerical simulation process of fracture propagation in hot dry rocks is completed.

2. The numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method according to claim 1, characterized in that, Step one specifically involves: According to d'Alembert's principle, the governing equations for the motion of the stress field in a solid are: In the formula, σ′ ij Let x be the total effective stress tensor, Pa; j Let the coordinates be m; b i For volume force, N·m -3 ; ρ s Density of rock, kg·m -3 ;u i t is displacement, m; t is time, s; c is damping coefficient, N·s·m -4 ; The constitutive relation of hot dry rock material is: Strain and displacement satisfy the geometric equation: In the formula, u i,j and u j,i All are first-order partial derivatives of displacement with respect to coordinates, and have no dimension.

3. The numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method according to claim 1, characterized in that, Step two specifically involves: Based on the principle of mass conservation, the governing equations of the fluid field are established as follows: In the formula, S f Pa is the fluid compressibility coefficient. -1 v is the fluid velocity, in m·s -1 ;q ap For the source term, s -1 ; For pore fluid flow, Darcy's law applies: In the formula, v i Let be the node velocity in the i-direction, m·s -1 ;p E The total pressure at the pore nodes is expressed in Pa and k. i Permeability in the i-direction, μm 2 μ is the fluid viscosity, Pa·s; k s The permeability is relative and dimensionless, calculated based on the average saturation of the pore units. In the formula, N e s represents the total number of nodes in the porous element, dimensionless; Ek Let be the saturation of the k-th node of the porosity element; it is dimensionless. According to Gaussian divergence theorem, the nodal seepage velocity can be expressed as: In the formula, V is the volume of the pore element, m 3 N m S is the total number of surfaces contained in the pore element, dimensionless. j Let m be the area of ​​the j-th face of the pore element. 2 ; p is the total pressure p at all nodes on the j-th face of the pore element. E The average value, Pa; n ij Let i be the component of the unit normal vector of the j-th face of the pore element in the i-th direction, which is dimensionless; The flow rate of the pore element was calculated as follows: In the formula, N d is the total number of faces associated with this node in the bulk element, dimensionless; v is the seepage velocity vector of the pore element, m / s; n j N is the unit outward normal vector of the j-th face of the pore element, dimensionless; j is the total number of nodes on the j-th face of the porosity element, which is dimensionless; The final calculated total fluid pressure at the node is: p E (t1)=p E (t0)+Δp E The calculation method for fracture seepage is similar to that for pore seepage, except that the pore elements in the above calculation process are replaced with fracture elements. For fracture seepage, its flow law follows the cubic law, so the nodal velocity of the fracture element is: In the formula, p F The total fluid pressure at the fracture element node is Pa; w e Let m be the aperture of the crack element node; The average pore pressure of a pore element is obtained by averaging the total fluid pressure at all nodes of the pore element. Similarly, the average pressure of the fracture element is obtained by averaging the total fluid pressure at all nodes on the fracture element. The fluid flow rate between the pore and fracture elements can be calculated using Darcy's law: In the formula, k E Permeability of pore units, in μm 2 ; d is the perpendicular distance between the centroid of the pore element and the surface of the fracture element, in meters; A is the contact area between the pore element and the fracture element, in square meters. 2 ; By Q EF These are used as flow boundary conditions for pore flow and fracture flow calculations, respectively, to achieve coupling between pores and fractures. Considering the influence of the solid stress field on pore flow, the porosity of hot dry rock is calculated by the following formula: φ=φ0+α B e v In the formula, φ0 is the initial porosity, which is dimensionless; ε v The volumetric strain of the block element is dimensionless.

4. The numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method according to claim 1, characterized in that, Step three specifically involves: Temperature field calculations are influenced by heat conduction, heat convection, and heat source terms. Based on the assumption of local thermal equilibrium, the governing equations for the temperature field of hot dry rock are obtained as follows: In the formula, (ρC) eff =ρ f C f (1-φ)+ρ s C s φ, λ eff =λ f (1-φ)+λ s φ;T m Q represents the reservoir rock temperature, in °C. fr As a heat source, W·m -3 ;ρ f For fluid density, kg·m -3 C f Specific heat capacity of the fluid, J·kg -1 ·℃ -1 ;ρ s Density of the rock skeleton, kg·m -3 C s Specific heat capacity of the rock skeleton, J·kg -1 ·℃ -1 v is the fluid velocity vector, in m·s -1 ; λ f The thermal conductivity of the fluid is W·m. -1 ·℃ -1 ; λ s The thermal conductivity of the rock skeleton is W·m. -1 ·℃ -1 ; Heat conduction follows Fourier's law of thermal conductivity. Based on Gaussian divergence theorem, the finite volume method is used to solve for the heat conduction process of a bulk element. First, the element heat flow rate is calculated based on the nodal temperatures. In the formula, q Ti Let be the heat flux velocity of the bulk element in the i-direction, W·m. -2 V represents the volume of a block unit, in meters. 3 N m S is the total number of faces contained in the block unit, dimensionless; j Let m be the area of ​​the j-th face of the block element. 2 ; Let T be the temperature of all nodes on the j-th face of the block element. m The average value, ℃; n ij Let i be the component of the unit normal vector of the j-th face of the block element in the i-th direction, which is dimensionless; Calculate the heat flow rate through each block element node based on the element heat flow rate: In the formula, N d q represents the total number of faces associated with this node in the block element, dimensionless; T Let W·m be the heat flux vector of the bulk element. -2 ;n j N is the unit outward normal vector of the j-th face of the block element, dimensionless; j The total number of nodes on the j-th face of the block element is dimensionless. When there are multiple units, the heat flux at the common node needs to be superimposed to obtain the total heat flux generated by heat conduction at a single node: In the formula, N c Q represents the total number of block elements connected to this node, dimensionless; Tj Let W be the heat flux at the j-th block unit connected to this node; The heat source term is the heat flux boundary at each block element node: Q app =V n Q fr In the formula, V n Let m be the volume of a block element node. 3； The fluid velocity in the thermal convection term is taken from the fluid seepage velocity calculated in the previous section, and then the heat flow rate generated by thermal convection is derived based on the Gaussian divergence theorem: In the formula, v x Let be the component of the pore seepage velocity in the x-direction, in m·s. -1 ;v y Let be the component of the pore flow velocity in the y-direction, in m·s. -1 ;v z Let be the component of the pore flow velocity in the z-direction, in m·s. -1 ;n xj Let n be the x-component of the unit normal vector of the j-th face of the block element, dimensionless; yj Let n be the dimensionless component of the unit normal vector of the j-th face of the block element in the y-direction; zj Let z be the component of the unit normal vector of the j-th face of the block element in the z-direction, which is dimensionless; Similarly, for multiple units, the convective heat flow at the common node is superimposed to obtain the total heat flow generated by thermal convection at a single node: In the formula, N n This represents the total number of nodes in the block element, which is dimensionless. The temperature change at a node is: The node temperature after time Δt is: T m (t+Δt)=T m (t)+ΔT m .

5. The numerical simulation method for hydraulic fractures in hot dry rock based on the continuous-discontinuous element method according to claim 1, characterized in that, Step four specifically involves: In the numerical simulation, the fracture behavior of the rock is described by the fracture of the interface contact spring. The relative displacement of the two ends of the contact spring at the interface and the spring force of the contact spring satisfy Hooke's law. Therefore, the normal and tangential test contact forces in the next time step are expressed by the incremental method: In the formula, F n F s These are the normal and tangential forces on the contact spring, respectively, in N; Δu n , Δu s Let m and K be the normal and tangential relative displacements of the two ends of the contact spring, respectively. n K s These are the normal and tangential stiffnesses, respectively, in N·m. -1 ; Calculate the normal and tangential contact forces of the contact spring based on the nodal displacement. When the normal contact force satisfies the maximum tensile stress criterion: -F n (t1)≥σ t (t0)A c At this point, tensile fracture occurs at the interface unit, and the normal contact force and tensile strength are corrected as follows: If the tangential contact force satisfies the Mohr-Coulomb criterion: At this point, shear failure occurs at the interface unit, and the tangential contact force and cohesion are corrected as follows: In the formula, σ t The tensile strength of the material is expressed in Pa and A. c m is the equivalent area of ​​the node. 2 ; is the internal friction angle of the material, °; c is the cohesive force of the material, Pa; t0 and t1 represent the current moment and the next moment, respectively; After fracture occurs, the crack aperture is obtained by the displacement difference between the nodes at both ends of the interface contact spring: w=|(u1-u2)·n 12 | In the formula, u1 and u2 are the displacements of the connecting nodes at both ends of the spring, m and n, respectively. 12 Let be the unit normal vector of the contact surface, which is dimensionless.

Images