A discrete element based method for hydraulic fracturing and microseismic / acoustic emission simulation

Abstract

The application discloses a hydraulic fracturing and microseismic / acoustic emission simulation method based on discrete elements, realizes quantitative design and in-situ characterization of natural discrete fracture network (DFN) based on a power-law distribution function, solves the problem that the existing method cannot accurately depict the guiding and controlling effect of DFN on hydraulic fracture propagation, constructs a fluid-solid coupling algorithm based on Dashi's law and the cubic law, realizes full-dynamic interactive simulation of fluid migration, pore pressure evolution, stress redistribution and rock mass bonding failure, solves the problem that the existing seepage and fracture evolution cannot be dynamically and synchronously coupled, constructs a virtual microseismic monitoring module which is synchronously operated with the hydraulic fracturing simulation, realizes multi-dimensional analysis of microseismic / acoustic emission events based on the moment tensor theory, solves the problem that the existing method cannot synchronously extract and quantitatively explain the microseismic / acoustic emission signals from the discrete element model, thereby realizing effective benchmarking of simulation results with laboratory test and field monitoring data, and improving the accuracy and construction safety of hydraulic fracturing design of deep reservoirs.

Description

Technical Field

[0001] This invention belongs to the field of deep resource extraction and rock mechanics technology, specifically involving numerical simulation and microseismic / acoustic emission response analysis of rock hydraulic fracturing process under multi-field coupling of heat, fluid, solid, and chemical fields. In particular, it relates to a simulation method based on discrete element method coupled with real-time monitoring of fluid seepage, rock mechanics fracturing, and microseismic / acoustic emission. This method can be used to predict and analyze fracture propagation, permeability evolution, and microseismic / acoustic emission event distribution in brittle reservoirs under high-pressure fluid injection in scenarios such as deep coalbed methane, shale gas, dry hot rock geothermal reservoir stimulation, and carbon dioxide geological storage. Background Technology

[0002] Unconventional resource development projects such as deep coalbed methane extraction, shale gas development, and hot dry geothermal reservoir stimulation all require the use of hydraulic fracturing technology to construct artificial fracture networks and effectively improve reservoir permeability. However, deep reservoirs exist in complex, multi-field environments characterized by high geostress, high pore pressure, and high temperature. The initiation, propagation, and connection behavior of hydraulic fractures are significantly influenced by the natural fracture network, resulting in a complex evolution process that is difficult to predict accurately.

[0003] Among existing numerical simulation methods, the continuum finite element method (FEM) and finite difference method (FDM) have limited ability to describe the discontinuous fracture behavior of rock masses, and cannot accurately characterize the dynamic propagation process of fractures and the guiding role of natural fractures. Although the traditional discrete element method (DEM) can effectively simulate the discontinuous failure of rock masses, it still has significant shortcomings in terms of fully coupled seepage-mechanics, accurate characterization of natural fracture networks, and simultaneous extraction and quantitative interpretation of microseismic / acoustic emission signals. Firstly, existing discrete element hydraulic fracturing models cannot accurately reproduce the in-situ geometry and forces of natural fractures. First, existing fluid-structure interaction algorithms cannot accurately quantify the impact of natural fracture networks on the propagation path of hydraulic fractures due to their inherent characteristics. Second, they cannot achieve full dynamic synchronous interaction of fluid seepage, pore pressure evolution, rock stress redistribution, and bond failure, resulting in insufficient accuracy in multi-field coupling. Third, existing technologies cannot achieve full-step synchronous operation of microseismic / acoustic emission simulation and hydraulic fracturing process, and can only perform statistical analysis of bond failure events afterward. They cannot complete real-time analysis and quantitative interpretation of the source mechanism, and it is difficult to effectively benchmark against microseismic / acoustic emission data from indoor tests and field monitoring.

[0004] Therefore, developing a comprehensive simulation method that can simultaneously achieve accurate characterization of natural fracture networks, full-process simulation of hydraulic fracturing dynamics, high-precision coupling of seepage-mechanics multi-fields, and synchronous extraction and quantitative interpretation of microseismic / acoustic emission signals is of great engineering significance and application value for guiding the safe and efficient exploitation of deep resources. Summary of the Invention

[0005] To address the problems existing in the prior art, this invention provides a discrete element method for hydraulic fracturing and microseismic / acoustic emission simulation. It achieves quantitative design and in-situ characterization of natural discrete fracture networks based on a power-law distribution function. A fluid-structure interaction algorithm based on Darcy's law and the cubic law is constructed, and a virtual microseismic monitoring module runs synchronously with the full-step hydraulic fracturing simulation. This enables effective benchmarking of simulation results with indoor experimental and field monitoring data, balancing the simulation requirements of laboratory rock sample scale and field engineering scale, and improving the accuracy and construction safety of deep reservoir hydraulic fracturing design.

[0006] To achieve the above objectives, the technical solution adopted by this invention is: a hydraulic fracturing and microseismic / acoustic emission simulation method based on discrete element method, specifically including the following steps: Step 1: Based on the bonded particle model of the Particle Flow Program (PFC), establish a discrete element fluid-solid coupled hydraulic fracturing model and calibrate the model input parameters.

[0007] Step 2: Design a natural discrete fracture network (DFN) based on the power-law distribution function and integrate the network into the calibrated discrete element fluid-structure coupled hydraulic fracturing model to characterize the in-situ geometry of rock fractures; the model supports dynamic adjustment of the relative modulus of fractures and the fracture criterion to describe the entire process of hydraulic fracture initiation, propagation and penetration.

[0008] Step 3: Based on the microscopic parameters after model calibration, construct a flow skeleton network composed of particle contact pores to clarify the correspondence between fluid flow channels and pore space.

[0009] Step 4: Set the model's ground stress boundary conditions and injection boundary conditions, start the iterative injection simulation program, use water as the fracturing fluid, and inject it into the pre-specified central borehole of the model according to the preset fixed flow rate.

[0010] Step 5: Calculate the fluid pore pressure in the flow region based on Darcy's law and the cubic law. Couple the pore pressure with the geostress and particle bonding stress to calculate the effective stress. Simultaneously update the pore space volume and fluid flow parameters to simulate the dynamic interaction process of fluid transport and rock mechanical response at the pore scale.

[0011] Step Six: Configure a virtual microseismic monitoring module that runs synchronously with the hydraulic fracturing simulation program throughout the entire process. This module monitors and records interparticle bonding failure events in real time during each calculation step, and simultaneously completes the source location, source mechanism analysis, and magnitude calculation of microseismic / acoustic emission events.

[0012] Step 7: When the model specimen completely fractures or reaches the preset simulation duration, terminate the simulation program and conduct a multi-dimensional analysis of the hydraulic fracturing fracture network morphology, seepage capacity evolution, and spatiotemporal distribution of microseismic / acoustic emission events.

[0013] Furthermore, the design of the natural discrete fracture network based on the power-law distribution function in step two specifically includes the following formula: S21. The number and size of fractures are defined by the power-law distribution function: In the formula, For the size range of [ , The number of cracks; denoted as the crack size; 'a' is the scaling index, which is the total crack density determined based on the crack size range. It is a proportionality constant that describes the relative proportion between small and large cracks.

[0014] S22, dimensions in [ , The number of cracks within the range is defined as: S23, the cumulative crack size density is defined as: Based on this formula, the cumulative surface area of ​​the fractures per unit area is calculated, thereby defining the density of the natural discrete fracture network in the model.

[0015] Furthermore, the pore pressure calculation and fluid-structure interaction process in step five specifically includes the following formulas: S51. Formula for calculating fluid velocity based on the cubic law: In the formula, Q is the fluid velocity, which is the volume of fluid flowing into or out of the pore space per unit time; e is the current hydraulic aperture. The pressure difference between adjacent pore spaces is η; the fluid dynamic viscosity is L. c This represents the length of the fluid flow channel.

[0016] S52. Dynamic calculation formula for hydraulic orifice diameter considering effective normal stress: In the formula, e is the current hydraulic aperture (or equivalent hydraulic opening). Residual hydraulic aperture refers to the equivalent hydraulic aperture of a crack under zero normal stress. The initial hydraulic aperture refers to the minimum hydraulic opening that a crack can achieve when the normal stress approaches infinity. This represents the effective normal stress acting on the crack surface.

[0017] S53. Formula for calculating the change in fluid pore pressure: In the formula, It represents the change in fluid pressure in the pore space of the fluid per time interval; Bulk modulus of the fluid; Where is the pore space volume; Q is the fluid velocity; For time intervals; This represents the change in pore space volume caused by mechanical load.

[0018] S54. Apply the calculated fluid pressure change to the particles and contacts around the flow region, and simultaneously execute the formulas in steps S51 to S54 to update the fluid pressure of each flow region hourly. Couple the fluid pressure with the total stress of the rock mass to simulate the fluid transport and rock deformation and failure process at the pore scale.

[0019] Furthermore, the operation method of the virtual microseismic monitoring module in step six specifically includes the following steps: S61. Monitor and record interparticle bonding failure events in real time during each calculation step. Define a single set of bonding failure events as a microseismic / acoustic emission source and record the time and spatial coordinate information of the source simultaneously.

[0020] S62. Calculate the moment tensor of microseismic / acoustic emission events based on moment tensor theory. The calculation formula is as follows: In the formula, This represents the i-th component of the change in contact force. is the j-th component of the distance between the contact point and the event centroid; s is the entire surface surrounding the event.

[0021] The calculated moment tensor is decomposed into isotropic components (ISO), dual-couple components (DC), and compensated linear vector dipole components (CLVD), and the proportion of each component is calculated. An R-value is introduced to characterize the proportion of the isotropic components, calculated using the following formula: In the formula, The trace of the moment tensor; Let i be the i-th partial eigenvalue. ;in The value of R is For a pure implosion event, the R value is... For a pure shearing event, the R value is For pure explosive events, a negative R value indicates an implosion-type seismic source mechanism, while a positive R value indicates an explosive seismic source mechanism.

[0022] S63. Calculate the scalar seismic moment M0 of microseismic / acoustic emission events based on the moment tensor. The calculation formula is as follows: In the formula, m j Let j represent the j-th eigenvalue of the moment tensor.

[0023] Moment magnitude (M) of an event estimated based on the Hanks-Kanamori empirical relation. w The calculation formula is: .

[0024] Furthermore, in step one, an empirical equation is established between the macroscopic mechanical properties of the rock and the microscopic input parameters of the model through indoor rock mechanical tests. The model input parameters are calibrated based on the empirical equation until the deviation between the simulation results and the indoor test results is less than a preset threshold. The indoor rock mechanical tests include uniaxial compressive strength tests, Brazilian splitting tests, and triaxial compression tests. The macroscopic mechanical properties include elastic modulus, Poisson's ratio, uniaxial compressive strength, tensile strength, cohesion, and internal friction angle. The preset threshold is 5%.

[0025] Furthermore, the discrete element fluid-solid coupled hydraulic fracturing model supports extended thermo-fluid-solid-chemical (THMC) multi-field coupling modules, which are used to simulate the influence of temperature field on rock mechanical strength and permeability, as well as the modification of fracture zone pore structure and permeability by chemical dissolution.

[0026] Furthermore, the multi-dimensional analysis in step seven supports the comparison and fitting of the simulated crack network morphology and spatiotemporal distribution characteristics of microseismic / acoustic emission events with indoor acoustic emission test data and on-site microseismic monitoring data, thereby correcting model parameters and improving simulation accuracy.

[0027] Compared with the prior art, the present invention has the following advantages: 1. This invention achieves quantitative design of natural discrete fracture networks based on power-law distribution functions, which can accurately restore the in-situ geometry, size distribution and density characteristics of reservoir fractures. By dynamically adjusting the relative modulus of fractures and the fracture criterion, it can accurately characterize the guiding and controlling role of natural fractures on the initiation, propagation and penetration of hydraulic fractures. The fracture propagation simulation is more in line with engineering practice, and solves the technical problem that existing models cannot accurately quantify the influence of natural fractures.

[0028] 2. This invention constructs a fully coupled algorithm based on Darcy's law and the cubic law, realizing the hourly synchronous update of fluid velocity, hydraulic aperture, and pore pressure. It fully couples the evolution of pore pressure with geostress and particle bonding stress, which can realistically simulate the dynamic interaction process of fluid transport and rock deformation and failure at the pore scale, effectively improving the accuracy of hydraulic fracturing process simulation under multi-field coupling conditions.

[0029] 3. This invention develops a virtual microseismic monitoring module that runs synchronously with hydraulic fracturing simulation. Using particle bonding failure events as the source, it realizes real-time decomposition of the source mechanism and quantitative calculation of scalar seismic moment and moment magnitude based on moment tensor theory. It can directly realize the real-time mapping between the dynamic process of crack propagation and microseismic activity, and solves the defects of existing technologies that can only perform post-event statistics and cannot quantitatively interpret microseismic signals.

[0030] 4. The simulation method of the present invention can simultaneously meet the simulation needs of laboratory rock sample scale and field engineering scale, and supports the benchmarking and fitting with indoor acoustic emission test and field microseismic monitoring data. It can be widely used in multiple fields such as coalbed methane, shale gas, hot dry rock geothermal development, and carbon dioxide geological storage, effectively improving the accuracy of deep reservoir hydraulic fracturing design and construction safety. Attached Figure Description

[0031] Figure 1 This is a schematic diagram of the overall process of the simulation system described in this invention.

[0032] Figure 2 This is a schematic diagram of the flow skeleton network in Embodiment 1 of the present invention.

[0033] Figure 3 This is a diagram of the discrete element numerical model of Embodiment 1 of the present invention.

[0034] Figure 4 This is a distribution diagram of acoustic emission events induced by hydraulic fracturing in Embodiment 1 of the present invention.

[0035] Figure 5 This is a fluid pressure cloud map at the end of hydraulic fracturing in Embodiment 1 of the present invention.

[0036] Figure 6 This is a displacement cloud diagram of the bonded particles after hydraulic fracturing in Embodiment 1 of the present invention.

[0037] Figure 7 This is a graph showing the change of wellbore injection pressure and acoustic emission event propagation distance over time in Embodiment 1 of the present invention. Detailed Implementation

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

[0039] Example 1: This example focuses on the deep coalbed methane extraction scenario. A laboratory-scale numerical simulation of hydraulic fracturing was conducted on coal samples containing natural fractures. Acoustic emission response analysis was performed simultaneously. The simulation results were compared with indoor test data to verify the accuracy of the system of this invention.

[0040] The specific implementation steps of the simulation system in this embodiment are as follows: Step 1: Discrete Element Model Construction and Parameter Calibration: A two-dimensional bonding particle model was established based on the particle flow program PFC2D. The model size was 100mm×100mm, with a particle size ratio of 1.66 and a minimum particle radius of 0.50mm. A total of 15296 bonding particles and 40277 effective particle contacts were generated.

[0041] Indoor uniaxial compressive strength tests and Brazilian splitting tests were conducted on coal samples to obtain the following macroscopic mechanical parameters: elastic modulus 7.89 GPa, Poisson's ratio 0.42, uniaxial compressive strength 5.99 MPa, and tensile strength 0.82 MPa. Based on the indoor test results, an empirical equation relating macroscopic mechanical properties to microscopic input parameters was established. The model's microscopic parameters were repeatedly calibrated until the deviation between the simulation and experimental results was less than 5%. The calibrated model microscopic parameters are shown in Table 1, and the comparison between the calibrated macroscopic mechanical parameters and experimental results is shown in Table 2.

[0042] Table 1: Input values ​​of model parameters after calibration Table 2: Comparison of macroscopic mechanical properties after model calibration with experimental test results Step 2: Design and Model Integration of Natural Discrete Fragment Network (DFN): The DFN model is designed based on the power-law distribution function. The fracture length range is set to 2%~10% of the model size (i.e., 2mm~10mm), the scaling exponent α is set to 2.2, and the scaling constant α is determined according to the preset fracture density. Finally, the DFN density of the model is set to 16.22 fractures / m. 2 With the fracture dip angle at 0°, the quantitative calculation of the number and size distribution of fractures is completed using the above corresponding formulas. The designed DFN is integrated into the calibrated discrete element model, and the relative modulus of fractures is set to 0.3 times that of the intact rock mass. The Mohr-Coulomb fracture criterion is dynamically adjusted to characterize the mechanical properties of natural fractures.

[0043] Step 3: Flow Skeleton Network Construction: Based on the microscopic parameters after model calibration, a flow skeleton network consisting of particle contact pores is constructed to clarify the correspondence between fluid flow channels and pore spaces; for example... Figure 2 As shown, a flow skeleton network is constructed based on the parameter input values ​​obtained from the model calibration process. The resulting model consists of sticky particles and flow channels. Figure 3 As shown, the model consists of 15,296 adhesive particles and 40,277 effective contacts. Black dots represent particles bonded together in the fissures, while the remaining particles are marked in gray. The blue grid depicts the flow channels of individual region structures.

[0044] Step 4: Set the model ground stress boundary conditions and injection boundary conditions: vertical stress 4MPa, horizontal stress 2MPa, to complete the initial ground stress balance of the model; set an injection borehole with a diameter of 4mm at the center of the model, set the injection boundary condition to fixed flow rate injection, injection rate 0.0001m³ / s, calculation time step 0.01s, and fluid related parameters are shown in Table 3.

[0045] Table 3: Parameters used to simulate fluid injection input Step 5: Iterative Injection Simulation and Fluid-Structure Coupling Calculation: Start the iterative injection simulation program according to the parameters in Table 3, and execute the following calculations step by step: S51. Formula for calculating fluid velocity based on the cubic law: In the formula, Q is the fluid velocity, which is the volume of fluid flowing into or out of the pore space per unit time; e is the current hydraulic aperture. The pressure difference between adjacent pore spaces is η; the fluid dynamic viscosity is L. c This represents the length of the fluid flow channel.

[0046] S52. Dynamic calculation formula for hydraulic orifice diameter considering effective normal stress: In the formula, e is the current hydraulic aperture (or equivalent hydraulic opening). Residual hydraulic aperture refers to the equivalent hydraulic aperture of a crack under zero normal stress. The initial hydraulic aperture refers to the minimum hydraulic opening that a crack can achieve when the normal stress approaches infinity. This represents the effective normal stress acting on the crack surface.

[0047] S53. Formula for calculating the change in fluid pore pressure: In the formula, It represents the change in fluid pressure in the pore space of the fluid per time interval; Bulk modulus of the fluid; Where is the pore space volume; Q is the fluid velocity; For time intervals; This represents the change in pore space volume caused by mechanical load.

[0048] S54. Apply the calculated fluid pressure change to the particles and contacts around the flow region, and simultaneously execute the formulas in steps S51 to S54 to update the fluid pressure of each flow region hourly. Couple the fluid pressure with the total stress of the rock mass to simulate the fluid transport and rock deformation and failure process at the pore scale.

[0049] Step Six: Synchronous Virtual Acoustic Emission Monitoring: Start the virtual acoustic emission monitoring module that runs synchronously with the entire injection simulation step, and perform the following operations: S61. Monitor and record interparticle bonding failure events in real time during each calculation step. Define a single set of bonding failure events as a microseismic / acoustic emission source and record the time and spatial coordinate information of the source simultaneously.

[0050] S62. Calculate the moment tensor of microseismic / acoustic emission events based on moment tensor theory. The calculation formula is as follows: In the formula, This represents the i-th component of the change in contact force. is the j-th component of the distance between the contact point and the event centroid; s is the entire surface surrounding the event.

[0051] The calculated moment tensor is decomposed into isotropic components (ISO), dual-couple components (DC), and compensated linear vector dipole components (CLVD), and the proportion of each component is calculated. An R-value is introduced to characterize the proportion of the isotropic components, calculated using the following formula: In the formula, The trace of the moment tensor; Let i be the i-th partial eigenvalue. ;in The value of R is For a pure implosion event, the R value is... For a pure shearing event, the R value is For pure explosive events, a negative R value indicates an implosion-type seismic source mechanism, while a positive R value indicates an explosive seismic source mechanism.

[0052] S63. Calculate the scalar seismic moment M0 of microseismic / acoustic emission events based on the moment tensor. The calculation formula is as follows: In the formula, m j Let j represent the j-th eigenvalue of the moment tensor.

[0053] Moment magnitude (M) of an event estimated based on the Hanks-Kanamori empirical relation. w The calculation formula is: .

[0054] Step 7: Simulation Termination and Result Analysis: When the simulation duration reaches 30 minutes, terminate the simulation program and perform multi-dimensional analysis of the simulation results: S71. Characteristics of acoustic emission event distribution: such as Figure 4 As shown, acoustic emission events were mainly concentrated in the DFN introduction area (30 mm to the left and right of the model center). No acoustic emission events were recorded outside this area. The figure shows that the DFN has a significant impact on the propagation of hydraulic fractures (especially for fractures extending downwards from the model center). Strong acoustic emission events (in...) Figure 4 The parts marked ① to ⑥ are all located at the tip of the natural fracture, which perfectly matches the stress concentration characteristics at the fracture tip. The natural fracture has a significant guiding and controlling effect on the propagation of hydraulic fractures. The hydraulic fractures deviate from the direction of the maximum principal stress and propagate along the natural fracture. After moving away from the fracture tip, they resume extending along the direction of the maximum principal stress.

[0055] S72. Fluid pressure distribution characteristics: such as Figure 5 As shown, the fluid pressure in the region corresponding to the strong acoustic emission event is significantly higher, indicating that the high-pressure fluid in the natural fissure is the core cause of the strong acoustic emission event. Some acoustic emission events far from the fluid flow trajectory are induced by stress disturbance caused by fluid injection.

[0056] S73, Rock mass deformation characteristics: such as Figure 6 As shown, to analyze the impact of hydraulic fracturing on model deformation, displacement contour plots of individual bonded particles in the model were plotted. The particle displacements are symmetrically distributed relative to the hydraulic fractures, with the largest displacements observed near the center of the wellbore, decreasing with distance from the wellbore. Specifically, a distribution zone is formed in the area where pressurized fluid is present. The displacement above this zone exceeds 0.2 mm, while the displacement below it abruptly decreases to less than 0.15 mm.

[0057] S74. Injection pressure and acoustic emission evolution: such as Figure 7 As shown, the system records the propagation distance of the model wellbore injection pressure and acoustic emission events away from the wellbore center. In the first half of the fracturing simulation, strong acoustic emission events occur after a sharp jump in wellbore injection pressure, for example, at 5.5 min and 12 min. In the second half, more strong acoustic emission events are generated, but their occurrence has almost no correlation with fluid pressure changes, because the fluid has moved a sufficiently far distance (25 mm) away from the wellbore center. In this embodiment, the fracture pressure (3.77 MPa) is the local pressure peak corresponding to the first acoustic emission event, marking the initiation of hydraulic fracture. The peak pressure (4.61 MPa) is the maximum pressure during the entire simulation, marking the termination of the indirectly induced fracture propagation caused by fluid injection.

[0058] The simulation results of this embodiment are in high agreement with the acoustic emission test results of hydraulic fracturing of indoor coal samples. The benchmark deviation of the crack propagation morphology and the spatiotemporal distribution of acoustic emission events is less than 10%, which verifies the accuracy of the system of the present invention.

[0059] 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 discrete element method-based simulation method for hydraulic fracturing and microseismic / acoustic emission, characterized in that, Includes the following steps: Step 1: Based on the cohesive particle model of the particle flow program, establish a discrete element fluid-solid coupled hydraulic fracturing model and calibrate the model input parameters; Step 2: Design a natural discrete fracture network based on the power-law distribution function, and integrate the network into the calibrated discrete element fluid-structure coupled hydraulic fracturing model to characterize the in-situ geometry of rock fractures; Step 3: Based on the microscopic parameters after model calibration, construct a flow skeleton network composed of particle contact pores to clarify the correspondence between fluid flow channels and pore space; Step 4: Set the model's ground stress boundary conditions and injection boundary conditions, start the iterative injection simulation program, and inject the fluid into the pre-specified central borehole of the model according to the preset fixed flow rate; Step 5: Calculate the fluid pore pressure in the flow region based on Darcy's law and the cubic law, couple the pore pressure with the ground stress and particle bonding stress to calculate the effective stress, and update the pore space volume and fluid flow parameters simultaneously to simulate the dynamic interaction process of fluid transport and rock mechanical response at the pore scale. Step 6: Configure a virtual microseismic monitoring module that runs synchronously with the hydraulic fracturing simulation program throughout the entire process. This module monitors and records interparticle bonding failure events in real time during each calculation step, and simultaneously completes the source location, source mechanism analysis, and magnitude calculation of microseismic / acoustic emission events. Step 7: When the model specimen completely fractures or reaches the preset simulation duration, terminate the simulation program and conduct a multi-dimensional analysis of the hydraulic fracturing fracture network morphology, seepage capacity evolution, and spatiotemporal distribution of microseismic / acoustic emission events.

2. The method according to claim 1, characterized in that, Step two involves designing a natural discrete fracture network based on a power-law distribution function, specifically including the following formula: S21. The number and size of fractures are defined by the power-law distribution function: In the formula, For the size range of [ , The number of cracks; denoted as the crack size; 'a' is the scaling index, which is the total crack density determined based on the crack size range. It is a proportionality constant that describes the relative proportion between small and large cracks; S22, dimensions in [ , The number of cracks within the range is defined as: S23, the cumulative crack size density is defined as: Based on this formula, the cumulative surface area of ​​the fractures per unit area is calculated, thereby defining the density of the natural discrete fracture network in the model.

3. The method according to claim 1, characterized in that, The pore pressure calculation and fluid-structure interaction process in step five specifically includes the following formulas: S51. Formula for calculating fluid velocity based on the cubic law: In the formula, Q is the fluid velocity, which is the volume of fluid flowing into or out of the pore space per unit time; e is the current hydraulic aperture. The pressure difference between adjacent pore spaces is η; the fluid dynamic viscosity is L. c The length of the fluid flow channel; S52. Dynamic calculation formula for hydraulic orifice diameter considering effective normal stress: In the formula, e is the current hydraulic aperture; Residual hydraulic aperture refers to the equivalent hydraulic aperture of a crack under zero normal stress. The initial hydraulic aperture refers to the minimum hydraulic opening that a crack can achieve when the normal stress approaches infinity. This refers to the effective normal stress acting on the crack surface; S53. Formula for calculating the change in fluid pore pressure: In the formula, It represents the change in fluid pressure in the pore space of the fluid per time interval; Bulk modulus of the fluid; Where is the pore space volume; Q is the fluid velocity; For time intervals; This refers to the change in pore space volume caused by mechanical load. S54. Apply the calculated fluid pressure change to the particles and contacts around the flow region, and simultaneously execute the formulas in steps S51 to S54 to update the fluid pressure of each flow region hourly. Couple the fluid pressure with the total stress of the rock mass to simulate the fluid transport and rock deformation and failure process at the pore scale.

4. The method according to claim 1, characterized in that, The operation method of the virtual microseismic monitoring module in step six specifically includes the following steps: S61. Monitor and record interparticle bonding failure events in real time during each calculation step, define a single set of bonding failure events as a microseismic / acoustic emission source, and synchronously record the time and spatial coordinate information of the source. S62. Calculate the moment tensor of microseismic / acoustic emission events based on moment tensor theory. The calculation formula is as follows: In the formula, This represents the i-th component of the change in contact force. is the j-th component of the distance between the contact point and the event centroid; s is the entire surface surrounding the event; The calculated moment tensor is decomposed into isotropic components, bicouple components, and compensated linear vector dipole components, and the proportion of each component is calculated. An R-value is introduced to characterize the proportion of the isotropic components, and the calculation formula is as follows: In the formula, The trace of the moment tensor; Let i be the i-th partial eigenvalue. ;in The value of R is For a pure implosion event, the R value is... For a pure shearing event, the R value is For pure explosive events; when the R value is negative, the source mechanism is determined to be implosion type; when the R value is positive, the source mechanism is determined to be explosive type. S63. Calculate the scalar seismic moment M0 of microseismic / acoustic emission events based on the moment tensor. The calculation formula is as follows: In the formula, m j This represents the j-th eigenvalue of the moment tensor; Moment magnitude (M) of an event estimated based on the Hanks-Kanamori empirical relation. w The calculation formula is: 。 5. The method according to claim 1, characterized in that, In step one, an empirical equation is established between the macroscopic mechanical properties of rock and the microscopic input parameters of the model through indoor rock mechanical tests. The model input parameters are calibrated based on the empirical equation until the deviation between the simulation results and the indoor test results is less than a preset threshold. The indoor rock mechanical tests include uniaxial compressive strength test, Brazilian splitting test, and triaxial compression test. The macroscopic mechanical properties include elastic modulus, Poisson's ratio, uniaxial compressive strength, tensile strength, cohesion, and internal friction angle.

6. The method according to claim 1, characterized in that, The discrete-element fluid-solid coupled hydraulic fracturing model supports an extended thermo-fluid-solid-chemical multi-field coupling module. This module is used to simulate the effects of temperature field on rock mechanical strength and permeability, as well as the modification of fracture zone pore structure and permeability by chemical dissolution.

7. The method according to claim 1, characterized in that, The multi-dimensional analysis in step seven supports the comparison and fitting of the simulated crack network morphology and spatiotemporal distribution characteristics of microseismic / acoustic emission events with indoor acoustic emission test data and on-site microseismic monitoring data, thereby correcting model parameters and improving simulation accuracy.

Images