A Quantitative Method for Enhancing Permeability of Rock Mass Hydraulic Fracturing Based on DEM

Abstract

This invention discloses a method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM (Discrete Image Model), belonging to the field of simulation analysis technology. The invention includes the following steps: obtaining the layered physical and mechanical parameters of the coal and rock mass using multi-source data fusion; constructing a layered discrete element particle aggregate-fracture dual-network model to generate a bpm model of the coal and rock mass sample; calculating the actual apparent area of ​​the fluid domain using the bpm model of the coal and rock mass sample; assigning seepage hydraulic parameters to the bpm model of the coal and rock mass sample, calculating the initial permeability, simulating the actual hydraulic fracturing process, and calculating the permeability of the bpm model of the coal and rock mass sample under constant pressure steady-state conditions. This invention, through the constructed dual-network model, can clearly present the particle distribution and fracture morphology of each layer. Combined with the coupling module, it restores the original fracture state of the coal and rock mass, and incorporates the influence of temperature on geostress and fluid viscosity, making the simulation environment highly matched with actual geological conditions, thereby improving the permeability enhancement effect of hydraulic fracturing.

Description

Technical Field

[0001] This invention belongs to the field of simulation analysis technology, and in particular relates to a method for quantifying the permeability enhancement effect of hydraulic fracturing of rock mass based on DEM, which solves the problem that traditional methods cannot accurately quantify the evolution of fracture networks. Background Technology

[0002] In the energy extraction sector, gas drainage from low-permeability coal and rock masses has always been an extremely challenging task. For coal and rock masses characterized by low permeability and high density, gas escapes effectively, significantly impacting the efficiency and safety of gas drainage. According to relevant statistics, most coal mining areas in my country have poor coal seam permeability and high gas content, resulting in a large proportion of mine accidents and causing serious casualties and economic losses. With the continuous increase in coal seam mining depth and the rise in outburst-prone coal seams, the difficulty of gas control further increases, seriously threatening the production safety of coal mines.

[0003] Hydraulic fracturing technology, as an effective method for depressurizing and enhancing the permeability of coal and rock masses, plays a crucial role in gas extraction. This technology injects high-pressure fluid into the coal and rock mass, causing fractures to form and propagate, thereby increasing the permeability of the coal and rock mass, expanding the effective range of the borehole, and improving gas extraction efficiency. Currently, hydraulic fracturing technology is widely used in oil and gas extraction and is gradually being promoted in coal mine gas extraction.

[0004] However, existing research on hydraulic fracturing mainly focuses on the development and propagation of hydraulic fractures, as well as the modes and areas of rock failure, aiming to explore the mechanism of fracture propagation in coal and rock masses under hydraulic fracturing. While these studies provide a theoretical basis for the application of hydraulic fracturing technology, for gas drainage operations, it is more crucial to evaluate the permeability enhancement effect of hydraulic fracturing on coal and rock masses, and the influence of different factors on this effect. Existing research has significant shortcomings in this regard and cannot directly provide effective guidance for production operations.

[0005] In the prior art, such as the application number 202310292717.7, the quantitative method for the permeability enhancement effect of hydraulic fracturing of coal and rock mass based on DEM still has the following shortcomings, such as: (1) only by testing coal and rock mass samples in the laboratory to obtain conventional physical and mechanical parameters such as Young's modulus and tensile strength, the sample is single and lacks spatial distribution consideration; (2) by giving parallel bonding contact force between particles to generate BPM model, the original fracture state of coal and rock mass is not fully restored; (3) the contact force parameter is calibrated by trial and error, the seepage hydraulic parameters are fixed and the actual working condition changes are not considered; (4) the model is simulated by injecting fluid into the model with a fixed flow rate, which can easily destroy the model; (5) the permeability enhancement effect is quantified by comparing the overall permeability before and after fracturing.

[0006] To fill this gap, this invention proposes a quantitative method for enhancing the permeability of hydraulic fracturing in rock masses based on the Discrete Element Method (DEM). The DEM has excellent applicability in the study of mechanical problems involving large deformations and discontinuous media, and can effectively simulate the mechanical behavior of discontinuous media such as coal and rock masses during hydraulic fracturing. This method allows for more accurate quantification of the permeability enhancement effect of hydraulic fracturing, in-depth analysis of the influence of different factors on the permeability enhancement effect, and provides scientific and reliable guidance for gas drainage operations, possessing significant engineering practical significance and economic value. Summary of the Invention

[0007] The purpose of this invention is to provide a quantitative method for the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM. Through innovative designs such as multi-source parameter acquisition, dual-network model construction, multi-field coupling simulation, multi-condition testing, and layered-overall collaborative evaluation, this method solves the problems of insufficient analysis of the permeability enhancement effect of hydraulic fracturing caused by the one-sidedness of existing quantification, inaccurate parameters, and model simplification.

[0008] To solve the above-mentioned technical problems, the present invention is achieved through the following technical solution:

[0009] This invention provides a method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM, comprising the following steps:

[0010] Step S1: Use a multi-source data fusion method of "on-site borehole sampling + laboratory precision testing + ground-penetrating radar inversion" to obtain the layered physical and mechanical parameters of coal and rock mass;

[0011] Step S2: Construct a hierarchical discrete element particle assembly-fracture dual-network model;

[0012] Step S3: Apply parallel bonding contact force to the particle aggregate model. The parallel bonding model enhances the connection strength between particles by introducing virtual bonds, generating a coal and rock mass sample bpm model.

[0013] Step S4: Define the fluid conduit as the contact part between particles in the particle aggregate, and the fluid domain as the polygonal closed area enclosed by the center points of adjacent and contacting particles. Calculate the actual apparent area of ​​the fluid domain using the bpm model of the coal and rock mass sample.

[0014] Step S5: Assign the seepage hydraulic parameters to the bpm model of the coal and rock mass sample and calculate the initial permeability;

[0015] Step S6: Remove the confining pressure, inlet pressure, and outlet pressure applied to the coal and rock mass sample bpm model in step S5, restoring the model to its initial pressureless state. Apply a fluid with a certain flow rate and pressure value to the coal and rock mass sample bpm model from the injection hole drilled in step S5 to simulate the actual hydraulic fracturing process.

[0016] Step S7: Remove the calculation results assigned to the bpm model of the coal and rock mass sample in Step S6, i.e., clear the temporary data and state changes generated during the hydraulic fracturing simulation, and restore the model to the initial bpm model state; repeat Step S6, assign a fixed confining pressure to the model again, drill the injection hole, and apply the inlet pressure. and export pressure After waiting for the average pore pressure and total flow rate of the model to reach a steady state, the permeability of the coal and rock mass sample under constant pressure steady-state condition (bpm model) is calculated. .

[0017] As a preferred technical solution, the specific process for obtaining the physical and mechanical parameters of coal and rock mass stratification in step S1 is as follows:

[0018] Step S11: Select 3-5 representative boreholes on site and take samples at 2m intervals to ensure coverage of different lithological sections of the coal and rock mass;

[0019] Step S12: The laboratory uses equipment such as a servo press and a direct shear tester to test the Young's modulus, tensile strength, shear strength, friction angle, friction coefficient and Poisson's ratio of each layer of specimens. Each parameter is tested three times and the average value is taken.

[0020] Step S13: Use ground-penetrating radar to scan the coal and rock mass within a 10m radius around the borehole, invert parameters such as layered porosity and fracture development density, cross-validate with laboratory data, correct parameter deviations, and finally establish a coal and rock mass layered parameter database to provide accurate layered parameter input for subsequent model construction.

[0021] As a preferred technical solution, the specific process for constructing the hierarchical discrete element particle assembly-fracture dual network model in step S2 is as follows:

[0022] Step S21: Based on the layering parameters obtained in step S1, set the domain in the PFC software according to the actual layering thickness of the coal and rock mass, and generate an independent wall space for each layer.

[0023] Step S22: For different layers, based on their porosity and particle size distribution characteristics, generate random rigid particles with differentiated diameters (particle size ratio is adjusted according to the lithology of the layer, 1:3 for sandstone layers and 1:2 for coal seam layers), to form the basic network of particle aggregates.

[0024] Step S23: The fracture development pattern is obtained by inverting the reference ground-penetrating radar. Fractures with different orientations, dip angles and lengths are pre-set in each layer. The fracture width is set according to the particle diameter at a ratio of 1:5 to 1:8 to form a particle aggregate-fracture dual network model.

[0025] As a preferred technical solution, the specific process for constructing the basic network of particle aggregates in step S22 is as follows:

[0026] Step S221: Extract key parameters of each layer from the database in step S1, determine the particle parameters of sandstone, coal seam and mudstone layers, set the particle contact model to linear contact by default, and unify the particle density to 2500 kg / m³.

[0027] Step S222: Activate the sandstone layer wall group, input the particle size range and calculate the number of particles through porosity, and generate particles in the sandstone layer domain;

[0028] Step S223: Repeat step S222 to generate particles in the coal seam and mudstone layer respectively, and determine the volume of particles in each layer;

[0029] Step S224: Gravity compaction is performed on each layer of generated particles until the stress of the particle system is stable; after compaction, the particle overlap rate is checked; if the particle overlap rate is greater than 3%, the particle positions are readjusted to ensure that the particle aggregate structure is stable and conforms to the actual dense state of coal and rock mass.

[0030] Step S225: Extract the fracture characteristics of each layer from the ground-penetrating radar inversion data in step S1, and determine the fracture parameters;

[0031] Step S226: Generate planar cracks based on the parameters of each layer of cracks;

[0032] Step S227: Observe the spatial distribution of the particle aggregate and the cracks, calculate the porosity of the model, and if the deviation from the target value is greater than 3%, adjust the crack density or width until the porosity meets the requirements, and finally form a particle geometry-crack dual network model.

[0033] As a preferred technical solution, the specific process for generating the bpm model of the coal and rock mass sample in step S3 is as follows:

[0034] Step S31: Perform simulation tests and record the stress-strain curve, peak strength, particle displacement field at failure and crack propagation trajectory in real time for each sub-Domain; set the loading rate to 0.001 m / s for each sub-Domain and the loading direction to be vertical.

[0035] Step S32: Simulate and generate simulated stress-strain curves for each layer, and extract simulated failure modes;

[0036] Step S33: Compare the simulation results with the uniaxial compression test results of the same layered specimen in the laboratory from Step S1, and calculate the deviation of key indicators. The formula for calculating the deviation of indicators is:

[0037] ;

[0038] Step S34: Set the parameter adjustment rules and re-copy the single-axis compression simulation according to the adjusted parameters;

[0039] The parameter adjustment rules are as follows:

[0040] When the peak strength deviation is >5%, the tensile strength / shear strength is corrected by 1.1 times the deviation rate; when the elastic modulus deviation is >5%, the Young's modulus is corrected by 1.05 times the deviation rate; if the failure modes do not match, the friction angle is adjusted.

[0041] Step S35: Recalculate the deviation. If the deviation of all indicators... If the calibration is successful, the calibration is complete; if not, the deviation calculation is repeated until the standard is met.

[0042] Step S36: Finally, accurate BPM models for each layer are formed.

[0043] As a preferred technical solution, in step S4, the actual apparent area of ​​the fluid domain is calculated using the bpm model of the coal and rock mass sample as follows:

[0044] Step S41: Based on the particle contact detection function of the DEM model, filter out particle pairs with valid contact conditions, such as contact type filtering, contact force threshold determination and output contact pair information.

[0045] Step S42: For the effective disconnection pairs after screening, assign physical properties to the fluid pipes, establish a contact-pipe mapping relationship, and assign a unique identifier to each fluid pipe to avoid confusion between pipe and particle contact in subsequent simulations;

[0046] Step S43: Randomly select a central particle C from the model, and extract all particles that have formed effective contact with the central particle C, denoted as the surrounding particle group D;

[0047] Step S44: Detect whether there are particles in contact with each other in the surrounding particle group D, and ensure that the central particle C and the surrounding particle group D can jointly form a closed area; if they cannot form a closed area, replace the central particle C.

[0048] Step S45: Determine the vertex coordinates of the polygon or polyhedron of the fluid domain based on the particle center coordinates and radius;

[0049] Step S46: Determine the flow volume or area based on the vertex coordinates of the polygon or polyhedron of the fluid domain, and mark the fluid pipes connected to each fluid domain;

[0050] Step S47: Extract particle coordinates and boundary correction parameters from the DEM model, and calculate the actual apparent area.

[0051] As a preferred technical solution, in step S45, when determining the vertex coordinates of the polygon or polyhedron of the fluid domain, firstly, the common external tangents of the central particle C and the surrounding particle D1 are drawn to obtain the intersection point P1 of the two common external tangents. Then, the intersection point P2 of the common external tangents of the central particle C and the surrounding particle D2 and the intersection point P3 of the common external tangents of the central particle C and the surrounding particle D3 are calculated in sequence. Next, the intersection point P4 of the common tangents of the surrounding particles D1 and D2 is calculated. Similarly, the intersection point P5 of the common tangents of the surrounding particles D2 and D3 and the intersection point P6 of the common tangents of the surrounding particles D3 and D1 are calculated. The intersection points P1-P6 are arranged in a clockwise or counterclockwise order to form the vertex coordinates of the closed polygon, ensuring that the vertices do not intersect after being connected in sequence.

[0052] As a preferred technical solution, the specific process of assigning seepage hydraulic parameters to the bpm model of the coal and rock mass sample in step S5 is as follows:

[0053] Step S51: Obtain the initial opening of the fluid conduit, the half-open compressive force, the apparent volume of the fluid domain, the residual opening magnification factor of the fracture, the bulk modulus of the fluid, and the fluid viscosity from the bpm model of the coal and rock mass sample;

[0054] Step S52: Perform correlation verification and correction on the parameters;

[0055] Step S53: Apply a fixed inlet pressure p1 to one side of the coal and rock mass sample bpm model and a fixed outlet pressure p2 to the other side, where p1>p2, to create a pressure difference and drive the fluid to flow in the model; wait for the average pore pressure and total flow rate of the coal and rock mass sample bpm model to reach a steady state. The criterion is that the rate of change of the average pore pressure and total flow rate is less than the set threshold within several consecutive calculation steps.

[0056] Step S54: When a steady state is reached, calculate the permeability k0 of the bpm model of the coal and rock mass sample according to Darcy's law.

[0057] As a preferred technical solution, the specific process of simulating the actual hydraulic fracturing process in step S6 is as follows:

[0058] Step S61: Before simulation, perform model pressureless reset and initial state calibration. Based on the injection hole drilled in step S5, redefine the boundary conditions for fracturing injection.

[0059] Step S62: Based on the on-site fracturing construction data, determine the injection parameters for fracturing and set dynamic control parameters;

[0060] Step S63: Use DEM software to dynamically simulate the fracturing process, adjust the injection state in real time, and ensure that the fracture propagation pattern is consistent with the field.

[0061] Step S64: Implement multi-dimensional data recording to ensure coverage of particle state, crack state, and fluid state.

[0062] As a preferred technical solution, in step S7, the permeability k0 and permeability k1 are compared to quantify the permeability enhancement effect of the bpm model hydraulic fracturing on the coal and rock mass sample; the permeability enhancement effect is measured by the rate of change of permeability, and the calculation formula is: Permeability enhancement = The higher the permeability enhancement rate, the more significant the permeability enhancement effect of hydraulic fracturing on coal and rock mass. By analyzing the influence of different factors (such as fracturing fluid pressure, flow rate, physical and mechanical parameters of coal and rock mass, etc.) on the permeability enhancement rate, we can gain a deeper understanding of the intrinsic mechanism of hydraulic fracturing permeability enhancement and provide scientific guidance for actual production operations.

[0063] The present invention has the following beneficial effects:

[0064] (1) This invention breaks through the limitations of traditional single-sample testing by multi-source fusion, not only covering different lithological sections of coal and rock mass, but also correcting parameter deviations through cross-validation, reducing the error of physical and mechanical parameters to within 5%, effectively solving the simulation distortion problem caused by insufficient parameter representativeness in traditional methods;

[0065] (2) The dual-network model constructed in this invention can clearly present the particle distribution and fracture morphology of each layer. When the contact force is calibrated in S3, the macroscopic mechanical response of the model (such as the uniaxial compressive strength of the coal seam of 15MPa, which deviates from the laboratory test value of 14.8MPa by only 1.3%) is highly consistent with the real coal and rock mass, laying a reliable foundation for subsequent multi-field coupling simulation and quantification of permeability enhancement effect;

[0066] (3) This invention classifies six core parameters into “pipeline attribute parameters, fluid attribute parameters, and crack response parameters”, and through four steps of “parameter value selection, model assignment, dynamic correlation, and verification correction”, it achieves accurate matching between parameters and the actual seepage characteristics of coal and rock mass, laying the foundation for subsequent hydraulic fracturing simulation and permeability calculation;

[0067] (4) Each parameter of the present invention is based on laboratory test or field monitoring data, avoiding the randomness of traditional empirical assignment and ensuring the credibility of the simulation results. The parameters are correlated through Darcy's law, contact mechanics model, etc. For example, the initial opening of the pipeline and the viscosity jointly determine the initial seepage velocity, and the residual opening amplification factor of the crack is linked with the pressure unloading triggering condition to restore the dynamic changes of the actual seepage.

[0068] (5) This invention avoids the problems of "sudden pressure rise leading to model collapse" or "insufficient flow leading to crack inability to expand" in traditional simulations by "pressure-flow dual control" and "dynamic switching threshold", and improves the simulation success rate to over 90%. At the same time, multi-dimensional and high-frequency data recording covers the entire system of "particle-crack-fluid", and the simulation results are highly matched with the field conditions, which can be directly used to optimize field fracturing parameters and reduce field test costs and risks.

[0069] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0070] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0071] Figure 1 This is a flowchart of a method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM, according to the present invention. Detailed Implementation

[0072] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0073] Furthermore, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0074] To make the purpose, technical solution, and advantages of this application clearer, the following description is provided in conjunction with the appendix. Figure 1 The present application will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of the application.

[0075] Please see Figure 1 As shown, this invention is a method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM, comprising the following steps:

[0076] Step S1: Use a multi-source data fusion method of "on-site borehole sampling + laboratory precision testing + ground-penetrating radar inversion" to obtain the layered physical and mechanical parameters of coal and rock mass;

[0077] Step S2: Construct a hierarchical discrete element particle assembly-fracture dual-network model;

[0078] Step S3: Apply parallel bonding contact force to the particle aggregate model. The parallel bonding model enhances the connection strength between particles by introducing virtual bonds, generating a coal and rock mass sample bpm model.

[0079] Step S4: Define the fluid conduit as the contact part between particles in the particle aggregate, and the fluid domain as the polygonal closed area enclosed by the center points of adjacent and contacting particles. Calculate the actual apparent area of ​​the fluid domain using the bpm model of the coal and rock mass sample.

[0080] Step S5: Assign the seepage hydraulic parameters to the bpm model of the coal and rock mass sample and calculate the initial permeability;

[0081] Step S6: Remove the confining pressure, inlet pressure, and outlet pressure applied to the coal and rock mass sample bpm model in step S5, restoring the model to its initial pressureless state. Apply a fluid with a certain flow rate and pressure value to the coal and rock mass sample bpm model from the injection hole drilled in step S5 to simulate the actual hydraulic fracturing process.

[0082] Step S7: Remove the calculation results assigned to the bpm model of the coal and rock mass sample in Step S6, i.e., clear the temporary data and state changes generated during the hydraulic fracturing simulation, and restore the model to the initial bpm model state; repeat Step S6, assign a fixed confining pressure to the model again, drill the injection hole, and apply the inlet pressure. and export pressure After waiting for the average pore pressure and total flow rate of the model to reach a steady state, the permeability of the coal and rock mass sample under constant pressure steady-state condition (bpm model) is calculated. .

[0083] In step S1, the specific process for obtaining the physical and mechanical parameters of coal and rock mass stratification is as follows:

[0084] Step S11: Select 3-5 representative boreholes on site and take samples at 2m intervals to ensure coverage of different lithological sections of the coal and rock mass;

[0085] Step S12: The laboratory uses equipment such as a servo press and a direct shear tester to test the Young's modulus, tensile strength, shear strength, friction angle, friction coefficient and Poisson's ratio of each layer of specimens. Each parameter is tested three times and the average value is taken.

[0086] In practice, before construction, based on the geological survey report of the work area (such as coal seam depth, lithological distribution, fault location, etc.), interfering areas such as fault fracture zones and abnormal water inflow areas are eliminated. On three different strikes (e.g., strikes of 15°, 45°, and 75°) of the target coal seam for gas extraction, one representative borehole point is selected for each (3-5 in total; this example uses 3). This ensures that the points cover complete lithological sections including the overlying strata (sandstone, mudstone), the target coal seam, and the underlying strata (limestone). The construction equipment used is an XY-44 core drilling rig (drilling depth ≥1000m, meeting the sampling requirements of medium-deep wells), equipped with a diamond drill bit (91mm diameter to reduce core breakage), double-layer core tubing (to protect core integrity), and tools such as a core box (with layer labels), a geological compass (to measure borehole inclination and strike), and a measuring tape (to measure core length).

[0087] During drilling and sampling, drill at the designed borehole inclination angle (e.g., 75°, consistent with the coal seam strike). Record the drilling pressure (controlled at 15-20 kN), rotation speed (200-300 r / min), and flushing fluid flow rate (50-60 L / min) every 50 m during drilling to avoid core breakage due to parameter fluctuations. When drilling reaches the predetermined lithological interface (e.g., from sandstone to coal seam, judge by changes in drilling speed: approximately 0.5 m / h in sandstone section, and 1.2 m / h in coal seam section), reduce the drilling pressure to 10-12 kN and drill slowly to ensure core integrity. When sampling by layers, layers are drilled at 2m intervals. Drilling is stopped every 2m, the core tube is pulled out, and the core is taken out. The strike and dip of the core are measured with a geological compass, and the length of the core is measured with a tape measure (the core recovery rate should be ≥85%. If it is less than 85%, additional drilling and sampling are required in that layer). The core is placed in the core box in the order of "top to bottom", and a label is attached (indicating the borehole number, layer depth: such as ZK1-5-7m, representing layer 5-7m of borehole No. 1, lithology: sandstone, and sampling date). At the same time, core photos are taken (including the scale, recording the appearance of the core and the development of fractures).

[0088] Three complete core samples were selected from each layer (size requirements: length ≥ 100 mm, diameter ≥ 80 mm, no obvious cracks or holes), wrapped in plastic wrap (to prevent moisture loss), placed in a sealed sample bag, and labeled "Test Sample - Drill Hole Number - Layer Depth"; the remaining core samples were used as backup samples, sealed and stored in the same way for subsequent supplementary testing or verification.

[0089] The specific implementation is as follows:

[0090] The test samples obtained on-site were sent to the laboratory, where a Q-2 core cutter (accuracy ±0.1mm) was used to cut the cores into standard specimens: uniaxial compression specimens (50mm in diameter, 100mm in height, height-to-diameter ratio 2:1), Brazilian splitting specimens (50mm in diameter, 25mm in height), and straight shear specimens (50mm×50mm×50mm). Three parallel samples were prepared for each type of specimen in each layer.

[0091] Sample pretreatment: Place the prepared samples in a vacuum drying oven (temperature 60℃, vacuum degree -0.09MPa) and dry for 24 hours to remove internal moisture. After drying, use a vernier caliper (accuracy 0.02mm) to measure the diameter, height, side length and other dimensions of each sample and record the data (e.g., uniaxial compression sample dimensions: diameter 50.02mm, height 100.05mm), ensuring that the dimensional deviation is ≤0.5%.

[0092] The specific test content is as follows:

[0093] 1. Young's modulus and Poisson's ratio test (uniaxial compression test):

[0094] Let a value be set to achieve this. The negative sign indicates that the radial strain and axial strain are in opposite directions. The average value of the ratio of radial strain to axial strain in the elastic stage is taken.

[0095] 2. Tensile strength test (Brazilian splitting test):

[0096] A servo press, similar to that used for uniaxial compression testing, is equipped with a Brazilian splitting clamp (steel, with cardboard pads 2mm thick). Place the Brazilian splitting specimen in the center of the clamp, ensuring the pads are aligned with the specimen axis; press... The loading rate was adjusted, and load and displacement data were recorded until the specimen fractured along the diameter direction. The maximum load at failure was recorded. ), and according to the formula Calculate the tensile strength (where d is the specimen diameter and h is the specimen height).

[0097] 3. Shear strength, friction angle, and friction coefficient tests (direct shear test):

[0098] A ZJ-type rock direct shear apparatus (maximum normal load 500kN, maximum shear load 300kN, displacement accuracy 0.001mm) was used. The direct shear specimen was placed in the shear box, and normal pressures (5MPa, 10MPa, and 15MPa were applied to cover the field stress range) were applied. Each normal pressure was maintained for 30 minutes to allow the specimen stress to stabilize. Then, a shear load was applied at a shear rate of 0.05mm / min, and the shear load and shear displacement data were recorded until the shear load reached its peak value (or the shear displacement reached 10% of the specimen side length), at which point loading was stopped.

[0099] According to different normal pressures ( The corresponding peak shear load () ), plot the τ-σ curve, the slope of the curve is the coefficient of friction ( The intercept is the cohesion ( ); Friction angle

[0100] After testing three parallel samples for each parameter (Young's modulus, Poisson's ratio, tensile strength, shear strength, friction angle, and friction coefficient), calculate the average value (e.g., if the Young's moduli of three uniaxial compression samples in a layer are 25.2 GPa, 24.8 GPa, and 25.0 GPa, the average value is 25.0 GPa). If the test data of a certain parallel sample deviates from the other two by more than 10% (e.g., if the Young's modulus of a certain sample is 22.0 GPa, and the deviation from the average value is 12%), then discard that data and calculate the average value using the remaining two data. If the deviation of all three data exceeds 10%, then prepare a new sample for retesting.

[0101] Step S13: Use ground-penetrating radar to scan the coal and rock mass within a 10m radius around the borehole, invert parameters such as layered porosity and fracture development density, cross-validate with laboratory data, correct parameter deviations, and finally establish a coal and rock mass layered parameter database to provide accurate layered parameter input for subsequent model construction.

[0102] Centered on each borehole, within a 10m radius around the borehole, three survey lines are laid out in two directions: parallel to the coal seam strike and perpendicular to the coal seam strike (a total of six survey lines). The survey lines are spaced 2m apart and are 20m long (covering 10m on each side of the borehole). Coordinate stakes are placed at the starting point, ending point, and turning point of the survey lines, and the coordinates are recorded using a GPS positioning instrument (accuracy ±0.5m) (e.g., starting point coordinates: X=3852100.5m, Y=521030.2m, Z=+1200.3m).

[0103] The scanning parameters were set as follows: sampling rate 2048 points / scan, time window length 100ns (corresponding to a detection depth of about 10m, calculated based on the electromagnetic wave propagation speed in coal and rock of 15cm / ns), and superposition times 64 times (to reduce noise interference); during scanning, the radar antenna moved at a constant speed (0.5m / s) along the survey line, and radar profile data was recorded every 0.1m of movement. At the same time, the camera was used to photograph the surrounding terrain and lithological outcrops to assist in data interpretation.

[0104] The radar data collected on-site was imported into ReflexW software for preprocessing: ① DC drift removal (using linear trend removal method); ② Filtering (using bandpass filtering, frequency range 50-150MHz, to remove high-frequency noise and low-frequency interference); ③ Gain adjustment (using exponential gain to compensate for electromagnetic wave attenuation with depth); ④ Inter-channel equalization (to make the signal amplitude of each channel consistent).

[0105] Based on the preprocessed radar profile, combined with the field borehole lithology records (e.g., strong radar wave amplitude and continuous phase axis in sandstone sections, and weak radar wave amplitude and discontinuous phase axis in coal seam sections), the lithological stratification interfaces corresponding to each survey line are identified (e.g., the sandstone-coal seam interface appears as a strong reflection interface on the radar profile, with a depth of about 6m), and the thickness of each stratum is determined (e.g., the coal seam stratification thickness is 2.5m).

[0106] Porosity Inversion: Based on the Correlation between Radar Wave Velocity and Porosity (Electromagnetic Wave Velocity in Coal and Rock) Where v0 is the electromagnetic wave velocity of the rock skeleton and φ is the porosity), the electromagnetic wave velocity of each layer of the rock skeleton is first obtained by drilling core tests (e.g., sandstone skeleton v0=20cm / ns, coal seam skeleton v0=20cm / ns). Then, extract the electromagnetic wave propagation velocity of each layer from the radar profile (calculated during travel along the phase axis). (where d is the layer depth and t is the travel time); finally, follow the formula... Calculate the porosity of layered structures (e.g., electromagnetic wave velocity in sandstone layers). Porosity (i.e., 43.75%).

[0107] Crack development density inversion: In the radar profile, cracks appear as discontinuities, distortions, or localized strong reflections of the reflection phase axis; using the "window counting method", each layer is divided into 1m×1m windows, and the number of crack reflections in each window is counted (each crack corresponds to 1 reflection). If the thickness of a coal seam is 2.5m and the number of reflections of a fracture within a certain window is 5, then the fracture development density = 5 / (1×1×2.5) = 2 fractures / m³.

[0108] Specifically, in step S13, the laboratory data cross-validation process is as follows:

[0109] The mechanical parameters and porosity were verified as follows: There is a negative correlation between porosity and Young's modulus within the same layer (the higher the porosity, the lower the Young's modulus). For example, the Young's modulus of a sandstone layer was 28 GPa in laboratory tests, and the porosity obtained from ground-penetrating radar was 35%. If the Young's modulus of an adjacent sandstone layer was 22 GPa (lower), but the obtained porosity was 30% (lower), then the radar data preprocessing process for that layer needed to be re-examined (e.g., whether the gain adjustment was reasonable), or an additional porosity laboratory test (using the helium replacement method, with an accuracy of ±0.1%) was performed to correct the inversion deviation (e.g., if the additional measured porosity was 38%, which is negatively correlated with the Young's modulus, the obtained porosity was corrected to 38%).

[0110] The relationship between fracture development density and tensile strength is verified as follows: the higher the fracture development density, the lower the tensile strength (fractures are stress concentration sources, which can easily lead to premature sample failure). For example, the tensile strength of a certain coal seam is 1.2 MPa in laboratory tests, and the fracture development density retrieved by radar is 8 fractures / m³. If the tensile strength of another coal seam is 0.8 MPa (lower), but the retrieved fracture development density is 5 fractures / m³ (lower), then it is necessary to rescan the radar survey lines corresponding to that layer (to check whether any fracture reflection signals are missed), or to conduct a direct shear test (to observe the number of fractures when the sample fails) to correct the fracture development density (e.g., if the retest finds that the sample has 7 fractures when it fails, the retrieved density is corrected to 7 fractures / m³).

[0111] When constructing the database, Excel or Access software is used to establish a "Coal and Rock Mass Layered Physical and Mechanical Parameter Database". The data table fields include: borehole number, layer depth (start-end), lithology, Young's modulus (mean ± deviation), Poisson's ratio (mean ± deviation), tensile strength (mean ± deviation), shear strength (cohesion, friction angle, mean ± deviation), porosity (inversion value ± correction deviation), fracture development density (inversion value ± correction deviation), test date, test personnel, etc. For example, the data record for the 5-7m layer of borehole No. 1 is: ZK1, 5-7m, sandstone, 25.0±0.5GPa, 0.25±0.02, 5.2±0.3MPa, c=3.5±0.2MPa, φ=30±1°, 38.0±1.0%, 3±0.5 lines / m³, 2025-XX-XX, XXX.

[0112] In step S2, the specific process for constructing the hierarchical discrete element particle assembly-fracture dual-network model is as follows:

[0113] Step S21: Based on the layering parameters obtained in step S1, set the domain in the PFC software according to the actual layering thickness of the coal and rock mass, and generate an independent wall space for each layer.

[0114] Specifically, determine the planar dimensions of the simulation area based on engineering requirements (usually taking the influence range of the fracturing borehole, such as 20m in the x direction and 20m in the y direction), and determine the z-direction dimensions based on the layer thickness (sandstone layer z: 0-3m, coal seam z: 3-8m, mudstone layer z: 8-10m). Input the Domain range through “Model→Domain→Set”: x∈[0,20]m, y∈[0,20]m, z∈[0,10]m;

[0115] Use the "Wall→Create→Box" command to generate independent wall boundaries for each layer:

[0116] For example, a sandstone wall: x∈[0,20], y∈[0,20], z∈[0,3], the wall property is set to "rigid", and the label is named "Sandstone_Wall";

[0117] Coal seam wall: x∈[0,20], y∈[0,20], z∈[3,8], labeled as “Coal_Wall”;

[0118] Mudstone wall: x∈[0,20], y∈[0,20], z∈[8,10], labeled as “Mudstone_Wall”;

[0119] By grouping the different layers of walls using "Wall→Group", it is easier to define specific areas when generating particles later, thus avoiding cross-mixing of particles from different layers.

[0120] Step S22: For different layers, based on their porosity and particle size distribution characteristics, generate random rigid particles with differentiated diameters (particle size ratio is adjusted according to the lithology of the layer, 1:3 for sandstone layers and 1:2 for coal seam layers), to form the basic network of particle aggregates.

[0121] Step S23: The fracture development pattern is obtained by inverting the reference ground-penetrating radar. Fractures with different orientations, dip angles and lengths are pre-set in each layer. The fracture width is set according to the particle diameter at a ratio of 1:5 to 1:8 to form a particle aggregate-fracture dual network model.

[0122] In step S22, the specific process for constructing the basic network of the particle assembly is as follows:

[0123] Step S221: Extract key parameters of each layer from the database in Step S1 (e.g., the No. 3 coal seam area of ​​a certain coal mine is divided into: sandstone layer (thickness 3m), coal seam (thickness 5m), mudstone layer (thickness 2m) from top to bottom), determine the particle parameters of the sandstone layer, coal seam and mudstone layer, establish a three-dimensional coordinate system with the center point of the bottom surface of the simulation area as the origin (0,0,0), the z-axis is perpendicular to the strike of the strata (representing the depth direction), the x-axis is parallel to the strike of the coal seam, and the y-axis is perpendicular to the strike of the coal seam (representing the horizontal direction); set the particle contact model to linear contact by default, and unify the particle density to 2500 kg / m³;

[0124] In practical implementation, the key parameters for each layer are as follows:

[0125] Sandstone layer: porosity 12%, particle size distribution is coarse sand (refer to the specification for a particle size range of 2-6mm), minimum particle size is set according to a 1:3 particle size ratio. =2mm, maximum particle size =6mm;

[0126] Coal seam: porosity 18%, particle size distribution is pulverized coal-fine coal (particle size range 1-2mm), set at a particle size ratio of 1:2. =1mm =2mm;

[0127] Mudstone layer: porosity 8%, particle size distribution is fine mudstone particles (particle size range 3-6mm), set at a particle size ratio of 1:2. =3mm =6mm;

[0128] Step S222: Activate the sandstone layer wall group, input the particle size range and calculate the number of particles through porosity, and generate particles in the sandstone layer Domain (z∈[0,3]);

[0129] Specifically, input particle size range =2mm =6mm, the number of particles is calculated by back-calculating the porosity (formula: total particle volume = domain volume × (1 - porosity), the volume of a single particle is calculated based on the average particle size, sandstone layer domain volume = 20 × 20 × 3 = 1200 m³, total particle volume = 1200 × (1 - 12%) = 1056 m³, average particle size 4mm, single particle volume ≈ The number of particles is approximately 3.15 × 10¹ 0 indivual.

[0130] Step S223: Repeat step S222 to generate particles in the coal seam and mudstone layer respectively, and determine the volume of particles in each layer;

[0131] Repeat the above operation to generate particles in the coal seam (z∈[3,8]) and mudstone layer (z∈[8,10]), respectively. The number of particles in the coal seam is approximately 1 / 2. (Average particle size 1.5mm, individual volume ≈) ), number of mudstone particles (Average particle size 4.5 mm, individual volume) .

[0132] Step S224: Gravity compaction is performed on each layer of generated particles until the stress of the particle system is stable; after compaction, the particle overlap rate is checked; if the particle overlap rate is greater than 3%, the particle positions are readjusted to ensure that the particle aggregate structure is stable and conforms to the actual dense state of coal and rock mass.

[0133] Step S225: Extract the fracture characteristics of each layer from the ground-penetrating radar inversion data in step S1, and determine the fracture parameters;

[0134] Specifically, the characteristics of each layer of fractures are as follows:

[0135] Sandstone layer: sparsely developed fractures (affected by the density of the rock layer), with the fractures trending nearly horizontally (parallel to the bedding plane), dip angle 0-10°, length 0.5-1.5m, and 5-8 fractures pre-installed per cubic meter;

[0136] Coal seam: densely developed fractures (superimposed primary fractures and tectonic fractures), with two groups of fractures (one group parallel to the coal seam strike, dip angle 10-20°; one group perpendicular to the coal seam strike, dip angle 70-80°), with a length of 0.3-1.0m, and 15-20 fractures pre-installed per cubic meter;

[0137] Mudstone layer: extremely sparsely developed fractures (strong plasticity, fractures are easy to close), random strike, dip angle 30-50°, length 0.2-0.8m, 3-5 fractures pre-placed per cubic meter;

[0138] The fracture width is set according to the ratio of particle diameter 1:5 to 1:8: for sandstone layers, the average particle diameter is 4mm and the fracture width is 0.5-0.8mm; for coal seams, the average particle diameter is 1.5mm and the fracture width is 0.2-0.3mm; and for mudstone layers, the average particle diameter is 4.5mm and the fracture width is 0.6-0.9mm.

[0139] Step S226: Generate planar cracks based on the parameters of each layer of cracks;

[0140] Specifically, this simulates the extension of actual two-dimensional fractures in three-dimensional space, taking coal seams as an example:

[0141] Generate fractures parallel to the coal seam strike: Set the fracture plane normal vector to (0,1,0) (perpendicular to the y-axis, i.e. parallel to the xz plane), with a dip angle of 15° (angle with the x-axis), and randomly distribute them within the coal seam domain (x∈[0,20], y∈[0,20], z∈[3,8]). Each fracture is 0.8m long and 0.25mm wide, cutting through the particles within the fracture range to form fracture space;

[0142] To generate fractures perpendicular to the coal seam strike: Set the fracture plane normal vector to (1,0,0) (perpendicular to the x-axis, i.e. parallel to the yz plane), dip angle 75°, length 0.6m, and width 0.25mm, and similarly cut particles to generate fractures;

[0143] Repeat the above operation for sandstone and mudstone layers, generating fractures according to their respective fracture parameters, and avoiding fracture overlap during the generation process (overlap rate controlled within 5%).

[0144] Step S227: Observe the spatial distribution of particle aggregates and cracks, calculate the porosity of the model. If the deviation from the target value is greater than 3%, adjust the crack density or width until the porosity meets the requirements, and finally form a particle geometry-crack dual network model.

[0145] Specifically, the porosity of the model (the ratio of particle volume + fracture volume to domain volume) is calculated. The porosity of sandstone layer should be approximately 12% (particle porosity) + 2% (fracture porosity) = 14%, coal seam approximately 18% + 5% = 23%, and mudstone layer approximately 8% + 1% = 9%. If the deviation from the target value is greater than 3%, the fracture density or width is adjusted until the porosity meets the requirements, ultimately forming a particle aggregate-fracture dual-network model. The constructed dual-network model can clearly present the particle distribution and fracture morphology of each layer. When the contact force is calibrated in S3, the macroscopic mechanical response of the model (such as the uniaxial compressive strength of coal seam of 15MPa, with a deviation of only 1.3% from the laboratory test value of 14.8MPa) is highly consistent with the real coal and rock mass, laying a reliable foundation for subsequent multi-field coupled simulation and quantification of permeability enhancement effect.

[0146] In practice, the coordinate system origin is set to (0,0,0), x∈[0,20]m (parallel to the coal seam strike), y∈[0,20]m (perpendicular to the coal seam strike), and z∈[0,10]m (depth direction, z=0-3m is sandstone, 3-8m is coal seam, and 8-10m is mudstone); three sets of rigid walls are generated, each corresponding to one of the three layers, with clear labels to ensure independent boundaries;

[0147] The particle aggregate is generated as follows:

[0148] Sandstone layer: =2mm =6mm, generating 3.15×10¹ 0Each particle, after compaction, has a particle overlap rate of 2.1%, which meets the requirements; Coal seam: =1mm =2mm, generating 3.91×10¹ 0 Each particle, with an overlap rate of 1.8% after compaction; mudstone layer: =3mm =6mm, generating 1.09×10¹ 0 Each particle has an overlap rate of 2.5% after compaction.

[0149] The cracks are pre-set as follows:

[0150] Sandstone layer: 6 fractures are pre-installed per cubic meter, horizontal in direction (dip angle 5°), 1.0m in length and 0.6mm in width, generating a total of 20×20×3×6=7200 fractures. After cutting the particles, the fracture porosity is 1.8%.

[0151] Coal seam: 18 fractures are pre-installed per cubic meter, divided into two groups (parallel strike: dip angle 15°, length 0.8m; perpendicular strike: dip angle 75°, length 0.6m), 9 fractures / m³ in each group, generating a total of 20×20×5×18=36000 fractures, with a width of 0.25mm and a fracture porosity of 4.7%.

[0152] Mudstone layer: 4 fractures are pre-installed per cubic meter, with random orientation (dip angle 40°), length 0.5m, width 0.7mm, generating a total of 20×20×2×4=3200 fractures, with a fracture porosity of 0.9%;

[0153] The final model's total porosity is 13.8% for sandstone, 22.7% for coal seam, and 8.9% for mudstone, with deviations from the target values ​​all less than 2%. The dual-network model has been successfully constructed.

[0154] In step S3, the specific process for generating the bpm model of the coal and rock mass sample is as follows:

[0155] Step S31: Conduct simulation tests and record the stress-strain curve, peak strength, particle displacement field at failure, and crack propagation trajectory in real time for each sub-Domain; set the loading rate to 0.001 m / s for each sub-Domain (the same as the laboratory uniaxial compression test rate), and the loading direction is vertical (in the direction of ground stress).

[0156] Step S32: Simulate and generate simulated stress-strain curves for each layer, and extract simulated failure modes;

[0157] Step S33: Compare the simulation results with the uniaxial compression test results of the same layered specimen in the laboratory from Step S1, and calculate the deviation of key indicators. The formula for calculating the deviation of indicators is:

[0158] ;

[0159] Specific implementation example: The peak strength in the test was 30 MPa, the simulated peak strength was 28 MPa, and the deviation rate was 6.7%; the elastic modulus in the test was 3.7 GPa, the simulated elastic modulus was 3.5 GPa, and the deviation rate was 5.4%.

[0160] Step S34: Set the parameter adjustment rules and re-copy the single-axis compression simulation according to the adjusted parameters;

[0161] The parameter adjustment rules are as follows:

[0162] When the peak strength deviation is >5%: the tensile strength / shear strength is corrected by 1.1 times the deviation rate (e.g., if the peak strength deviation of the upper layer is 6.7%, the tensile strength is adjusted from 2.8MPa to 2.8×(1+6.7%×1.1)=3.0MPa); when the elastic modulus deviation is >5%: the Young's modulus is corrected by 1.05 times the deviation rate (e.g., if the elastic modulus deviation of the upper layer is 5.4%, the Young's modulus is adjusted from 3.5GPa to 3.5×(1+5.4%×1.05)=3.7GPa); if the failure modes do not match, the friction angle is adjusted (e.g., if the simulation is pure tensile failure and the test is tensile-shear failure, the friction angle is increased by 3°-5°).

[0163] Step S35: Recalculate the deviation. If the deviation of all indicators... If the calibration is successful, the calibration is complete; if not, the deviation calculation is repeated until the standard is met.

[0164] Step S36: Finally, accurate BPM models for each layer are formed.

[0165] In practice, the initial parameters were Young's modulus 3.5 GPa, tensile strength 2.8 MPa, and shear strength 5.2 MPa; the laboratory test results showed a peak strength of 30 MPa, an elastic modulus of 3.7 GPa, and a mixed tensile-shear failure mode; the initial simulation results showed a peak strength of 28 MPa (deviation 6.7%), an elastic modulus of 3.5 GPa (deviation 5.4%), and a failure mode predominantly tensile (deviation).

[0166] The iterative calibration process is as follows:

[0167] First adjustment: Young's modulus increased to 3.7 GPa, tensile strength increased to 3.0 MPa, and friction angle increased from 32° to 34°; First simulation results: peak strength 29.5 MPa (deviation 1.7%), elastic modulus 3.7 GPa (deviation 0%), failure mode tensile-shear hybrid (matched); deviation rate <5%, calibration completed, and the BPM model parameters of the upper coal seam determined.

[0168] In step S4, the actual apparent area of ​​the fluid domain is calculated using the bpm model of the coal and rock mass sample as follows:

[0169] Step S41: Based on the particle contact detection function of the DEM model, filter out particle pairs with valid contact conditions, such as contact type filtering, contact force threshold determination and output contact pair information.

[0170] Specifically, the contact type screening is as follows: only parallel bonding contact (corresponding to the cementation between coal and rock particles) and direct contact (corresponding to the interparticle gaps in the original pores of coal and rock) are retained, while sliding contact (contacts where relative displacement has occurred and there is no stable flow channel) is excluded.

[0171] The contact force threshold is determined as follows: Set a lower limit threshold for contact force (e.g., 0.1N, calculated based on the compressive strength of coal and rock mass obtained from S1 to ensure the stability of the contact structure), and retain only particle pairs with contact force ≥ the threshold to avoid false fluid channels caused by tiny contact gaps;

[0172] Output contact pair information. Using built-in model commands (such as contact.extract() in PFC3D), output the particle number of the effective contact pair (e.g., particle A: ID=101, particle B: ID=102), contact point coordinates (e.g., X=5.2mm, Y=3.8mm, Z=2.1mm), and contact gap width (e.g., 0.05mm, calculated from the difference between particle radius and center distance).

[0173] Step S42: For the effective disconnection pairs after screening, assign physical properties to the fluid pipes, establish a contact-pipe mapping relationship, and assign a unique identifier to each fluid pipe to avoid confusion between pipe and particle contact in subsequent simulations;

[0174] Specifically, establishing the contact-pipe mapping relationship includes the pipe inner diameter, pipe length, and initial pipe permeability value. The pipe inner diameter is directly related to the particle contact gap width, and is taken as 1.2 times the contact gap width (correcting for the actual channel reduction caused by particle surface roughness transition, verified by back-calculation based on the coal and rock mass porosity in step S1). For example, if the contact gap is 0.05 mm, then the pipe inner diameter = 0.06 mm. The pipe length is calculated as the length of the line connecting the centers of the two contacting particles (e.g., particle A center coordinates (5.0 mm, 3.5 mm, 2.0 mm), particle B center coordinates (5.4 mm, 4.1 mm, 2.2 mm), then... The initial value of pipeline permeability is calculated based on the Kozeny-Carman equation, and the formula is as follows: ,in, The inner diameter of the pipe. The porosity corresponding to this contact area (take the average porosity of the coal and rock mass layers in step S1, such as 8% for the coal seam). , hour, .

[0175] Each fluid pipe is assigned a unique identifier using the code "contact pair particle ID combination + contact sequence number". For example, the first effective contact between particles 101 and 102 is identified as "Pipe_101_102_001". In the model, the centers of the two contacting particles are connected by a red line segment, and the thickness of the line segment corresponds to the inner diameter of the pipe (e.g., an inner diameter of 0.06 mm corresponds to a line segment width of 0.03 mm), which makes it easy to visually view the pipe distribution and density.

[0176] Step S43: Randomly select a central particle C from the model, and extract all particles that have formed effective contact with the central particle C, denoted as the surrounding particle group D;

[0177] Step S44: Detect whether there are particles in contact with each other in the surrounding particle group D, and ensure that the central particle C and the surrounding particle group D can jointly form a closed area; if they cannot form a closed area, replace the central particle C.

[0178] In specific implementation, the particle range of the fluid domain is defined using a central particle plus surrounding contacting particles as the basic unit: particles are randomly selected from the model (particles with ≥3 contacting particles are preferred to ensure that a closed area can be formed), such as selecting a central particle C (ID=201, radius r=1.5mm, center coordinates (10.0mm, 8.0mm, 6.0mm)); all particles that have effective contact with the central particle C (i.e., contact pairs of the central particle C included in step S41) are extracted and denoted as the surrounding particle group D (e.g., particles D1: ID=202, D2: ID=203, D3: ID=204, a total of 3 particles); it is checked whether there are particles in mutual contact in the surrounding particle group D (e.g., D1 and D2, D2 and D3, D3 and D1 are all in effective contact), to ensure that the central particle C and the surrounding particles D can jointly form a closed area. If the number of contacting particles is <3 or no mutual contact is formed, the central particle is replaced.

[0179] Step S45: Based on the particle center coordinates and radius, determine the vertex coordinates of the polygon (2D) or polyhedron (3D) of the fluid domain;

[0180] Step S46: Determine the flow volume or area based on the vertex coordinates of the polygon or polyhedron of the fluid domain, and mark the fluid pipes connected to each fluid domain;

[0181] Specifically, in order to assign flow-related properties to fluid domains and establish their association with fluid pipes, the calculation is temporarily based on the geometric boundaries of polygons (or polyhedra) (to be corrected in subsequent step S47). For example, the area of ​​polygons P1-P6 in the 2D model is temporarily denoted as S_temp. The fluid pipes connected to each fluid domain are marked (e.g., fluid domain F1 connects to Pipe_201_202_001, Pipe_201_203_001, and Pipe_201_204_001) to ensure that the fluid can flow continuously between the pipes and the domains. In the model, the polygon boundaries of the fluid domains are filled with blue semi-transparent surfaces with a transparency of 50% to create a visual distinction from the red fluid pipes, which facilitates the verification of the connectivity between the domains and the pipes.

[0182] Step S47: Extract particle coordinates and boundary correction parameters from the DEM model, and calculate the actual apparent area;

[0183] Specifically, the center coordinates (e.g., particle C: (10.0, 8.0, 6.0) mm, D1: (9.0, 7.5, 6.0) mm) and radii (e.g., C: 1.5 mm, D1: 1.4 mm) of all particles within the output fluid domain are calculated. Based on the scanning electron microscope (SEM) image analysis of the coal and rock mass in step S1, the proportion of micropores (diameter < 0.01 mm) is statistically analyzed and used as the blockage coefficient. (e.g., coal seam extraction) =0.05, meaning 5% of the area is blocked and cannot flow); the particle surface roughness correction value is calculated by measuring the surface roughness of coal and rock particles using atomic force microscopy (AFM) and converting it into an area correction factor. (like =0.95, meaning the actual effective area is 95% of the geometric area.

[0184] In practical implementation, taking a 2D model as an example, the geometric area is calculated using the "shoelace formula," and then the actual apparent area is obtained by combining it with a correction factor. The formula and steps are as follows:

[0185] Step 1: List the coordinates of the vertices of the fluid domain (in order);

[0186] Let the vertices of the fluid domain be in order. ,and (Closed polygon), in this embodiment n=6, coordinates are as follows: , , , , , , ;

[0187] The formula for calculating the geometric area of ​​shoelace is as follows: Substituting the data, the calculation is as follows:

[0188] ;

[0189] ;

[0190] ;

[0191] ;

[0192] ;

[0193] ;

[0194] Summation: ;

[0195] ;

[0196] Let the corrected actual apparent area be... The specific formula is as follows: In the formula, =0.05, =0.95; Substituting into the formula, we get: =2.625×(1-0.05)×0.95≈2.625×0.95×0.95≈2.36mm²;

[0197] In the 3D model, the fluid domain is a polyhedron. The "tetrahedral decomposition method" is used to decompose the polyhedron into multiple tetrahedra. The volume of each tetrahedron is calculated and then summed to obtain the geometric volume. Then through = The actual apparent volume is obtained by multiplying (1-f) by k, which is then used for subsequent 3D seepage simulation.

[0198] In step S45, when determining the vertex coordinates of the polygon or polyhedron in the fluid domain, first, draw the common external tangents of the central particle C and the surrounding particle D1 to obtain the intersection point P1 of the two common external tangents. Then, calculate the intersection point P2 of the common external tangents of the central particle C and the surrounding particle D2, and the intersection point P3 of the common external tangents of the central particle C and the surrounding particle D3. Next, calculate the intersection point P4 of the common tangents of the surrounding particles D1 and D2, and similarly calculate the intersection point P5 of the common tangents of the surrounding particles D2 and D3, and the intersection point P6 of the common tangents of the surrounding particles D3 and D1. Arrange the intersection points P1-P6 in a clockwise or counterclockwise order to form the vertex coordinates of the closed polygon (e.g., ...). , , , , , This ensures that the vertices are connected sequentially without crossing.

[0199] In step S5, the specific procedure for assigning seepage hydraulic parameters to the bpm model of the coal and rock mass sample is as follows:

[0200] Step S51: Obtain the initial opening of the fluid conduit, the half-open compressive force, the apparent volume of the fluid domain, the residual opening magnification factor of the fracture, the bulk modulus of the fluid, and the fluid viscosity from the bpm model of the coal and rock mass sample;

[0201] Specifically:

[0202] 1. Initial opening of fluid pipeline

[0203] The distribution of pore throats in coal and rock mass samples was tested using laboratory mercury intrusion porosimetry to obtain the average pore throat width (denoted as ) of the target layer (e.g., coal seam No. 3). ), assuming the test result is =0.04mm; Since the particles in the DEM model are rigid spheres, but the actual coal and rock mass particles have micro-protrusions on their surfaces, a correction factor needs to be introduced. ( =1.1-1.3 (α takes a larger value if the porosity is high), this process takes... =1.2; therefore, the calculation formula is: Substituting the data, we get: Meanwhile, in the PFC3D software, the fish function is used to traverse all fluid pipes (pipe identifiers marked with S4), and the pipe.openness parameter is assigned a value of 0.048mm.

[0204] 2. Compression force at half-opening angle

[0205] The half-open compressive force reflects the change in opening degree when a fluid pipe is subjected to confining pressure or particle compression. Down to The compressive force required in the (partially open) state reflects the pipe's ability to resist extrusion and deformation, and directly affects the change in seepage resistance.

[0206] The specific value acquisition process is as follows: Extract the elastic modulus of coal seam No. 3 from the layered parameter database in step S1. Poisson's ratio Based on Hertzian contact theory, the relationship between pipe compressibility and opening change is as follows: ,in, The average radius of the contacting particles (the average radius of the particles in coal seam No. 3 in step S2 is r = 1.5 mm). The change in opening ( = -0.5 =0.024mm);

[0207] Substitute the data, In the model, a value is assigned to the "Contact Mechanical Parameters" module for each fluid pipe. The contact.pipe.half_open_force is set to 54.2N to ensure that the pipe opening changes with the compressive force during subsequent confining pressure loading in accordance with the actual mechanical laws.

[0208] 3. Apparent volume of the fluid domain

[0209] The effective volume of fluid that can be contained in all fluid domains in the model takes into account both fluid storage capacity and flow channel continuity. In the 3D model, it is volume, and in the 2D model, it is area (which needs to be multiplied by the model thickness to convert to volume).

[0210] The specific value retrieval process is as follows: Call the total actual apparent area of ​​the No. 3 coal seam model calculated by S4. (2D model, size 50mm×50mm); If it is a 3D model, the model thickness needs to be set. (Determined based on the core drilling diameter, e.g., h=50mm, consistent with the model's planar dimensions), perform thickness correction for the 3D model; the calculation formula is as follows: ,in, The correction factor for pore connectivity (based on laboratory CT scans of pore connectivity, taken from coal seam No. 3). (excluding isolated pores) to obtain The total apparent volume is determined through the "fluid domain management module" of the model. The fluid volume is allocated to each fluid domain (based on the area ratio of a single fluid domain; for example, if the area of ​​a fluid domain is 2.36 mm², the ratio is 2.36 / 486.2 ≈ 0.48%), then its volume...

[0211] 4. Residual crack aperture magnification factor

[0212] When a fracture closes (unloads) after hydraulic fracturing, the ratio of the residual aperture to the initial aperture reflects the ability of the coal and rock mass to "maintain" the fracture due to plastic deformation. >1 indicates that the residual aperture is greater than the initial aperture, and the anti-reflective effect is long-lasting.

[0213] The fracture state of the No. 3 coal seam after fracturing was observed using a downhole fracture inspection instrument, and the initial fracture aperture (right after fracturing) at 5 observation points was recorded. Residual opening (after 24 hours of unloading) ;

[0214] Pick The average value is 0.72. Considering the idealization of crack propagation in the model (without complex geological interference in the field), a correction factor is introduced. (Slightly magnified to match actual communication retention capability), then ;

[0215] Model assignment operation: In the "Fractured Mechanical Parameters" module, set fracture.residual_opening_coeff to 0.76, and set the trigger condition: when the pressure on the fracture drops to 30% of the initial fracturing pressure (unloading threshold), automatically update the fracture aperture to "current aperture ×". ".

[0216] 5. Fluid bulk modulus

[0217] The ability of a fluid to resist volume compression under pressure is expressed by the formula: The larger the value, the more difficult the fluid is to compress, and the more stable the pressure transmission is during the seepage process.

[0218] The specific value retrieval process is as follows:

[0219] In-situ fracturing, clean water (20℃) was used. Referring to the *Handbook of Fluid Mechanics*, the standard value of the bulk modulus of clean water at 20℃ is 2.1 GPa. In step S1, the geothermal temperature of coal seam No. 3 is 35℃ (calculated based on a geothermal gradient of 3℃ / 100m and a burial depth of 500m). For every 1℃ increase in temperature, the bulk modulus of clean water decreases by 0.004 GPa. The correction formula is as follows: ;

[0220] Model assignment operation: In the "Fluid Properties" module, directly set fluid.bulk_modulus to 2.04GPa. This parameter is a global parameter and is shared by all fluid domains and pipelines.

[0221] 6. Fluid viscosity

[0222] The internal friction between fluid molecules determines the ease or difficulty of fluid flow, as shown in the formula: ( For shear stress, (for velocity gradient) The larger the value, the greater the flow resistance.

[0223] Value retrieval process:

[0224] Basic testing: The viscosity of the fracturing fluid (water + 0.5% drag reducer) at 20℃ was tested using a rotational viscometer. Test results... ;

[0225] Temperature and pressure correction:

[0226] Temperature correction: At a ground temperature of 35℃, the viscosity of clean water increases with temperature according to... Calculate and substitute to get ;

[0227] Pressure correction: The coal seam is buried at a depth of 500m, with a pressure of approximately 5MPa. For every 1MPa increase in pressure, the viscosity increases by 0.02mPa·s. After correction... ;

[0228] Model assignment operation: In the "Fluid Properties" module, set fluid.viscosity to 1.433 mPa・s, which works together with the fluid bulk modulus on the seepage equation.

[0229] Step S52: Perform correlation verification and correction on the parameters;

[0230] After assigning values ​​to individual parameters, it is necessary to verify the synergy between the parameters to ensure that they conform to the actual seepage patterns and avoid simulation distortion caused by parameter inconsistencies. The verification process is as follows:

[0231] Initial seepage velocity verification: based on Darcy's law ,in Calculated from the initial pipe opening ( (Derivation of the formula for flow from a parallel flat plate). Take 5MPa (inlet pressure in step S6). Take a model length of 50mm and substitute the parameters:

[0232] ;

[0233] ;

[0234] Compare the initial seepage velocity of the in-situ borehole water injection test , Within reasonable limits, the verification passed.

[0235] Seepage verification after fracture closure: Assuming an initial fracture aperture of 0.3 mm after fracturing, the residual aperture after unloading... Calculate residual permeability Compare the residual permeability tested after in-situ fracturing. The verification passed.

[0236] Correction mechanism: If the verification deviation exceeds 15% (e.g., in calculating seepage velocity) The scene is If the deviation is 50%, then prioritize adjusting the fluid viscosity (which has the highest impact weight) or the initial pipe opening until the deviation is less than 10%.

[0237] Step S53: Apply a fixed inlet pressure p1 to one side of the coal and rock mass sample bpm model and a fixed outlet pressure p2 to the other side, where p1>p2, to create a pressure difference and drive the fluid to flow in the model; wait for the average pore pressure and total flow rate of the coal and rock mass sample bpm model to reach a steady state. The criterion is that the rate of change of the average pore pressure and total flow rate is less than the set threshold within several consecutive calculation steps.

[0238] Step S54: When a steady state is reached, calculate the permeability k0 of the coal and rock mass sample bpm model according to Darcy's law;

[0239] Specifically, when a steady state is reached, the permeability k0 of the coal-rock sample bpm model is calculated according to Darcy's law. The expression for Darcy's law is:

[0240] ;

[0241] in, The flow rate under steady-state conditions can be calculated using a model; This represents the length of the fluid flow, i.e., the distance between the inlet and outlet in the model; The fluid viscosity has been assigned in step S5; The flow area of ​​the fluid can be calculated using the geometric information of the fluid pipes and fluid domains in the model; The pressure difference between imports and exports, i.e. By substituting the values ​​of these parameters, the initial penetration rate of the model can be calculated. .

[0242] The core of hydraulic fracturing simulation is to construct a dynamically coupled process of "fluid injection - stress transfer - particle movement - fracture propagation" using DEM software. This requires overcoming the limitations of traditional "fixed parameter injection" and achieving a combination of "pressure-flow coordinated control" and "multi-state real-time monitoring." This process follows the main thread of "model depressurization reset → engineering parameter mapping → dynamic injection control → full-dimensional data recording," ensuring a high degree of consistency between the simulation process and the timing and mechanical response laws of on-site fracturing operations, providing complete process data support for subsequent permeability enhancement effect analysis.

[0243] In step S6, the specific process of simulating the actual hydraulic fracturing process is as follows:

[0244] Step S61: Before simulation, perform model pressureless reset and initial state calibration. Based on the injection hole drilled in step S5, redefine the boundary conditions for fracturing injection.

[0245] Specifically, the confining pressure set in S6 (e.g., low confining pressure condition: vertical 10MPa, horizontal 8MPa) is reduced to 0MPa using the DEM software in a linear unloading manner, with the unloading rate set to 0.5MPa / hour (1 hour step). To avoid violent vibration of model particles caused by instantaneous unloading; close the inlet pressure source and outlet pressure boundary set in S5, and reduce the initial pore pressure of all fluid domains in the model to atmospheric pressure (0.1MPa) through the fluid pressure release command (such as fluid.pressure.release() in PFC3D) to simulate the initial pore pressure state of the formation; after unloading, run the model for 500 time steps, and check the average velocity of the particles through the software. When the average velocity is <0.01mm / time step, it is determined that the model has reached a pressureless equilibrium state to avoid residual displacement affecting the fracturing simulation;

[0246] Based on the drilled injection hole (e.g., 10mm diameter, center coordinates (25mm, 25mm, 25mm), 3D model), the boundary conditions for fracturing injection are redefined: the injection hole range is marked in the model, the particles inside the injection hole are set as removable particles (simulating coal and rock mass removed during on-site drilling), and the particles on the injection hole wall (ID range: 1000-1200) are retained as the contact boundary for fluid injection; the contact area of ​​the particles on the injection hole wall is set as the fluid injection surface, and a "pressure-flow dual control" attribute is assigned, that is, during the injection process, it can automatically switch between "constant pressure injection" and "constant flow rate injection" according to the model response; the injection hole is filled with a fluid consistent with the S5 parameters (viscosity 1.433mPa・s), and the filling volume is 80% of the injection hole volume (to avoid initial overfilling leading to a sudden pressure rise). After filling, the pressure inside the hole is monitored to ensure that the initial pressure is stable at 0.1MPa (consistent with the formation pore pressure).

[0247] Step S62: Based on the on-site fracturing construction data, determine the injection parameters for fracturing and set dynamic control parameters;

[0248] Specifically, referring to the on-site fracturing construction data of the No. 3 coal seam in a certain mine (daily water injection volume 80m³, single-hole injection time 12h), the injection volume per unit time is calculated as follows: Combined with model size ( The scaling ratio of the actual volume on site () ), calculate the injected flow of the model: Converted to commonly used units within the model (mm³ / s): Finally, the model injection traffic was set to be (Reserve a fluctuation range);

[0249] Based on the tensile strength (1.2 MPa) of coal seam No. 3 in S1, and considering the need to overcome in-situ stress and coal-rock mass strength during fracturing, the upper limit of the injection pressure is set to three times the tensile strength (3.6 MPa) to avoid excessive pressure causing instantaneous model failure and failing to reflect the fracture propagation process. The fracturing fluid parameters set in S5 (clean water + 0.5% drag reducer, viscosity) are retained. bulk modulus This ensures that the fluid properties in the simulation are consistent with those in reality.

[0250] When the maximum pore pressure in the model reaches 90% of the upper limit of the injection pressure (3.24 MPa), it automatically switches from "constant flow rate injection" to "constant pressure injection" to avoid the model from collapsing due to the pressure exceeding the upper limit; when the crack extends to the model boundary (forming a through crack) and the pore pressure drops to 2.0 MPa, it switches back to "constant flow rate injection" to maintain stable crack extension.

[0251] Simulation duration setting: Refer to the on-site fracturing fracture propagation rate ( ), combined with model size ( ), calculate the time required for the crack within the model to propagate to the boundary: Considering the time required for particle movement and crack stabilization, the total simulation duration is set to 500 seconds (corresponding to the model time step). (Time step).

[0252] Step S63: Use DEM software to dynamically simulate the fracturing process, adjust the injection state in real time, and ensure that the fracture propagation pattern is consistent with the field.

[0253] The fracturing process is dynamically simulated using the "fluid-particle coupling module" of DEM software, and the injection state is adjusted in real time to ensure that the fracture propagation law is consistent with the field. The steps are as follows:

[0254] Step 1, Initial Injection Stage (0-100s, Constant Flow Rate Injection): Fluid is injected into the injection hole at a constant flow rate of 1.85 mm³ / s. The software automatically guides the fluid through the contact channels between the particles on the injection hole wall and the model interior. The average and maximum pore pressure within the model are recorded every 10s. In the initial stage (0-50s), the pressure rises slowly (from 0.1 MPa to 1.5 MPa). During this stage, the fluid mainly fills the initial pores within the model, and no new cracks are generated. The stress distribution of the particles is observed through the particle stress cloud diagram. When the maximum tensile stress of the particles around the injection hole (ID=1000-1200) reaches 1.0 MPa (close to the tensile strength of 1.2 MPa), the crack initiation stage begins.

[0255] Step 2, Crack Initiation and Propagation Stage (100-300s, Pressure-Flow Switching): When the maximum pore pressure rises to 3.24MPa (90% of the upper limit of the injection pressure), the software automatically switches to a constant pressure injection of 3.6MPa. At this time, the contact adhesion force of the particles (ID=1500) on the right side of the injection hole reaches the threshold (based on the parallel adhesion strength calibrated by S3), and the adhesion breaks, generating the first microcrack (2mm in length and 0.08mm in width). During crack propagation (150-300s), constant pressure injection continues, and the fluid penetrates along the first microcrack, causing the particles at the crack tip to be stressed and triggering a new adhesion fracture. The crack propagates to the right side of the model at a speed of 0.1mm / s. During this period, the crack status (length, width, direction) is recorded every 50s. When the crack propagates to the model boundary (50mm), the pore pressure inside the model drops from 3.6MPa to 2.0MPa, and the software automatically switches back to constant flow injection (1.85mm³ / s) to maintain a stable crack width (increasing from 0.2mm to 0.3mm).

[0256] Step 3, Fracture Stabilization Stage (300-500s, Constant Flow Rate Injection): Inject continuously at a flow rate of 1.85 mm³ / s. The fluid flows along the penetrating fracture. The pore pressure in the model stabilizes at 2.0-2.2 MPa, and the fracture width stabilizes at 0.3-0.32 mm. No new fractures are generated. When the simulation duration reaches 500s, or the fracture width remains unchanged for 100s (fluctuation < 0.01 mm), stop the injection to complete the fracturing simulation.

[0257] Step S64: Implement multi-dimensional data recording to ensure coverage of particle state, crack state, and fluid state.

[0258] Specifically, at 10-second intervals (corresponding to 1×10 5 The system records global data, including the model's average pore pressure, total injection flow rate, average particle displacement, and total crack length. Event-triggered recording: When key events such as "new crack generation" (bonding fracture count ≥ 5 times), "crack propagation direction change" (propagation angle change ≥ 15°), and "pressure drop" (pore pressure drop ≥ 0.5 MPa / s) occur, high-frequency recording (1-second interval) is automatically triggered to record local data at the time of the event (such as particle ID, crack coordinates, and pressure distribution in the event area).

[0259] In step S7, the permeability k0 and permeability k1 are compared to quantify the permeability enhancement effect of the bpm model hydraulic fracturing on the coal and rock mass sample. The permeability enhancement effect is measured by the rate of change of permeability, and the calculation formula is: Permeability enhancement = The higher the permeability enhancement rate, the more significant the permeability enhancement effect of hydraulic fracturing on coal and rock mass. By analyzing the influence of different factors (such as fracturing fluid pressure, flow rate, physical and mechanical parameters of coal and rock mass, etc.) on the permeability enhancement rate, we can gain a deeper understanding of the intrinsic mechanism of hydraulic fracturing permeability enhancement and provide scientific guidance for actual production operations.

[0260] Specifically, after the fracturing simulation, the model retains temporary data such as particle displacement, fracture morphology, and fluid pressure. A systematic reset operation is required to ensure the model retains only the permanent fractures formed by fracturing (reflecting the actual rock mass state after fracturing), while maintaining all other conditions identical to S6, thus avoiding interference with the permeability retest results. It is also crucial to ensure that the test conditions (confining pressure, pressure difference, stability criterion) for k1 and k0 are completely consistent to guarantee data comparability. The steps are as follows:

[0261] The nine test conditions in S6 are fully reused (3 sets of confining pressure × 3 sets of inlet pressure), and the parameters are shown in the table below (taking coal seam No. 3 as an example):

[0262] serial number Containment pressure (horizontal / vertical) Import pressure P1 Export pressure P2 Pressure difference ΔP(P1-P2) 1 10MPa / 8MPa (low) 5MPa 0.1MPa 4.9MPa 2 10MPa / 8MPa (low) 10MPa 0.1MPa 9.9MPa 3 10MPa / 8MPa (low) 15MPa 0.1MPa 14.9MPa 4 20MPa / 15MPa (Medium) 5MPa 0.1MPa 4.9MPa 5 20MPa / 15MPa (Medium) 10MPa 0.1MPa 9.9MPa 6 20MPa / 15MPa (Medium) 15MPa 0.1MPa 14.9MPa 7 30MPa / 25MPa (High) 5MPa 0.1MPa 4.9MPa 8 30MPa / 25MPa (High) 10MPa 0.1MPa 9.9MPa 9 30MPa / 25MPa (High) 15MPa 0.1MPa 14.9MPa

[0263] A confining pressure of 20 MPa vertically and 15 MPa horizontally is applied to the model at a rate of 0.5 MPa / hour step. After loading is completed, the model is stabilized for 200 hours step to ensure that the confining pressure is uniformly transferred to the particles.

[0264] An inlet pressure of 10 MPa is applied to the left side of the model (x=0 mm), and an outlet pressure of 0.1 MPa is applied to the right side (x=50 mm), creating a pressure difference of 9.9 MPa. ;

[0265] The model's average pore pressure and total flow rate are monitored in real time. When the average pore pressure fluctuation is <2% and the total flow rate fluctuation is <3% within 300 consecutive time steps, a steady state is considered reached (e.g., in condition 5, the steady-state average pore pressure is 5.05 MPa and the steady-state total flow rate is...). );

[0266] The formula for calculating k1 using Darcy's law is as follows: In the formula: For steady-state total flow ( ); For fluid viscosity ( (Same as step S5) The model flow path length is 50mm = 0.05m, x-direction dimension. The actual apparent area of ​​the fluid domain (calculated in step S4) ); For the pressure difference (9.9MPa=9.9×10), 6 Pa); Substitute into the calculation: Permeability tests were completed under nine different operating conditions, and the post-fracture permeability matrix was obtained. (Some data is shown in the table below):

[0267] serial number Pre-fracking permeability Permeability after fracturing 1 5 9

[0268] Transparency enhancement using formula Calculate the permeability enhancement rate for each working condition, reflecting the extent to which fracturing improves permeability. Take working condition 5 as an example:

[0269] ;

[0270] ;

[0271] Under operating condition 5, hydraulic fracturing increases the permeability of the coal and rock mass by approximately 190%.

[0272] A "reflection uniformity coefficient C" is introduced to evaluate the stability of the reflection enhancement effect under different working conditions. The formula is as follows: ,in: The standard deviation of the light transmission rate for the nine working conditions; The average value of the anti-reflection rate under 9 working conditions;

[0273] The calculation example is as follows:

[0274] The light transmission enhancement rates for the nine working conditions were 178%, 185%, 192%, 165%, 190%, 195%, 152%, 160%, and 170%, respectively; the average value was... Standard deviation ;but ( The closer it is to 1, the better the uniformity of antireflection.

[0275] It is worth noting that the various units included in the above system embodiments are only divided according to functional logic, but are not limited to the above division, as long as the corresponding functions can be achieved; in addition, the specific names of each functional unit are only for easy differentiation and are not used to limit the scope of protection of the present invention.

[0276] Furthermore, those skilled in the art will understand that all or part of the steps in the methods of the above embodiments can be implemented by a program instructing related hardware, and the corresponding program can be stored in a computer-readable storage medium.

[0277] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM, characterized in that, Includes the following steps: Step S1: Use a multi-source data fusion method of "on-site borehole sampling + laboratory precision testing + ground-penetrating radar inversion" to obtain the layered physical and mechanical parameters of coal and rock mass; Step S2: Construct a hierarchical discrete element particle assembly-fracture dual-network model; Step S3: Apply parallel bonding contact force to the particle aggregate model. The parallel bonding model enhances the connection strength between particles by introducing virtual bonds, generating a coal and rock mass sample bpm model. Step S4: Define the fluid conduit as the contact part between particles in the particle aggregate, and the fluid domain as the polygonal closed area enclosed by the center points of adjacent and contacting particles. Calculate the actual apparent area of ​​the fluid domain using the bpm model of the coal and rock mass sample. Step S5: Assign the seepage hydraulic parameters to the bpm model of the coal and rock mass sample and calculate the initial permeability; Step S6: Remove the confining pressure, inlet pressure, and outlet pressure applied to the coal and rock mass sample bpm model in step S5, restoring the model to its initial pressureless state. Apply a fluid with a certain flow rate and pressure value to the coal and rock mass sample bpm model from the injection hole drilled in step S5 to simulate the actual hydraulic fracturing process. Step S7: Remove the calculation results assigned to the bpm model of the coal and rock mass sample in Step S6, i.e., clear the temporary data and state changes generated during the hydraulic fracturing simulation, and restore the model to the initial bpm model state; repeat Step S6, assign a fixed confining pressure to the model again, drill the injection hole, and apply the inlet pressure. and export pressure After waiting for the average pore pressure and total flow rate of the model to reach a steady state, the permeability of the coal and rock mass sample under constant pressure steady-state condition (bpm model) is calculated. .

2. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S1, the specific process for obtaining the physical and mechanical parameters of coal and rock mass stratification is as follows: Step S11: Select 3 to 5 representative boreholes on site and take samples at 2m intervals in layers to ensure that different lithological sections of the coal and rock mass are covered; Step S12: The laboratory uses a servo press and a direct shear tester to test the Young's modulus, tensile strength, shear strength, friction angle, friction coefficient and Poisson's ratio of each layer of specimens. Each parameter is tested three times and the average value is taken. Step S13: Use ground-penetrating radar to scan the coal and rock mass within a 10m radius around the borehole, invert the layered porosity and fracture development density, cross-validate with laboratory data, correct parameter deviations, and finally establish a coal and rock mass layered parameter database.

3. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S2, the specific process for constructing the hierarchical discrete element particle assembly-fracture dual-network model is as follows: Step S21: Based on the layering parameters obtained in step S1, set the domain in the PFC software according to the actual layering thickness of the coal and rock mass, and generate an independent wall space for each layer. Step S22: For different layers, based on their porosity and particle size distribution characteristics, generate random rigid particles with differentiated diameters to form the basic network of the particle aggregate. Step S23: The fracture development pattern is obtained by inverting the reference ground-penetrating radar. Fractures with different orientations, dip angles and lengths are pre-set in each layer. The fracture width is set according to the particle diameter at a ratio of 1:5 to 1:8 to form a particle aggregate-fracture dual network model.

4. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 3, characterized in that, In step S22, the specific process for constructing the basic network of the particle aggregate is as follows: Step S221: Extract key parameters of each layer from the database in step S1, determine the particle parameters of sandstone, coal seam and mudstone layers, and set the particle contact model to linear contact by default. Step S222: Activate the sandstone layer wall group, input the particle size range and calculate the number of particles through porosity, and generate particles in the sandstone layer domain; Step S223: Repeat step S222 to generate particles in the coal seam and mudstone layer respectively, and determine the volume of particles in each layer; Step S224: Gravity compaction of each layer of generated particles; After compaction, check the particle overlap rate; if the particle overlap rate is greater than 3%, readjust the particle positions. Step S225: Extract the fracture characteristics of each layer from the ground-penetrating radar inversion data in step S1, and determine the fracture parameters; Step S226: Generate planar cracks based on the parameters of each layer of cracks; Step S227: Observe the spatial distribution of the particle aggregate and the cracks, calculate the porosity of the model, and if the deviation from the target value is greater than 3%, adjust the crack density or width until the porosity meets the requirements, and finally form a particle geometry-crack dual network model.

5. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S3, the specific process for generating the bpm model of the coal and rock mass sample is as follows: Step S31: Perform simulation tests and record the stress-strain curve, peak strength, particle displacement field at failure and crack propagation trajectory in real time for each sub-Domain; Step S32: Simulate and generate simulated stress-strain curves for each layer, and extract simulated failure modes; Step S33: Compare the simulation results with the uniaxial compression test results of the same layered specimen in the laboratory in Step S1, and calculate the deviation of key indicators. Step S34: Set the parameter adjustment rules and re-copy the single-axis compression simulation according to the adjusted parameters; Step S35: Recalculate the deviation. If the deviation of all indicators... If the calibration is successful, the calibration is complete; if not, the deviation calculation is repeated until the standard is met. Step S36: Finally, accurate BPM models for each layer are formed.

6. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S4, the actual apparent area of ​​the fluid domain is calculated using the bpm model of the coal and rock mass sample as follows: Step S41: Based on the particle contact detection function of the DEM model, filter out particle pairs with effective contact conditions; Step S42: For the filtered valid disconnection pairs, assign physical properties to the fluid pipes, establish a contact-pipe mapping relationship, and assign a unique identifier to each fluid pipe; Step S43: Randomly select a central particle C from the model, and extract all particles that have formed effective contact with the central particle C, denoted as the surrounding particle group D; Step S44: Detect whether there are particles in contact with each other in the surrounding particle group D, and ensure that the central particle C and the surrounding particle group D can jointly form a closed area; if they cannot form a closed area, replace the central particle C. Step S45: Determine the vertex coordinates of the polygon or polyhedron of the fluid domain based on the particle center coordinates and radius; Step S46: Determine the flow volume or area based on the vertex coordinates of the polygon or polyhedron of the fluid domain, and mark the fluid pipes connected to each fluid domain; Step S47: Extract particle coordinates and boundary correction parameters from the DEM model, and calculate the actual apparent area.

7. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S45, when determining the vertex coordinates of the polygon or polyhedron of the fluid domain, firstly, draw the common external tangents of the central particle C and the peripheral particle D1 to obtain the intersection point P1 of the two common external tangents. Then, calculate the intersection point P2 of the common external tangents of the central particle C and the peripheral particle D2, and calculate the intersection point P3 of the common external tangents of the central particle C and the peripheral particle D3. Next, calculate the intersection point P4 of the common tangents of the peripheral particles D1 and D2. Similarly, calculate the intersection point P5 of the common tangents of the peripheral particles D2 and D3, and the intersection point P6 of the common tangents of the peripheral particles D3 and D1. Arrange the vertices P1-P6 in clockwise or counterclockwise order to form the vertex coordinates of the closed polygon.

8. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S5, the specific process for assigning seepage hydraulic parameters to the bpm model of the coal and rock mass sample is as follows: Step S51: Obtain the initial opening of the fluid conduit, the half-open compressive force, the apparent volume of the fluid domain, the residual opening magnification factor of the fracture, the bulk modulus of the fluid, and the fluid viscosity from the bpm model of the coal and rock mass sample; Step S52: Perform correlation verification and correction on the parameters; Step S53: Apply a fixed inlet pressure p1 to one side of the coal and rock mass sample bpm model and a fixed outlet pressure p2 to the other side, where p1>p2, to create a pressure difference and drive the fluid to flow in the model; wait for the average pore pressure and total flow rate of the coal and rock mass sample bpm model to reach a steady state. The criterion is that the rate of change of the average pore pressure and total flow rate is less than the set threshold within several consecutive calculation steps. Step S54: When a steady state is reached, calculate the permeability k0 of the bpm model of the coal and rock mass sample according to Darcy's law.

9. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S6, the specific process of simulating the actual hydraulic fracturing process is as follows: Step S61: Before simulation, perform model pressureless reset and initial state calibration. Based on the injection hole drilled in step S5, redefine the boundary conditions for fracturing injection. Step S62: Based on the on-site fracturing construction data, determine the injection parameters for fracturing and set dynamic control parameters; Step S63: Use DEM software to dynamically simulate the fracturing process, adjust the injection state in real time, and ensure that the fracture propagation pattern is consistent with the field. Step S64: Implement multi-dimensional data recording to ensure coverage of particle state, crack state, and fluid state.

10. The method for quantifying the permeability enhancement effect of hydraulic fracturing in rock masses based on DEM according to claim 1, characterized in that, In step S7, the permeability k0 and permeability k1 are compared to quantify the permeability enhancement effect of the hydraulic fracturing model on the coal and rock mass sample. The permeability enhancement effect is measured by the rate of change of permeability, and the calculation formula is: Permeability enhancement = .

Images