A geomechanically-constrained dynamic modeling method for geothermal reservoir fracture networks

Abstract

The present application belongs to the field of geothermal energy development and geological engineering technology, and discloses a geomechanics-constrained geothermal reservoir fracture network dynamic modeling method. The method comprises the following steps: first, obtaining multi-source monitoring data of the geothermal reservoir, including rock mass surface deformation data, micro-damage spatiotemporal distribution data and acoustic emission rupture signal data; constructing a micro-macro correlation model based on the bonded particle model and the tensor theory to obtain the eigenvalues of rock mass damage evolution; positioning the three-dimensional rupture source of the acoustic emission rupture signal; constructing a cubic cluster model based on the positioning results to generate the initial topological structure of the fracture network; dynamically adjusting the geometric shape and connectivity of the fracture network through the random walk algorithm and the fracture intersection correction mechanism according to the evolution eigenvalues and the geomechanics constraint conditions; and finally predicting the fracture expansion trend and damage evolution law under different geomechanics conditions. The present application significantly improves the accuracy and reliability of fracture network evolution prediction.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geothermal energy development and geological engineering, in particular to a method for dynamic modeling of geothermal reservoir fracture network based on geomechanical constraints. BACKGROUND

[0002] As a clean and renewable energy, geothermal energy has important strategic significance in global energy transformation. As the core carrier of geothermal energy development, the distribution characteristics and dynamic evolution law of the fracture network of geothermal reservoirs directly affect the efficiency and economy of resource exploitation. Due to the fact that geothermal reservoirs are usually located in deep high-temperature and high-pressure environments, the rock mass structure is complex, and the formation and evolution of the fracture network are strongly affected by geomechanical constraint conditions (such as stress field distribution, rock mass physical property heterogeneity, and fluid-solid coupling). Therefore, accurately modeling the dynamic evolution process of the fracture network and revealing the damage evolution law are key technical requirements for improving the efficiency of geothermal resource development and preventing and controlling engineering risks.

[0003] There are many bottlenecks in the existing technology in the field of dynamic modeling of geothermal reservoir fracture network. First, in terms of micro-macro correlation, there is a lack of systematic quantitative model between micro physical parameters (such as particle bond strength, contact stiffness, and friction coefficient) and macro mechanical parameters (such as stress-strain relationship and fracture strength), which makes it difficult to accurately reveal the evolution mechanism of micro damage inducing macro fracture, limits the prediction accuracy of reservoir damage evolution, and further affects the optimization of drilling design and stimulation measures. Secondly, in terms of fracture source positioning, the existing acoustic emission signal positioning methods (such as the traditional simplex method) have low accuracy, especially in reservoirs with significant acoustic velocity heterogeneity. For example, under high temperature and high pressure conditions, the inhomogeneity of rock acoustic velocity may cause a fracture source positioning error of up to tens of meters, directly affecting the accuracy of the spatial and temporal distribution characteristics of the fracture network.

[0004] In addition, there are obvious deficiencies in the fracture network modeling method. The existing topological structure construction process lacks comprehensive quantitative analysis of fracture source density, spacing and main direction. For example, in high stress concentration areas, the description of fracture connectivity is often too idealized, ignoring the geometric heterogeneity and directional distribution characteristics of the actual fracture network, resulting in a large deviation between the model and the real reservoir fracture distribution law. Further, the defects of fracture dynamic expansion simulation are particularly prominent: the existing methods often ignore the details of geomechanical constraints, such as not fully considering the dynamic changes of stress concentration coefficient and fracture surface friction resistance, and lack accurate description of fracture intersection behavior (such as fracture merging or termination conditions). This makes it difficult for the simulation results to truly reflect the reservoir evolution behavior, for example, in the process of hydraulic fracturing, the prediction deviation of fracture propagation path may cause stimulation failure or even seismic risk.

[0005] Meanwhile, the distortion problem of acoustic wave propagation model needs to be solved urgently. In the prior art, when the acoustic emission signal propagation time is calculated, the acoustic velocity is usually assumed to be uniformly distributed, and the complex influence of the acoustic velocity heterogeneity of the reservoir rock mass on the propagation path is ignored. For example, in a reservoir containing a multi-phase medium (such as a fluid-saturated fracture), the theoretical propagation path deviates significantly from the actual measured value, which seriously reduces the reliability of the fracture source positioning. Finally, the geometric feature statistical method is too simple. The fracture surface normal distribution is often based only on the direct average of the main direction of the fracture source, and lacks a detailed analysis of the direction distribution density. For example, in an anisotropic stress field, the statistical result is difficult to reflect the true fracture surface orientation characteristics, affecting the geometric expression accuracy of the fracture network model.

[0006] In summary, the prior art has significant deficiencies in micro-macro correlation, fracture source positioning, fracture network modeling, dynamic expansion simulation, and geometric feature description, which leads to weak prediction ability of geothermal reservoir damage evolution rules, and it is difficult to meet the high-precision needs of engineering development optimization and risk prevention and control.

[0007] Therefore, the present application provides a geomechanically constrained dynamic modeling method for geothermal reservoir fracture networks. SUMMARY

[0008] In order to overcome the defects of the prior art and achieve the above-mentioned purpose, the present application provides a geomechanically constrained dynamic modeling method for geothermal reservoir fracture networks. The method first acquires rock mass surface deformation data, mesoscopic micro-crack spatiotemporal distribution data, and acoustic emission fracture signal data of the geothermal reservoir under different geomechanical constraint conditions; constructs a correlation model of microphysical property parameters and macroscopic mechanical parameters based on the bonded particle model and the tensor of inertia, and obtains the evolution characteristic values of micro-damage-induced macro-fracture. Further, according to the mesoscopic micro-crack distribution data, the improved simplex method is used to perform three-dimensional fracture source positioning on the acoustic emission fracture signal, and the spatiotemporal distribution characteristics of the internal fracture source of the geothermal reservoir are extracted; based on the spatiotemporal distribution characteristics, a cubic cluster model is constructed to generate an initial topological structure of the fracture network.

[0009] After obtaining the initial topological structure and evolution characteristic values, the geometric shape and connectivity of the fracture network are dynamically adjusted in combination with the geomechanical constraint conditions to construct a dynamic evolution model. Specifically,

[0010] Based on the dynamic evolution model, the fracture network expansion trend and damage evolution law of the geothermal reservoir under different geomechanical constraints can be predicted. In constructing the correlation model of microphysical parameters and macroscopic mechanical parameters to obtain the evolution characteristic value, the microphysical parameters of the rock mass are first discretized, and the particle bond strength, contact stiffness and friction coefficient of each micro unit are extracted. Based on the bond particle model, the stress and strain tensor of each unit under the constraint condition is calculated to obtain the damage accumulation value. According to the decomposition of the damage accumulation value based on the moment tensor theory, the damage main direction and damage strength of the micro unit are extracted. The damage main direction distribution frequency and damage strength spatial heterogeneity of all units are counted to construct the spatial distribution characteristic matrix of micro damage. Based on the matrix and the macroscopic stress-strain relationship, the evolution characteristic value of the macroscopic fracture induced by micro damage is obtained through nonlinear mapping, including the critical stress threshold and the fracture expansion direction angle.

[0011] In the three-dimensional fracture source positioning using the improved simplex method, the time-frequency analysis of the acoustic emission fracture signal is performed to extract the arrival time, signal amplitude and main frequency characteristics. According to the arrival time, a time difference matrix between different monitoring points is constructed. The singular value decomposition of the matrix is performed to obtain the principal characteristic vector as the initial reference direction. Based on the simplex method, the reference direction is iteratively optimized. A tetrahedral search grid is constructed with the initial fracture source position as the center. The grid vertex position is dynamically adjusted according to the heterogeneity of acoustic wave velocity until the time difference residual is less than the preset threshold, and the accurate three-dimensional coordinates are output. Finally, based on all the fracture source coordinates and occurrence times, the spatiotemporal distribution characteristics are constructed, covering the dynamic changes of fracture source density distribution, expansion speed and direction.

[0012] In constructing the cubic cluster model, the fracture sources are divided into multiple local clusters (regions with density greater than a preset threshold) according to the spatiotemporal distribution characteristics. For each cluster, a cubic grid is constructed with the density center as the reference, and the number of fracture sources, spacing and main direction within the grid are counted. Based on the negative exponential decay function of fracture source spacing and the weighted calculation of cosine similarity of the main direction, the fracture connectivity probability is calculated. According to the probability, the cubic grid topology is connected to generate an initial topology structure containing the fracture geometry, connectivity and fracture plane normal distribution.

[0013] Further, the dynamic adjustment of the geometry and connectivity of the fracture network obtains a dynamic evolution model of the geothermal reservoir fracture network, including: determining the fracture expansion driving force of the geothermal reservoir fracture network under different geomechanical constraints according to the evolution characteristic value, the fracture expansion driving force including the stress concentration coefficient and the fracture surface friction resistance, the fracture expansion driving force calculation formula being:

[0014] ; wherein, is the fracture expansion driving force, is the stress concentration coefficient, for the crack length, for the friction coefficient, for the crack surface normal stress;

[0015] Based on the initial topology, a random walk algorithm is adopted: starting from the crack end point, preferentially selecting the direction with larger stress concentration coefficient and smaller friction resistance according to the probability distribution of driving force; introducing an intersection correction mechanism during the expansion process: when two crack paths intersect, if the normal angle is less than a preset threshold and the stress concentration coefficient difference is greater than a preset threshold, the crack merging is triggered; updating the network geometry and connectivity according to the expansion result, and outputting the dynamic evolution model.

[0016] The calculation of the acoustic wave velocity heterogeneity weighted time difference includes: generating an acoustic wave velocity distribution map by interpolating acoustic logging data; determining the theoretical propagation path from the rupture source to the monitoring point based on the map using the shortest path algorithm; calculating the path propagation time in segments and accumulating to obtain the theoretical arrival time; calculating the difference between the theoretical value and the actual arrival time as the weighted time difference. The acquisition of the crack surface normal distribution is obtained by projecting the main direction of the rupture source on the spherical surface in the cubic grid; generating a probability density function based on the density using kernel density estimation; determining the maximum probability direction as the crack surface normal; and constructing the crack network normal distribution by statistically all grid normals.

[0017] The technical effects and advantages of the geomechanics-constrained dynamic modeling method for geothermal reservoir crack network based on the present application are:

[0018] The present application effectively solves the traditional problem of disconnection between micro-mechanism and macro-performance through multi-scale data fusion and cross-scale mechanical correlation, significantly improves the accuracy and reliability of crack network evolution prediction. The dynamic evolution model constructed accurately describes the response characteristics of geothermal reservoirs under complex geomechanical environment, providing a scientific basis for efficient resource development. The method significantly improves the reservoir evaluation accuracy and reliability by accurately capturing the spatio-temporal evolution law of cracks, and reduces the exploration and development risk; its adaptive optimization characteristics can flexibly adapt to different geological conditions and engineering needs, significantly enhancing the adaptability and robustness of practical application. BRIEF DESCRIPTION OF DRAWINGS

[0019] Figure 1 A geomechanics-constrained dynamic modeling method for geothermal reservoir crack network based on the present application is shown in the figure. DETAILED DESCRIPTION

[0020] The embodiments of the present application are described clearly and completely in conjunction with the drawings, and the described embodiments are only part of examples and not all. Based on the embodiments of the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection.

[0021] The application provides a geomechanical constraint-based geothermal reservoir fracture network dynamic modeling method. The method obtains multi-source data of a geothermal reservoir under different geomechanical constraint conditions, constructs a correlation model of microphysical property parameters and macroscopic mechanical parameters, locates a breakage source by using an improved simplex method, constructs a cubic cluster model to generate an initial topological structure, and predicts a fracture network expansion trend and damage evolution law in combination with dynamic evolution characteristics. The method has high precision and multi-scale fusion characteristics, can accurately depict the dynamic evolution process of a fracture network under geomechanical constraints, and provides key technical support for geothermal energy development.

[0022] In the embodiment of the application, the detailed implementation steps of the geomechanical constraint-based geothermal reservoir fracture network dynamic modeling method include:

[0023] First, the rock mass surface deformation data, mesoscopic microcrack spatiotemporal distribution data and acoustic emission breakage signal data are obtained. The surface deformation data are collected in real time by a geological monitoring device and contain key indicators such as rock mass surface displacement and strain under pressure / temperature conditions; the mesoscopic microcrack data are obtained by high-resolution imaging technology and microscopic analysis and cover microcrack position, density distribution and geometric characteristics; and the acoustic emission breakage signal data are collected by an acoustic emission monitoring system and reflect elastic wave information released by internal breakage of the rock mass. The above data constitute the basis for a fracture network evolution model and ensure the accuracy and applicability of the model.

[0024] A micro-macro correlation model is constructed based on the bonded particle model and the moment tensor theory to obtain evolution characteristic values of micro-damage-induced macro-breakage. The model establishes a quantitative relationship between micro-rock mass physical property changes and macro-mechanical performance, reveals the cross-scale damage evolution mechanism, and quantifies the critical conditions for macro-breakage caused by micro-damage accumulation to provide key parameters for fracture expansion prediction. The model construction process combines theoretical analysis and experimental verification to ensure scientificity and accuracy. Specifically, the model construction process includes:

[0025] Based on the mesoscopic microcrack spatiotemporal distribution data, an improved simplex method is used to locate three-dimensional breakage sources of acoustic emission breakage signals to obtain the spatiotemporal distribution characteristics of internal breakage sources of the geothermal reservoir. The improved simplex method is a high-precision optimization algorithm and can effectively solve the nonlinear problem in breakage source positioning. The extracted spatiotemporal distribution characteristics reflect the evolution law of internal breakage activities of the reservoir and cover key information such as spatial density distribution, time activity frequency and expansion direction of the breakage sources, thereby providing a data basis for subsequent construction of a fracture network topological structure.

[0026] Based on the spatiotemporal distribution characteristics, the initial topological structure of the fracture network is generated by constructing a cubic cluster model. The cubic cluster model converts the discrete fracture source data into a continuous fracture network structure through a spatial clustering method. The initial topological structure completely describes the geometric shape, connectivity and spatial distribution characteristics of the fracture network, laying a foundation for dynamic evolution simulation. The model construction process fully considers the density and direction characteristics of the fracture source, ensuring the representativeness and integrity of the generated structure.

[0027] Based on the initial topological structure, evolution characteristic values and geomechanical constraint conditions, the geometric shape and connectivity of the fracture network are dynamically adjusted to construct a dynamic evolution model. This model accurately depicts the temporal evolution process of the fracture network under geomechanical constraints, including fracture expansion, closure and intersection, etc. The model takes the initial topological structure as the starting point, simulates the response of fractures to different geomechanical conditions through a physical driving evolution algorithm, and realizes dynamic prediction of the development of the fracture network.

[0028] Based on the dynamic evolution model, the expansion trend and damage evolution law of the fracture network under different geomechanical constraint conditions are predicted. The prediction results cover key information such as the spatial development direction of the fracture, the connectivity change and the evolution of the fracture density, providing a scientific basis for engineering decision-making in geothermal resource development. The prediction process uses a multi-scenario analysis method to evaluate the response characteristics of the fracture network under different geomechanical constraints, ensuring the comprehensiveness and applicability of the prediction results.

[0029] In the embodiment of the present application, the correlation model for constructing the microphysical property parameters and macroscopic mechanical parameters of the geothermal reservoir rock mass, and obtaining the evolution characteristic values of the micro-damage induced macro-fracture, comprises:

[0030] First, the microphysical property parameters are discretized to extract the particle bond strength, contact stiffness and friction coefficient of each micro unit. The discretization process converts the continuous rock mass medium into a set of discrete micro units, facilitating micro-mechanism simulation and calculation. The particle bond strength represents the bonding ability between particles, the contact stiffness reflects the deformation response characteristics, and the friction coefficient describes the resistance to relative motion of particles. These parameters directly determine the rock mass fracture characteristics and energy release mode, and constitute the basic data support for micro-modeling. Discretization uses multi-scale analysis technology to balance the representativeness and calculation efficiency of micro units.

[0031] Based on the bonded particle model, the stress tensor and strain tensor of each micro unit under geomechanical constraint conditions are calculated to obtain the damage accumulation value of the micro unit. The stress tensor and strain tensor represent the mechanical state of the micro unit, and the damage accumulation value quantifies the damage degree of the micro structure. The bonded particle model can accurately describe the non-continuous deformation and fracture process of rock mass, and the calculation process considers stress field, temperature field and fluid pressure constraints, etc., to ensure the authenticity and integrity of the simulation results.

[0032] According to the tensor decomposition of the damage accumulation value based on the moment tensor theory, the damage principal direction and the damage intensity of the micro unit are extracted. The damage principal direction represents the dominant propagation path of micro cracks, and the damage intensity quantifies the damage degree of the micro structure. The moment tensor theory simplifies the complex damage state into the combination of direction and intensity by decomposing the anisotropic damage through the eigenvalue analysis method, which is convenient for subsequent analysis and application.

[0033] The distribution frequency of the damage principal direction and the spatial heterogeneity of the damage intensity of all micro units are counted, and a spatial distribution characteristic matrix of micro damage is constructed. The distribution frequency describes the concentration degree of the damage principal direction, and the spatial heterogeneity reflects the change characteristics of the damage intensity. The matrix serves as a statistical expression of the micro damage state and characterizes the distribution law of the damage in space. The matrix construction adopts a spatial statistical method to integrate discrete damage information into continuous spatial characteristics, providing data support for micro-macro correlation analysis.

[0034] Based on the spatial distribution characteristic matrix and the macro stress-strain relationship, the evolution characteristic values of micro damage induced macro fracture are obtained by using a nonlinear mapping method. The nonlinear mapping method captures the complex correlation between micro damage and macro fracture, and the evolution characteristic values include a critical stress threshold and a fracture propagation direction angle. The critical stress threshold defines the critical condition for micro damage to trigger macro fracture, and the direction angle describes the directional characteristics of fracture propagation. The characteristic values are obtained by combining theoretical analysis and experimental data fitting, establishing a quantitative bridge between micro and macro scales, and providing key parameters for the dynamic evolution of the fracture network.

[0035] In the embodiment of the present application, the improved simplex method is used to perform three-dimensional fracture source positioning on the acoustic emission fracture signals, and the spatio-temporal distribution characteristics of the fracture sources inside the geothermal reservoir are obtained, including:

[0036] Firstly, time-frequency analysis is performed on the acoustic emission fracture signals to extract the arrival time, amplitude and main frequency characteristics of each signal. As a basic step of signal processing, time-frequency analysis extracts time domain and frequency domain characteristics through techniques such as short-time Fourier transform and wavelet analysis: the arrival time is a key basis for positioning the fracture source, the signal amplitude represents the fracture energy intensity, and the main frequency characteristics are related to the fracture mechanism type. The extracted characteristics provide complete data support for subsequent positioning algorithms.

[0037] Based on the arrival time, a time difference matrix between different monitoring points is constructed. As the core data structure for positioning, the matrix records the time difference information of the arrival of sound waves at each monitoring point, and its construction process fully considers the geometric layout of the monitoring network: for each pair of monitoring points, the arrival time difference of sound waves is calculated, and a complete matrix is formed to reflect the geometric characteristics of the sound wave propagation path, which provides a direct basis for subsequent positioning.

[0038] The singular value decomposition is performed on the time difference matrix to obtain a principal eigenvector as a reference direction for initial fracture source positioning. The singular value decomposition extracts the main variation direction of the data through mathematical dimension reduction, and the principal eigenvector reveals the spatial trend of the time difference distribution, which is directly related to the fracture source position, significantly improving the efficiency and accuracy of subsequent optimization.

[0039] The reference direction is corrected by iterative optimization based on the simplex method: a tetrahedral search grid is constructed with the initial fracture source position as the center; the time difference is dynamically adjusted according to the heterogeneity of the acoustic velocity of the grid vertices; until the time difference residual is less than the preset threshold, the accurate three-dimensional coordinates are output. As an efficient nonlinear optimization algorithm, the simplex method adapts to complex acoustic environments with its tetrahedral search structure, and effectively compensates for the effects of medium non-uniformity with its acoustic velocity heterogeneity weighting mechanism, gradually approaching the true fracture source position in the iterative process.

[0040] Integrate all three-dimensional coordinates of the fracture sources and the corresponding occurrence times to construct a spatiotemporal distribution feature representing the evolution of internal fracture activity in the reservoir: the density distribution describes the spatial concentration of the fracture sources, the extension velocity quantifies the activity intensity in the time dimension, and the extension direction reflects the fracture development trend. This feature connects micro-fracture events with macro-fracture structures through spatiotemporal clustering analysis, providing a core basis for the construction of the fracture network.

[0041] In the embodiment of the present application, the cubic cluster model of the acoustic emission fracture signal is constructed to obtain the initial topological structure of the geothermal reservoir fracture network, which comprises:

[0042] When constructing the cubic cluster model, first, the fracture sources inside the geothermal reservoir are divided into multiple local fracture source clusters according to the spatiotemporal distribution feature. A local fracture source cluster is a spatial region with a fracture source density greater than a preset density threshold, which essentially is a set of spatially associated fracture events, reflecting the local characteristics of fracture development. The division process uses a density clustering algorithm to identify high-density potential fracture development areas, and the preset density threshold is set to the average density plus twice the standard deviation based on statistical analysis, ensuring that the cluster structure has statistical significance. This process converts discrete fracture source data into a cluster structure with clear physical meaning, laying the foundation for topological construction.

[0043] For each local fracture source cluster, a cubic grid is constructed with the fracture source density center as the reference. The cubic grid is a regular spatial discrete structure that facilitates topological analysis, and the density center is a geometric reference point. The number of fracture sources, the distance between fracture sources, and the main direction of fracture sources are obtained through spatial statistical methods: the number of fracture sources represents the fracture activity intensity, the distance between fracture sources describes the distribution density, and the main direction of fracture sources indicates the fracture development trend. These characteristics quantify the local area fracture characteristics and provide a basis for fracture connectivity analysis.

[0044] The fracture connectivity probability of each cubic grid is calculated based on the distance between fracture sources and the main direction, and the calculation formula is:

[0045] ; wherein, is the crack connectivity probability, is the distance between crack sources, is the attenuation coefficient, is the included angle between the main directions of the crack sources.

[0046] The probability quantifies the possibility of adjacent crack sources forming continuous cracks, which is a key parameter for topology construction. The negative exponential decay function in the formula reflects the influence of distance on connectivity, and the closer the distance, the higher the probability. The cosine similarity function reflects the promoting effect of direction consistency, and the more consistent the direction, the higher the probability. The weighted combination of the two comprehensively considers the spatial geometry and direction characteristics to generate a reasonable probability estimate.

[0047] Based on the crack connectivity probability, all cubic grids are topologically connected: when the connectivity probability between grids exceeds a certain threshold, a connection is established to form a preliminary network skeleton. This process integrates discrete grids into a continuous crack network, generating an initial topological structure that includes geometric morphology, connectivity, and crack surface normal distribution. Geometric morphology describes the spatial shape and size of the crack, connectivity reflects the channel characteristics of the network, and crack surface normal distribution represents the directionality of the crack space. This structure is extracted through graph theory and network analysis methods and serves as the starting model for dynamic evolution, which is the core link of crack network modeling.

[0048] In the embodiment of the present application, the dynamic adjustment of the geometric morphology and connectivity of the crack network to obtain the dynamic evolution model of the geothermal reservoir crack network comprises:

[0049] According to the evolution characteristic value, the crack propagation driving force of the geothermal reservoir crack network under different geomechanical constraint conditions is determined, and the crack propagation driving force includes a stress concentration coefficient and a crack surface friction resistance.

[0050] The crack propagation driving force calculation formula is:

[0051] ; wherein, is the crack propagation driving force, is the stress concentration coefficient, is the crack length, is the friction coefficient, is the crack surface normal stress.

[0052] The driving force, as the core mechanism of crack dynamic evolution, determines the expansion rate and direction: KI KIThe stress amplification effect at the fracture tip is reflected, promoting fracture propagation; frictional resistance characterizes the sliding resistance at the fracture surface, inhibiting fracture activity. The process of determining the driving force combines evolutionary characteristic values ​​with geomechanical theory, integrating factors such as rock mass properties, stress state, and temperature field, to provide a physical driving basis for propagation simulation.

[0053] Based on the initial topology, a random walk algorithm is used to simulate the dynamic propagation process of the crack. The random walk algorithm includes: starting from the crack endpoint, determining the probability distribution of the propagation direction according to the crack propagation driving force, and preferentially selecting the direction with a larger stress concentration coefficient and smaller crack surface friction resistance for propagation.

[0054] The formula for calculating the probability distribution of the expansion direction in the random walk algorithm is:

[0055] ;

[0056] in, For direction The expansion probability, The normalization constant is This is the stress concentration factor in that direction. This represents the frictional resistance in that direction.

[0057] Random walk algorithms are efficient path generation methods suitable for simulating the random and directional propagation process of cracks. The crack endpoints are locations of stress concentration and are also the most likely starting points for crack propagation. The probability distribution reflects the likelihood of propagation in different directions, transforming physical driving forces into mathematical probabilities and achieving stochastic simulation of the physical process. The priority selection principle ensures that the simulation process conforms to the laws of mechanics, maintaining randomness while reflecting the characteristics of physical driving forces, making the simulation results more realistic and reliable.

[0058] A fracture convergence correction mechanism is introduced during fracture propagation. This mechanism includes: when two fracture paths intersect, the behavior is determined based on the difference between their normal angle and stress concentration factor. If the normal angle is less than a preset angle threshold and the stress concentration factor difference is greater than a preset difference threshold, fracture merging is triggered. This mechanism simulates the topological changes in fracture network interactions: the normal angle characterizes the spatial directional relationship of fractures, and the difference in stress concentration factor reflects the imbalance of propagation driving forces. When two fracture directions are close and the driving forces differ significantly, the fracture with the stronger driving force dominates the propagation process, triggering merging, significantly improving the physical realism and prediction accuracy of the model.

[0059] The dynamic evolution model is generated by updating the geometry and connectivity of the fracture network according to the extension result. The model completely describes the process of the fracture network evolution over time, covers the fracture length growth, new connection formation and topological structure change, and records the network state at different time points to analyze the evolution law. The model combines the initial topological structure and the physical driving extension process to form a complete description of the fracture network spatio-temporal evolution, and provides a core theoretical basis for geothermal reservoir evaluation and development.

[0060] In the embodiment of the present application, the method for calculating the acoustic velocity heterogeneity weighted time difference comprises:

[0061] The calculation of the acoustic velocity heterogeneity weighted time difference first needs to obtain the acoustic velocity distribution map of the geothermal reservoir rock mass. The distribution map is obtained by interpolating acoustic logging data, and the acoustic logging data is obtained by the logging tool during drilling. The acoustic velocity distribution map is a spatial expression of the acoustic characteristics of the rock mass, reflecting the influence of medium heterogeneity on acoustic propagation. The interpolation process uses geostatistical methods, such as Kriging interpolation, to convert discrete logging data into continuous velocity field distribution, while considering the spatial correlation and anisotropy characteristics of the formation. The resolution of the acoustic velocity distribution map usually reaches the meter level, which can effectively capture the velocity variation characteristics in the rock mass and provide basic data for subsequent acoustic path analysis.

[0062] Based on the above acoustic velocity distribution map, the theoretical propagation path of the acoustic emission fracture signal from the fracture source to each monitoring point is calculated. The theoretical propagation path reflects the actual propagation trajectory of the acoustic wave in the non-uniform medium, and is determined by using the shortest path algorithm. The algorithm is based on Fermat's principle and aims to find the path with the shortest propagation time in a non-uniform velocity field. In the implementation process of the algorithm, the continuous velocity field is discretized into a grid structure, the propagation time between grid nodes is calculated, and the Dijkstra algorithm or A* algorithm is used to search for the global optimal path. This method fully considers the influence of medium heterogeneity on acoustic propagation, overcomes the limitations of the traditional straight-line path assumption, and significantly improves the accuracy of path calculation.

[0063] Subsequently, according to the determined acoustic velocity distribution map and the theoretical propagation path, the propagation time of the acoustic emission fracture signal in different path segments is calculated. The propagation time of all path segments is accumulated, and the theoretical arrival time of the acoustic emission fracture signal to each monitoring point is obtained. The theoretical arrival time is the theoretical reference value for acoustic source positioning. The specific calculation process is to decompose the theoretical propagation path into several path segments, calculate the segment propagation time according to the length of each path segment and the acoustic velocity corresponding to the position of the path segment, and finally perform accumulation summation. This segmented calculation method fully considers the spatial variation of acoustic velocity in the propagation process, improves the accuracy of time calculation, and lays a reliable foundation for subsequent time difference calculation.

[0064] Finally, the difference between the theoretical arrival time and the actual arrival time is calculated, and this difference is defined as the sound wave velocity heterogeneity-weighted time difference. This weighted time difference is a corrected time difference that takes into account the influence of medium inhomogeneity, and it can more accurately reflect the actual characteristics of sound wave propagation. By directly comparing the theoretical value and the measured value, the difference calculation quantifies the degree of influence of velocity heterogeneity on propagation time. This weighted time difference, as an input parameter of the positioning algorithm, replaces the traditional direct time difference, effectively compensating for the error caused by medium inhomogeneity, thereby significantly improving positioning accuracy. The introduction of the weighted time difference is an important innovation of this method, which enables the positioning algorithm to adapt to complex geological environments and has a wider range of applicability.

[0065] In this embodiment of the invention, the method for obtaining the normal distribution of the fracture surface includes:

[0066] First, the principal directions of fracture origins within each cubic mesh are spherically projected to obtain their distribution density on the sphere. Spherical projection, as an effective tool for analyzing three-dimensional directional data, can intuitively display the spatial distribution characteristics of directions. The principal directions of fracture origins reflect the dominant directions of microcracks and are directly related to the directions of macroscopic fracture surfaces. The projection process maps three-dimensional direction vectors onto a unit sphere to form a spherical distribution. The distribution density is obtained by counting the number of points in different regions on the sphere; this density reflects the concentration of the directional distribution. This step transforms discrete directional data into a continuous density distribution, providing a data foundation for subsequent quantitative analysis and normal determination.

[0067] Subsequently, based on this distribution density, the probability density function of the main direction of the rupture source is obtained using the kernel density estimation method. Kernel density estimation is a non-parametric density estimation technique that can extract continuous distribution characteristics from data. The probability density function quantifies the likelihood of occurrence in different directions and is a mathematical expression of the directional distribution. The estimation process uses the Fisher kernel function, which is particularly suitable for spherical directional data analysis and can effectively capture the central tendency and discrete characteristics of the directional distribution. The smoothing parameter of the probability density function is adaptively adjusted according to the amount of data and distribution characteristics to ensure the reliability and representativeness of the estimation results. The formula for calculating the probability density function of the main direction of the rupture source is:

[0068] ;in, For direction The probability density at that location, For the sample size, For lumped parameters, For the first The direction vector of the rupture source.

[0069] This formula transforms discrete projection points into a continuous probability distribution, providing a theoretical basis for determining the direction of maximum probability.

[0070] Next, according to the probability density function, the maximum probability direction of the main direction of the fracture source is determined, and is taken as the normal of the fracture surface in the corresponding cubic grid. The maximum probability direction is the global maximum point of the probability density function, representing the dominant trend of the direction distribution. The normal determination process uses a numerical optimization method to search for the maximum point of the probability density function on the sphere to obtain the corresponding direction vector. If there are multiple local maxima, the globally maximum direction is selected as the main normal, and the significant local maximum can be recorded as the secondary normal. This statistical-based normal determination method comprehensively considers all fracture source information, has high stability and representativeness, and can accurately reflect the spatial direction characteristics of the fracture surface.

[0071] Finally, the normals of the fracture surfaces of all cubic grids are counted, and the normal distribution of the fracture surface of the geothermal reservoir fracture network is constructed. The normal distribution of the fracture surface is the expression of the directionality characteristics of the entire fracture network, reflecting the spatial direction regularity of the reservoir fractures. The statistical process analyzes the normal data of all grids, and calculates the characteristic parameters such as the concentration of the main direction and the dispersion of the direction by using the method of directional statistics. The distribution results are usually visualized in the way of polar azimuthal projection or equal-area projection, directly showing the spatial direction characteristics of the fracture surface. This distribution information is of great significance for understanding the stress state, tectonic background and fluid flow characteristics of the reservoir, and is a key reference data for geothermal resource assessment and development.

[0072] The present application realizes high-precision dynamic modeling of geothermal reservoir fracture network through the technical route of multi-source data acquisition, micro-macro correlation model construction, improved simplex positioning, cubic cluster modeling and dynamic evolution adjustment. Its geomechanical constraint characteristics can accurately reflect the response law of the geothermal reservoir under different conditions, and provide key technical support for reservoir evaluation, injection-production optimization and production prediction in geothermal energy development and utilization. The above formulas are based on dimensionless numerical calculation, and are obtained by software simulation fitting of a large amount of data. The preset parameters and thresholds in the formulas are set by the person skilled in the art according to the actual situation. The above is the preferred embodiment of the present application, and the protection scope is not limited to the above examples. Any technical solution within the idea of the present application is protected. Improvements and refinements of the person skilled in the art without departing from the principles of the present application should also be considered as the protection scope of the present application.

Claims

1. A dynamic modeling method for geothermal reservoir fracture networks based on geomechanical constraints, characterized in that, include: Acquire data on rock surface deformation, spatiotemporal distribution of microcracks, and acoustic emission fracture signals of geothermal reservoirs under different geomechanical constraints. Based on the bonded particle model and moment tensor theory, a correlation model between the microscopic physical parameters and macroscopic mechanical parameters of the geothermal reservoir rock mass is constructed, and the evolution characteristic values ​​of macroscopic fracture induced by microscopic damage are obtained. According to the spatiotemporal distribution data of the micro-cracks, the improved simplex method is used to locate the three-dimensional fracture source of the acoustic emission fracture signal, so as to obtain the spatiotemporal distribution characteristics of the fracture source inside the geothermal reservoir. Based on the aforementioned spatiotemporal distribution characteristics, a cubic cluster model is constructed to generate the initial topology of the fracture network. Combining the initial topology, evolutionary characteristic values, and geomechanical constraints, a dynamic evolution model is obtained by dynamically adjusting the geometry and connectivity of the fracture network. Based on the dynamic evolution model, the expansion trend and damage evolution law of the fracture network under different geomechanical constraints are predicted. The model establishing the correlation between the microscopic physical parameters and macroscopic mechanical parameters of the geothermal reservoir rock mass is used to obtain the evolutionary characteristic values ​​of macroscopic fracturing induced by microscopic damage, including: The microscopic physical parameters of the geothermal reservoir rock mass are discretized to obtain the particle bonding strength, interparticle contact stiffness, and interparticle friction coefficient of each micro-unit. Based on the bonded particle model, the stress tensor and strain tensor of each micro-unit under different geomechanical constraints are calculated to obtain the cumulative damage value of the micro-unit. The cumulative damage value is decomposed according to the moment tensor theory to extract the principal damage direction and damage intensity of the micro-unit. The spatial heterogeneity of the distribution frequency of the principal damage direction and the damage intensity of all micro-units are statistically analyzed to construct a spatial distribution feature matrix of micro-damage. Based on the stress-strain relationship between the spatial distribution feature matrix and the macroscopic mechanical parameters, a nonlinear mapping method is used to obtain the evolutionary feature values ​​of micro-damage-induced macroscopic fracture. The evolutionary feature values ​​include the critical stress threshold for damage-induced fracture and the deflection angle of the fracture propagation direction.

2. The method for dynamic modeling of geothermal reservoir fracture networks based on geomechanical constraints according to claim 1, characterized in that, The improved simplex method is used to locate the three-dimensional fracture source of acoustic emission fracture signals and obtain the spatiotemporal distribution characteristics of fracture sources inside the geothermal reservoir. The location includes: Time-frequency analysis is performed on acoustic emission rupture signals to obtain the arrival time, signal amplitude, and dominant frequency characteristics of each signal. Based on the arrival time, a time difference matrix is ​​constructed between different monitoring points for the acoustic emission rupture signals. Singular value decomposition is performed on the time difference matrix to obtain its principal eigenvector, which serves as the reference direction for initial rupture source localization. Based on the simplex method, an iterative optimization approach is used to correct the reference direction. This iterative optimization includes: constructing a tetrahedral search grid centered on the initial rupture source location; dynamically adjusting the tetrahedral vertex positions based on the time difference weighted by the acoustic velocity heterogeneity of the tetrahedral vertices until the time difference residual is less than a preset residual threshold, thus obtaining accurate three-dimensional coordinates of the rupture source. Based on the three-dimensional coordinates of all rupture sources and the corresponding acoustic emission rupture signal occurrence times, the spatiotemporal distribution characteristics of rupture sources within the geothermal reservoir are constructed. These spatiotemporal distribution characteristics include the dynamic changes in the density distribution, propagation velocity, and propagation direction of the rupture sources.

3. The method for dynamic modeling of geothermal reservoir fracture networks based on geomechanical constraints according to claim 1, characterized in that, The process of constructing a cube cluster model to generate the initial topology of the fracture network includes: Based on the spatiotemporal distribution characteristics, the fracture sources within the geothermal reservoir are divided into multiple local fracture source clusters, where each local fracture source cluster is defined as a spatial region where the fracture source density is greater than a preset density threshold. For each local fracture source cluster, a cubic grid is constructed with the fracture source density center as the reference, and the number of fracture sources, the distance between fracture sources, and the main direction of the fracture sources within each grid are obtained. Based on the distance between fracture sources and the main direction of the fracture sources, the fracture connectivity probability of each cubic grid is calculated using a weighted average of a negative exponential decay function and cosine similarity. Based on the fracture connectivity probability, all cubic grids are topologically connected to generate the initial topology of the geothermal reservoir fracture network. The initial topology includes fracture geometry, connectivity, and fracture surface normal distribution.

4. The method for dynamic modeling of geothermal reservoir fracture networks based on geomechanical constraints according to claim 1, characterized in that, The method of obtaining a dynamic evolution model by dynamically adjusting the geometry and connectivity of the fracture network includes: The driving force for fracture propagation is determined based on the aforementioned evolutionary characteristic values. This driving force includes the stress concentration factor and the frictional resistance of the fracture surface. The formula for calculating the driving force for fracture propagation is as follows: ;in, As the driving force for crack propagation, The stress concentration factor is... The crack length is... The coefficient of friction, The stress is the normal stress on the fracture surface; Based on the initial topology, a random walk algorithm is used to simulate the dynamic propagation of fractures. The random walk algorithm includes: starting from the fracture endpoint, determining the probability distribution of the propagation direction based on the propagation driving force; preferentially selecting the direction with a larger stress concentration coefficient and smaller frictional resistance for propagation; introducing a fracture convergence correction mechanism during the propagation process, which includes: when the propagation paths of two fractures intersect, determining whether fracture merging or fracture termination occurs based on the difference between the normal angle and the stress concentration coefficient of the two fractures. The determination criteria are that fracture merging occurs when the normal angle is less than a preset angle threshold and the difference in stress concentration coefficient is greater than a preset difference threshold; and updating the geometry and connectivity of the fracture network based on the dynamic propagation results of the fractures to obtain a dynamic evolution model of the geothermal reservoir fracture network.

5. The method for dynamic modeling of geothermal reservoir fracture networks based on geomechanical constraints according to claim 2, characterized in that, The calculation of the sound wave velocity heterogeneity weighted time difference includes: A sonic velocity distribution map of the geothermal reservoir rock mass is generated by interpolating sonic logging data. The sonic velocity distribution map is obtained by interpolating sonic logging data. Based on the sonic velocity distribution map, the theoretical propagation path from the fracture source to each monitoring point is determined by the shortest path algorithm. The propagation time of the acoustic emission signal is calculated segment by segment along the theoretical propagation path, and the theoretical arrival time is obtained by accumulating the segments. The difference between the theoretical arrival time and the actual arrival time is calculated as the sonic velocity heterogeneity weighted time difference.

6. The method for dynamic modeling of geothermal reservoir fracture networks based on geomechanical constraints according to claim 3, characterized in that, The method for obtaining the normal distribution of the fracture surface includes: For each cubic grid, the principal direction of the fracture source is spherically projected to obtain the distribution density of the principal direction of the fracture source on the sphere. Based on the distribution density, the kernel density estimation method is used to obtain the probability density function of the principal direction of the fracture source. Based on the probability density function, the direction with the highest probability of the principal direction of the fracture source is determined as the normal of the fracture surface in the corresponding cubic grid. The fracture surface normals of all cubic grids are statistically analyzed to construct the fracture surface normal distribution of the geothermal reservoir fracture network.

Images