A method of grout injection control coupling fracture aperture variation with slurry particle deposition

Abstract

The application discloses a grouting control method coupling crack opening degree change and slurry particle deposition, and relates to the technical field of geotechnical engineering grouting. Firstly, the current crack opening degree is acquired, the slurry diffusion flow field and the space-time viscosity distribution are estimated based on a three-dimensional discrete crack network model and a space-time viscosity evolution model of slurry, then the slurry pressure is applied to the crack wall surface as a normal load, the crack dynamic deformation under the fluid-structure coupling is calculated, then the deposition particle volume fraction is determined based on the flow velocity, the viscosity and the deformation, and the crack network permeability and the effective opening degree are updated accordingly, finally, the iteration cycle is performed by judging whether the evolved permeability meets the plugging threshold value, and the optimal grouting duration is determined. The application can truly reflect the coupling effect of crack deformation and particle deposition in the grouting process, and realizes the accurate prediction and control of the grouting plugging effect.

Description

Technical Field

[0001] This application relates to the field of grouting technology in geotechnical engineering, and more specifically, to a grouting control method that couples changes in fracture aperture with grout particle deposition. Background Technology

[0002] In geotechnical engineering fields such as tunnel excavation, dam foundation treatment, and groundwater sealing, fractured rock masses are widespread and often serve as major channels for groundwater seepage, easily inducing geological disasters such as water and mud inrushes, seriously threatening engineering safety and structural stability. Grouting technology, as an effective means of reinforcing rock masses, sealing fractures, and reducing permeability, is widely used in engineering practice. The core objective of grouting engineering is to effectively fill and seal fractures through the diffusion and solidification of grout in the fracture network. Accurately predicting the diffusion law and sealing effect of grout in complex fracture networks is a key prerequisite for formulating scientific grouting plans and determining the optimal grouting duration.

[0003] However, existing grouting control and simulation methods still have significant limitations when dealing with grouting processes involving complex fracture networks. Traditional grouting theoretical models or numerical simulation methods are usually based on the assumption of rigid fractures, which assumes that the fracture aperture remains constant during grouting. However, under actual high-pressure grouting conditions, the grout pressure acts on the fracture walls, producing splitting or expansion effects, leading to dynamic deformation of the fracture aperture. This fluid-structure interaction effect significantly alters the geometric characteristics of the grout flow channels. Ignoring the dynamic changes in fracture aperture often results in significant deviations in the prediction of grout diffusion range and grout volume, making it difficult to accurately reflect the mechanical response during the grouting process. Summary of the Invention

[0004] To address the aforementioned problems in existing technologies, this application provides a grouting control method that couples fracture aperture changes with grout particle deposition, comprising: S1: obtaining the current fracture aperture; S2: estimating the grout diffusion flow field and spatiotemporal viscosity distribution based on a three-dimensional discrete fracture network model and the current fracture aperture, combined with a grout spatiotemporal viscosity evolution model, to obtain the overall grout pressure distribution, grout velocity distribution, and real-time viscosity field; S3: applying the overall grout pressure distribution as a normal load to the fracture wall boundary and calculating the dynamic deformation of the fracture under fluid-structure interaction to obtain the fracture... S4: Determine the volume fraction of deposited particles based on the slurry flow velocity distribution, real-time viscosity field, and fracture aperture deformation; S5: Update the fracture network permeability and effective aperture according to the volume fraction of deposited particles to obtain the evolved fracture network permeability and the corrected fracture aperture, wherein the corrected fracture aperture is used as the updated current fracture aperture; S6: Determine whether the evolved fracture network permeability meets the preset sealing threshold. If it does not meet the threshold, the corrected fracture aperture is fed back to steps S2 to S5 for iterative looping until the sealing threshold is met to determine the optimal grouting time.

[0005] This application also provides a grouting control system that couples fracture aperture changes with slurry particle deposition, comprising: a fracture aperture acquisition module for acquiring the current fracture aperture; a slurry flow field estimation module for estimating the slurry diffusion flow field and spatiotemporal viscosity distribution based on a three-dimensional discrete fracture network model and the current fracture aperture, combined with a slurry spatiotemporal viscosity evolution model, to obtain the overall slurry pressure distribution, slurry velocity distribution, and real-time viscosity field; and a fracture deformation calculation module for applying the overall slurry pressure distribution as a normal load to the fracture wall boundary and calculating the dynamic deformation of the fracture under fluid-structure interaction to obtain the fracture aperture deformation; and particle deposition. The deposition analysis module determines the volume fraction of deposited particles based on slurry velocity distribution, real-time viscosity field, and fracture aperture deformation. The permeability update module updates the fracture network permeability and effective aperture based on the volume fraction of deposited particles to obtain the evolved fracture network permeability and corrected fracture aperture, with the corrected fracture aperture serving as the updated current fracture aperture. The plugging judgment module determines whether the evolved fracture network permeability meets a preset plugging threshold. If not, it feeds back the corrected fracture aperture to the slurry flow field estimation module to trigger iterative loops in each module until the plugging threshold is met to determine the optimal grouting duration.

[0006] Compared with existing technologies, the grouting control method provided in this application, which couples fracture aperture changes with grout particle deposition, firstly, effectively overcomes the calculation errors caused by traditional models. It can capture in real time the dynamic competitive relationship between fracture expansion caused by grout pressure and channel blockage caused by particle deposition during the grouting process, thereby accurately simulating the nonlinear evolution trajectory of fracture network permeability. Secondly, by introducing a spatiotemporal viscosity model, it solves the problem that a single temporal viscosity model cannot reflect spatial flow field differences, improving the fidelity of flow field simulation. Finally, it can accurately determine the critical state of grouting and sealing based on the convergence of evolving permeability, thereby determining the optimal grouting duration. This ensures the sealing effect while effectively avoiding waste of grout materials and potential formation disturbances caused by over-grouting, significantly improving the level of refined control in grouting projects. Attached Figure Description

[0007] The above and other objects, features and advantages of this application will become more apparent from the more detailed description of the embodiments of this application in conjunction with the accompanying drawings.

[0008] Figure 1 This is a flowchart of a grouting control method for coupling fracture aperture change and grout particle deposition according to an embodiment of this application.

[0009] Figure 2 This is a schematic diagram of the data flow of the grouting control method for coupling fracture opening change and grout particle deposition according to an embodiment of this application.

[0010] Figure 3 This is a schematic diagram illustrating the construction process of a three-dimensional discrete fracture network model in the grouting control method for coupled fracture aperture change and grout particle deposition according to an embodiment of this application.

[0011] Figure 4 This is a schematic diagram illustrating the construction process of the spatiotemporal viscosity evolution model of grout in the grouting control method for coupled fracture aperture change and grout particle deposition according to an embodiment of this application.

[0012] Figure 5 This is a block diagram of a grouting control system for coupled fracture opening variation and grout particle deposition according to an embodiment of this application.

[0013] Figure 6 This is a schematic diagram illustrating the principle of the spatiotemporal distribution of slurry viscosity and the deformation of cracks and pores according to an embodiment of this application.

[0014] Figure 7 This is a schematic diagram illustrating the mechanism by which particle deposition leads to a decrease in the permeability of a fracture network according to an embodiment of this application.

[0015] Figure 8 This is a simulation result of the spatiotemporal distribution of slurry viscosity under specific working conditions according to the embodiments of this application.

[0016] Figure 9 This is a simulation result diagram of the spatial distribution of slurry viscosity and deposited particles under specific working conditions according to the embodiments of this application.

[0017] Figure 10 This is a simulation result diagram of the spatial distribution of the crack aperture change under specific working conditions according to the embodiments of this application.

[0018] Figure 11 This is a simulation result diagram showing the evolution of crack aperture change over time under specific working conditions according to an embodiment of this application.

[0019] Figure 12 This is a cumulative permeability decrease curve for different fracture pore sizes according to the embodiments of this application.

[0020] Figure 13 This is a cumulative decrease curve of permeability under different grouting pressures according to the embodiments of this application.

[0021] Figure 14 This is a cumulative decrease curve of permeability under different water-cement ratios according to the embodiments of this application.

[0022] Figure 15 This is a cumulative decrease curve of permeability under different fracture network densities according to the embodiments of this application.

[0023] Figure 16This is a statistical chart showing the optimal grouting time under different fracture geometry characteristics according to the embodiments of this application.

[0024] Figure 17 This is a statistical chart showing the optimal grouting time under different grout characteristics according to the embodiments of this application. Detailed Implementation

[0025] The embodiments of this application will now be described in more detail with reference to the accompanying drawings. It should be understood that the drawings and embodiments of this application are for illustrative purposes only and are not intended to limit the scope of protection of this application.

[0026] This application proposes a grouting control method that couples the changes in fracture aperture with the deposition of grout particles. Figure 1 This is a flowchart of a grouting control method for coupling fracture aperture variation and grout particle deposition according to an embodiment of this application. Figure 2 This is a schematic diagram of the data flow of a grouting control method for coupling fracture aperture variation and grout particle deposition according to an embodiment of this application. Figure 1 and Figure 2 As shown, the grouting control method for coupling fracture aperture change and grout particle deposition according to an embodiment of this application includes:

[0027] S1: Obtain the current fracture aperture. It should be understood that fracture aperture, as the most critical parameter controlling the hydraulic properties of fractured rock masses, directly determines the flow resistance, diffusion velocity, and pressure distribution pattern of the grout in the fracture network. It is also a decisive geometric factor affecting whether suspended particles deposit and become blocked. Therefore, accurately obtaining the fracture aperture at the current moment is a prerequisite for achieving high-fidelity grouting simulation. Of course, it should be known that the current fracture aperture is a physical quantity that dynamically evolves with the grouting process, encompassing both the geometric aperture at the initial moment of grouting and the effective aperture after correction by physical field coupling during the grouting process.

[0028] Specifically, step S1 is divided into two levels based on the different stages of the grouting simulation: initial construction and iterative update. In the initial stage of the grouting simulation, obtaining the current fracture aperture mainly relies on the statistical processing of geological exploration data, which is directly related to the construction of the three-dimensional discrete fracture network (DFN) model. The construction of the DFN model is a geometric reconstruction process based on probabilistic statistics, aiming to transform limited geological exploration data into virtual fracture entities in three-dimensional space.

[0029] In one specific implementation of this application, such as Figure 3As shown, the construction process of the three-dimensional discrete fracture network model includes: firstly, acquiring geological statistical data, including parameters such as fracture radius, aperture, azimuth, and fracture density. This data can be obtained through on-site geological surveys. Specific implementation methods include using borehole television imaging technology to scan deep rock fractures, or laying survey lines in outcrop areas of the rock mass for manual statistical analysis, thereby obtaining raw data such as fracture attitude (dip, dip angle), fracture width, trace length, and linear density.

[0030] Then, the mean and variance parameters of the fracture radius in the geological statistics data are processed based on the Monte Carlo method and the log-normal distribution function, and the fracture azimuth parameter in the geological statistics data is processed in combination with the Fisher distribution function module to obtain the fracture geometric attribute set. This attribute set is a set of parameters describing the spatial characteristics of each independent fracture disk, specifically including the equivalent radius of the fracture disk, the initial fracture aperture, and the direction vector describing the spatial attitude of the fracture.

[0031] Specifically, for the equivalent radius of the rift disk and initial fracture aperture The acquisition of these two scalar parameters relies on a combination of the Monte Carlo method and the log-normal distribution function. The Monte Carlo method is a computational unit based on random sampling theory; its core function is to generate a sequence of random numbers conforming to a specific probability distribution to simulate the randomness in geological formation. The log-normal distribution function is a mathematical computational unit that encapsulates the log-normal probability density function and its inverse cumulative distribution function, used to describe the long-tailed distribution law generally followed by fracture scale parameters (such as radius and aperture) in nature. In practical applications, the Monte Carlo method first generates a series of pseudo-random numbers uniformly distributed in the interval [0,1], representing sampling points in the probability space. Subsequently, these random numbers are input into the log-normal distribution function, which utilizes a preset mean... and standard deviation For example, the fracture aperture follows , The log-normal distribution is used to map abstract probability values ​​to specific physical parameter values ​​through inverse transform sampling techniques. For example, regarding the crack radius... The module is based on the following formula:

[0032]

[0033] The determined distribution pattern converts the input random numbers into specific radius values. This process ensures that the generated fracture network contains both a large number of average-sized and medium-sized fractures, as well as a small number of large-scale fractures that play a controlling role in seepage.

[0034] For the orientation vector describing the spatial orientation of the fracture, the system converts the dip angle and dip direction into spherical coordinates and uses the Fisher distribution function module to generate the fracture orientation vector to simulate the concentration and dispersion of the fracture orientation. The Fisher distribution function is:

[0035]

[0036] In the formula, The polar angle of the fracture. This is a concentration parameter.

[0037] After determining the geometric properties of individual fractures, the system uses fracture density parameters from geostatistical data. The total number of cracks in the simulated region is calculated using the Poisson distribution function. :

[0038]

[0039] In the formula, This means that it is generated exactly within a given simulation region. The probability of a crack; The fracture density represents the expected number of fractures per unit volume of rock mass. This parameter is directly derived from the statistical analysis of geological exploration data. Represents the total volume of the current numerical simulation region; Represents the specific number of cracks (takes a non-negative integer); is the base of the natural logarithm. Subsequently, the system instantiates in three-dimensional space... The model consists of several disk-shaped fractures. The inter-fracture and interconnection relationships between fractures are addressed, isolated and invalid fractures are removed, and the resulting network is assembled into a three-dimensional discrete fracture network model capable of fluid computation. The initial aperture value generated during this process is used as the time step. The current fracture aperture at that time.

[0040] Meanwhile, to support subsequent flow field calculations, the system also needs to construct a spatiotemporal viscosity evolution model for the slurry. This model aims to address the problem of traditional models neglecting spatial rheological differences in the slurry; its key lies in establishing a mapping relationship between diffusion distance, residence time, and viscosity. For example... Figure 4 As shown, the construction process of the spatiotemporal viscosity evolution model of slurry includes: obtaining the initial rheological parameters of slurry through laboratory slurry rheological experiments (such as using a rotational viscometer), including the initial viscosity at different water-cement ratios. and yield stress and the rate constant of viscosity increase with time The system incorporates the diffusion distance of slurry particles in the fracture network. Combined with fracture porosity diffusion cross-sectional area and grouting flow rate Calculate the equivalent residence time of the grout as it travels from the grouting hole to the current location. Based on this, the purely time-dependent viscosity formula is transformed into a space-time coupled function:

[0041]

[0042] In the formula, These are rheological model constants. It's worth noting that in the above model... , , , These are all fixed parameters that can be measured in advance through a testing process. The model defines any time... and any position The viscosity of the slurry at a certain point was measured, enabling a refined characterization of the slurry's rheological properties.

[0043] As the grouting process progresses, it enters the time step. After the iterative calculation phase, obtaining the current fracture aperture becomes a feedback extraction of the calculation results from the previous time step. At this point, the fracture aperture is no longer merely a static geometric property, but a corrected value that incorporates stress deformation and particle deposition effects. Specifically, the system receives the updated current fracture aperture calculated in the subsequent step (i.e., step S5) and transmits it as input data to step S2 to update the solution domain boundary of the fluid dynamics equations, thereby forming a closed-loop dynamic simulation.

[0044] S2: Based on a three-dimensional discrete fracture network model and the current fracture opening, and combined with a slurry spatiotemporal viscosity evolution model, the slurry diffusion flow field and spatiotemporal viscosity distribution are estimated to obtain the overall slurry pressure distribution, slurry velocity distribution, and real-time viscosity field. It should be understood that grouting slurries are typically non-Newtonian fluids (such as Bingham fluids), and their rheological properties exhibit significant time-varying characteristics. Furthermore, in complex fracture networks, the viscosity of slurry at different locations exhibits high spatial non-uniformity due to differences in transport paths and residence times. Using only a single time-dependent viscosity model cannot accurately reflect the velocity reduction and pressure accumulation phenomena caused by increased viscosity at the far end of the fracture. Therefore, through coupled iteration of the flow field and viscosity field, the hydrodynamic state across the entire field at each moment is accurately calculated, providing an accurate mechanical and kinematic basis for subsequent fracture deformation calculations and particle deposition analysis.

[0045] In one specific implementation of this application, step S2 includes: First, using the initial three-dimensional discrete fracture network model as the solution domain, the current fracture aperture and the spatiotemporal viscosity evolution model of the slurry are introduced into the Navier-Stokes equations and the continuity equation to obtain the fluid dynamics governing equations. That is, using the three-dimensional discrete fracture network model obtained in step S1 as the geometric basis, the physical boundary gap of the fluid flow is defined using the current fracture aperture. Within this solution domain, the system treats the slurry as an incompressible fluid and establishes a set of fluid dynamics governing equations including spatiotemporal viscosity terms. Specifically, a dynamic viscosity function is introduced into the dissipation term of the Navier-Stokes equations and described in conjunction with the continuity equation. The mathematical expression of its momentum conservation equation is:

[0046]

[0047] Simultaneously satisfies the continuity equation:

[0048]

[0049] In the formula, Indicates the density of the slurry. Represents the slurry flow velocity vector. Indicates time, Indicates slurry pressure, Represents volume force. This indicates dependence on dwell time. The spatiotemporal distribution viscosity.

[0050] After constructing the governing equations, the system enters the core numerical solution stage, which involves using a finite element method (FEM) and pressure inlet boundary conditions to perform coupled iterative solutions of the fluid dynamics governing equations for velocity and pressure. The specific implementation of this process is as follows: The system first uses the finite element method (FEM) to spatially discretize the aforementioned partial differential equations, transforming the continuous physical field into a system of algebraic equations with a finite number of nodes. Regarding the boundary conditions, the pressure at the grouting hole is set as a Dirichlet boundary condition (i.e., grouting pressure). The fracture wall is set as a no-slip boundary condition (i.e., the flow velocity at the wall is zero). This is because the viscosity field $\mu(t_\mu)$ and the velocity field... There is a strong coupling relationship: the flow velocity determines the residence time of the slurry at a certain point, thus determining the viscosity, and the viscosity, in turn, determines the flow resistance and flow velocity, making a direct solution impossible in one step. The system employs a pressure-velocity coupled algorithm (such as PISO or SIMPLE) for iterative calculation: first, based on the viscosity field of the previous time step or an initial guess, the momentum equation is solved to obtain the predicted velocity field; then, a pressure correction equation is constructed using the continuity equation to correct the pressure and velocity fields to satisfy mass conservation; subsequently, streamline tracing is performed using the corrected velocity distribution to calculate the velocity at each spatial location. The cumulative diffusion path length of grout particles after leaving the injection hole, combined with the fracture porosity. Flow area and traffic Convert path length into effective dwell time Next, the specific viscosity value at that point is calculated using the slurry spatiotemporal viscosity evolution model:

[0051]

[0052] In the formula, The initial viscosity, The viscosity growth rate constant is related to the water-cement ratio. The rheological constant is used. The calculated new viscosity value is mapped back to the grid to generate an updated real-time viscosity field, which is then substituted into the momentum equation for the next round of solving. This process is repeated until the residuals of the pressure field, velocity field, and viscosity field are all less than the preset convergence criterion. Finally, the full-field slurry pressure distribution, slurry velocity distribution, and real-time viscosity field are output.

[0053] S3: The entire field grout pressure distribution is applied as a normal load to the fracture wall boundary, and the dynamic deformation of the fracture under fluid-structure interaction is calculated to obtain the fracture aperture deformation. It should be understood that in a high-pressure grouting environment, the fractured rock mass is not an absolutely rigid medium; the grout pressure directly acts on the fracture wall, causing the fracture to open or close elastically. This deformation significantly changes the geometric aperture of the fracture, which in turn affects the flow resistance and diffusion range of the grout. Ignoring this dynamic deformation effect will lead to significant errors in the prediction of grout volume and diffusion distance. Therefore, it is necessary to accurately calculate the rock mass displacement caused by fluid pressure and quantify the dynamic changes in fracture aperture to achieve high-fidelity grouting simulation.

[0054] In one specific implementation of this application, step S3 includes: first, identifying the contact interface between the grout and the fracture wall, and extracting the nodal pressure values ​​from the full-field grout pressure distribution calculated in step S2. Based on the assumptions and magnitude analysis of the Goodman model, the slight influence of fluid shear stress on the opening deformation is ignored, and the fluid pressure is converted into normal stress acting on the rock mass surface. According to the formula... Construct fluid-solid interface stress loads, where For the rock mass surface stress tensor, This is the unit normal vector of the fracture wall. The load is directed along the normal to the fracture wall and its magnitude is equal to the fluid pressure, forming the boundary condition for the solid mechanics calculation. Subsequently, the system treats the rock mass as a linear elastic material and calls the rock density... Elastic modulus Compared to Poisson By applying the aforementioned fluid-solid interface stress load to the solid mechanics governing equations and assuming equal rock mechanics parameters, a rock mass deformation governing equation based on momentum conservation is established:

[0055]

[0056] In the formula, The rock mass displacement vector. The force is a volume force. The equation is solved using a finite element method (FEM) to calculate the full-field displacement response of the rock mass under fluid pressure, generating a full-field rock mass displacement field. Finally, based on the full-field rock mass displacement field, the fracture aperture deformation is determined.

[0057] Regarding how to determine the fracture aperture deformation based on the whole-field rock mass displacement field, this application provides two specific implementation methods, which are applicable to different geological scenarios with different accuracy requirements.

[0058] In the first implementation, the system employs the assumption of symmetrical deformation or the principle of superposition of displacements on both sides of the fracture. The system extracts the displacement vector data at the fracture wall. Since the rock masses on both sides of the fracture undergo opposing separation displacements under internal pressure, the system calculates the projected component of the wall displacement in the normal direction and directly multiplies it by a coefficient of 2 to represent the total aperture change caused by the symmetrical displacement of the fracture walls. The calculation formula is as follows:

[0059]

[0060] In the formula, This represents the deformation amount of the crack aperture.

[0061] However, research has shown that in the discrete fracture network (DFN) grouting simulation technical scheme, the calculation accuracy of fracture aperture deformation is relatively low due to idealized assumptions, directly affecting the accuracy of grouting and sealing prediction under complex geological conditions. Specifically, the calculation logic in the first implementation method adopts a simplified processing based on the assumption of symmetrical deformation or the superposition of displacements on both sides of the wall, directly using formulas... Solving for the opening increment has a double defect in terms of physical reality.

[0062] First, this implementation ignores the stiffness asymmetry caused by the fracture network topology. In a real DFN scenario, fractures act like random cutting lines, slicing the rock mass into blocks of varying shapes and volumes. One side of a fracture may be surrounded by multiple intersecting fractures, forming isolated, weakly constrained, small critical blocks with low stiffness and high degrees of freedom; while the other side may be an uncut, continuous, semi-infinite rock mass with high stiffness and strong constraints. Under the same grouting pressure, the displacement responses of the two rock walls are drastically different. The original solution simply multiplies the displacement on one side by a coefficient of 2, implicitly assuming that the two rock masses are symmetrical half-spaces, completely losing the anisotropy of the mechanical response determined by the topological cutting method. This results in the prediction of aperture changes near fracture intersections only applicable to idealized scenarios.

[0063] Secondly, this implementation method overlooks the shear dilatation effect induced by fluid pressure. The original formula only considers elastic expansion caused by normal stress unloading, which is consistent with the characteristics of smooth fractures. However, fracture surfaces in actual engineering are usually rough and undulating. During high-pressure grouting, the grout pressure significantly counteracts the effective normal stress acting on the fracture surface, leading to a sharp decrease in the fracture's shear strength and inducing minute shear displacement on both sides of the rock wall. For rough fracture surfaces, this displacement is inevitably accompanied by a climbing effect, resulting in irreversible volume expansion, i.e., shear dilatation. The original mechanism lacks this plastic increment component, leading to an underestimation of the actual conductivity of rough fractures under high pressure.

[0064] To address the aforementioned shortcomings, the second implementation proposes a refined computational process based on topological and constitutive coupling, aiming to reconstruct the physical reality of fracture aperture evolution.

[0065] Specifically, in the second implementation, the fracture aperture deformation is determined based on the full-field rock mass displacement field. This includes: First, generating a topological asymmetry correction factor based on discrete fracture network topology data and the full-field rock mass displacement field. Addressing the issue of the first implementation failing to perceive fracture asymmetry constraints, this step first introduces discrete fracture network topology data and the full-field rock mass displacement field as inputs, aiming to quantify the stiffness mismatch caused by the geometric differences in cutting on both sides of the fracture. By extracting the volumetric characteristics of local rock blocks and the distance information from nodes to the nearest intersection zone, a mathematical model capable of dynamically correcting the symmetry assumption is constructed. The core of the calculation lies in constructing the topological asymmetry correction factor. The physical essence of this factor is a dynamic scaling of coefficient 2 based on the specific context.

[0066] The specific calculation formula is as follows:

[0067]

[0068] in, This is the correction factor for the output; and These represent the equivalent volumes of local rock blocks on the positive and negative sides of the fracture, respectively. The normalized ratio of their difference directly reflects the degree of difference in local constraint stiffness. The volume difference sensitivity coefficient is used to adjust the influence weight of volume mismatch on stiffness. It can be obtained by fitting numerical experiments under different rock block size ratios. For example, when one side of the fracture is an isolated small rock block cut by multiple fractures and the other side is intact bedrock, the value can be 0.5. is the Euclidean distance from the current computation node to the nearest crack intersection line, and The preset characteristic length of rock mass damage is determined based on the average width of the fracture zone observed in the field geological survey; the two constitute an exponential term. The physical phenomenon was simulated that the local compliance of the rock mass increases exponentially as it approaches the cross fracture zone. The enhancement factor is determined based on the elastic modulus reduction rate of the rock mass in the cross zone during rock mechanics testing.

[0069] This step gives the numerical model the ability to perceive the geometry of the DFN: near isolated rock blocks or intersections, the calculated... A value significantly greater than 1 means that the model can predict larger opening deformations than the traditional symmetry assumption, thereby accurately capturing the instability characteristics of key blocks under high pressure, realizing physical correction of deformation prediction in areas with weakened local stiffness, and achieving the goal of improving the accuracy of opening calculation under complex topological structures.

[0070] Then, based on the overall grout pressure distribution, fracture surface roughness parameter set, and initial in-situ stress parameters, the shear dilatation aperture component is determined. Given the lack of a plastic dilatation component in the first implementation method, this step introduces the overall grout pressure distribution, fracture surface roughness parameter set (including JRC and JCS), and initial in-situ stress parameters. Based on Barton-Bandis rock mechanics theory, the shear dilatation aperture component induced by the reduction in effective stress is calculated. The core logic of this step lies in quantifying the grouting pressure. Ground stress The offsetting effect is used to deduce the volume increase of the rough surface climbing caused by the reduction in shear strength.

[0071] The calculation formula is as follows:

[0072]

[0073] In the formula, This is the obtained shear expansion component; This represents the local micro-shear displacement that occurs during the grouting process. This value is usually set to one-thousandth of the characteristic size of the crack or determined by a shear test. The joint roughness coefficient is obtained through standard comparison of the fracture surface profile in the field or 3D scanning; for example, a value of 10 is used. The compressive strength of the joint wall was obtained through on-site testing using a Schmidt rebound hammer; for example, a value of 50 MPa was taken. The initial normal geostress was obtained through field measurements using the hydraulic fracturing method or the stress relief method. This represents the current slurry fluid pressure. The logarithmic term in the formula... Dynamically describes the change in pore pressure The increase in shear dilatation leads to a decrease in the effective normal stress at the fracture surface, resulting in a nonlinear mechanical process that increases the shear dilatation angle.

[0074] This step represents a leap from simple elastic mechanics to elastoplastic coupling mechanics. Its physical essence explains why high-pressure grouting often leads to the permanent expansion of fracture channels in hard but rough fractures. By introducing a roughness parameter, this step enables the model to distinguish the differences in grouting response between smooth and rough fractures, correcting the underestimation of opening caused by single elastic predictions. This results in a more realistic reflection of the physical process of high-pressure grouting opening rough fractures.

[0075] Finally, a non-uniform stiffness correction is performed on the entire rock mass displacement field using a topological asymmetric correction factor, and a multi-mechanism deformation superposition is performed using the shear expansion aperture component to obtain the fracture aperture deformation. As the final comprehensive processing step, this step fuses the topological asymmetric correction factor and shear expansion aperture component generated in the previous steps with the entire rock mass displacement field. The final fracture aperture deformation is synthesized using a nonlinear superposition algorithm. This calculation no longer simply relies on the doubling of normal displacement, but instead utilizes… The unilateral elastic displacement is weighted and amplified, and the shear dilatation increment is explicitly superimposed, as shown below:

[0076]

[0077] in, Crack aperture deformation; Based on the elastic displacement of the fracture wall in the normal direction; The term implements a weighted correction to topological asymmetry (when...) (Time-regression symmetry assumption) Shear expansion opening component.

[0078] This step elevates the fluid-solid coupling to a fluid-solid-structure-morphology coupling, ensuring that the final output aperture data not only responds to the macroscopic fluid pressure field but also incorporates the unique geometric cutting characteristics and microscopic surface morphology information of the fracture network. This provides an extremely accurate geometric boundary condition that conforms to the true geomechanical values ​​for subsequent permeability updates, thereby significantly improving the adaptability and predictive reliability of the entire grouting control model to complex geological environments.

[0079] The second approach, through the implementation of the aforementioned technical means, fundamentally solves the problem of insufficient response to the complexity of geological structures in the original scheme during DFN grouting simulation. By introducing topological asymmetry correction, the mechanism can accurately locate and quantify the local stiffness weakening region caused by rock block cutting, eliminating the prediction deviation of the cross-zone aperture caused by the symmetry assumption. By coupling a shear expansion model, the mechanism successfully captures the plastic volume expansion of rough fractures under high pressure, correcting the underestimation of grouting volume by the traditional elastic model. This series of improvements enables the numerical model to output highly accurate fracture aperture evolution data even under extreme geological conditions (such as fracture zones and high geostress areas), significantly improving the scientific rigor and accuracy of grout diffusion range and sealing efficiency predictions.

[0080] S4: Determine the volume fraction of deposited particles based on slurry velocity distribution, real-time viscosity field, and fracture aperture deformation. It should be understood that slurry flow in fractures is not a simple fluid transport process, but rather involves the dynamic deposition of suspended particles onto the fracture walls or throats. This deposition behavior is jointly controlled by flow velocity, viscosity, and fracture geometry. Deposited particles occupy effective pore space, directly leading to a decrease in permeability. Therefore, only by comprehensively considering fluid dynamics parameters and geometric deformation parameters and accurately calculating the volume fraction of deposited particles can the true clogging effect and permeability evolution during the grouting process be reflected.

[0081] In one specific implementation of this application, step S4 includes: First, the system uses the current crack aperture obtained in step S1. The crack aperture deformation amount calculated in step S3 Adding them together gives the instantaneous dynamic fracture aperture through which the fluid actually passes. .

[0082] Subsequently, based on the pore media filtration theory, a deposition coefficient calculation model was constructed. This model reflects the physical mechanisms that higher flow velocities introduce more particles, smaller apertures facilitate interception, and higher viscosity provides stronger confinement. Based on the overall slurry flow velocity distribution, real-time viscosity field, and instantaneous dynamic aperture, the dynamic deposition coefficient for the entire field was calculated point-by-point. The calculation formula is as follows:

[0083]

[0084] In the formula, The current fracture porosity, The velocity modulus of the slurry. The viscosity is the real-time spatiotemporal viscosity.

[0085] Next, a mass conservation equation for the slurry particles is established, including convection and transport terms and sedimentation sink terms. Assuming the transformation from suspended particles to sedimentary particles is unidirectional and irreversible, the calculated dynamic sedimentation coefficient is substituted into the equation as the source term coefficient:

[0086]

[0087] Simultaneously solve the sedimentary facies evolution equation:

[0088]

[0089] In the formula, The mass concentration of slurry particles. This represents the volume fraction of suspended particles. This represents the volume fraction of the deposited particles. For fluid density, For the density of water, Let be the density of the sedimentary rock mass. By solving the above set of partial differential equations, the system updates the volume fraction of suspended particles at the current moment and calculates the cumulative volume fraction of sedimentary particles.

[0090] S5: The permeability and effective aperture of the fracture network are updated based on the volume fraction of deposited particles to obtain the evolved fracture network permeability and the corrected fracture aperture, which serves as the updated current fracture aperture. It should be understood that as the grouting process continues, solid particles deposited on the fracture walls and in the pores gradually occupy the flow space originally belonging to the fluid, leading to a decrease in effective porosity and consequently a nonlinear decay in permeability. If these hydraulic parameters are not updated in a timely manner, subsequent flow field calculations will be unable to detect the resistance changes caused by blockage, resulting in distorted simulation results. Therefore, it is necessary to dynamically correct the permeability and aperture parameters based on the accumulated amount of deposited particles and feed the corrected aperture back to the flow field calculation stage to form a realistic physical coupling cycle.

[0091] In one specific implementation of this application, step S5 includes: first, based on the porous medium volume-filling model, using the volume fraction of deposited particles calculated in step S4. Calculate the volume occupied by the sediment. Starting from the initial porosity... The corrected effective porosity is obtained by subtracting the normalized sedimentary particle volume fraction from the value. The calculation formula is as follows:

[0092]

[0093] In the formula, For fluid density, This represents the density of the rock matrix.

[0094] Subsequently, the system introduces the Kozeny-Carman (KC) equation to establish a nonlinear mapping relationship between permeability and porosity. The effective porosity is then... Volume fraction of sedimentary particles and suspended particulate volume fraction Substituting into the modified KC equation, the evolutionary fracture network permeability reflecting the current clogging state is calculated. :

[0095]

[0096] In the formula, The initial permeability is given. This formula reflects the cubic sensitivity of permeability to changes in porosity.

[0097] Finally, the system utilizes the inverse relationship of the Cubic Law for flow in flat plate fractures to calculate the evolved permeability. Converted to equivalent modified fracture aperture The calculation formula is:

[0098]

[0099] In the formula, The length of the fractured rock mass element is given. The generated corrected fracture aperture is then assigned the new current fracture aperture.

[0100] In one specific embodiment, the initial permeability of a certain fracture element is The initial porosity was 0.1. After a period of grouting, the volume fraction of deposited particles reached 0.5. The system first calculated that the effective porosity decreased to 0.05. Substituting into the KC equation, the calculated evolutionary permeability decreased sharply to... Subsequently, the system uses the cubic law to deduce that the corresponding corrected fracture aperture decreases from the initial 1 mm to 0.5 mm. This updated 0.5 mm aperture value will be transmitted back to the flow field calculation module to simulate higher flow resistance in the next moment, thus realistically reflecting the grouting and sealing process.

[0101] S6: Determine if the permeability of the evolving fracture network meets the preset sealing threshold. If not, correct the fracture aperture and feed it back to steps S2 through S5 for iterative cycling until the sealing threshold is met to determine the optimal grouting time. It should be understood that grouting projects need to find the optimal balance between ensuring sealing effectiveness and controlling project costs. If the grouting time is too short, the fracture network will not be effectively filled, and the permeability will not drop to a safe level, leading to project failure. If the grouting time is too long, it not only wastes grout material, but the continuous high pressure may also cause formation uplift or damage. Therefore, an automatic determination mechanism based on quantitative indicators is established to scientifically determine the grouting endpoint by monitoring the evolution of permeability.

[0102] In one specific implementation of this application, step S6 includes: first, receiving the evolutionary fracture network permeability output in step S5. And call the initial penetration rate. Calculate the current penetration rate reduction rate. This ratio directly reflects the sealing effectiveness produced by grouting, and its calculation formula is as follows:

[0103]

[0104] Subsequently, the system will calculate the results. With preset blocking threshold In comparison, the preset sealing threshold can be set according to project requirements. For example, a grouting project may require reducing rock permeability by more than 98%. Below This indicates that the current grouting volume has not yet reached the sealing requirements, and the system determines that the grouting process is incomplete, thus triggering the feedback mechanism. At this point, the system increases the current grouting time by one time step. The corrected fracture aperture generated in step S5 is assigned to the current fracture aperture, and used as a new geometric boundary condition to be fed back to step S2. This feedback action restarts the computational chain from step S2 to step S5: flow field calculation - deformation solution - deposition analysis - parameter update, forming a closed-loop iterative cycle. Conversely, if Higher than or equal to If the permeability decline curve enters a plateau phase (the slope approaches zero), the system determines that the grouting has achieved the expected effect, immediately terminates the iterative calculation, locks the current grouting time, and marks it as the optimal grouting duration.

[0105] To more intuitively understand the physical mechanism of the grouting control method for coupling fracture aperture change and grout particle deposition described in this application, please refer to [reference needed]. Figure 6 and Figure 7 .like Figure 6As shown, this application fully considers that the viscosity of the slurry is not constant during diffusion within a fracture network, but rather exhibits spatiotemporal distribution characteristics depending on the residence time (aging time) and spatial location of the slurry particles in the flow channel (the shades of color in the figure represent viscosity). Simultaneously, the slurry pressure acting on the fracture wall causes dynamic deformation of the fracture pores, which in turn affects the slurry flow path. For example... Figure 7 As shown, as grouting progresses, suspended particles gradually transform into deposited particles, occupying the original fracture space and narrowing the effective flow channels, thus causing a nonlinear decrease in the permeability of the fracture network. This application achieves precise control of the grouting process based on coupled calculations of the above mechanism.

[0106] Furthermore, in order to verify the effectiveness of the grouting control method for coupled fracture aperture change and slurry particle deposition in the embodiments of this application, a numerical simulation experiment was conducted on a typical three-dimensional fracture network based on the grouting control method for coupled fracture aperture change and slurry particle deposition described in this application.

[0107] Figures 8 to 11 This demonstrates the evolution of the grouting process under a specific working condition (e.g., a water-cement ratio of 1.0 and a grouting pressure of 0.5 MPa). Figure 8 As shown, the grout viscosity is lowest near the grouting pipe and gradually increases with increasing diffusion distance and time, exhibiting significant spatiotemporal heterogeneity. Figure 9 The spatial correspondence between grout viscosity and deposited particles was further demonstrated at different times during grouting (e.g., 12 min and 18 min). It can be seen that high viscosity areas are often accompanied by more significant particle deposition.

[0108] Figure 10 and Figure 11 The dynamic response law of fracture aperture was revealed. For example... Figure 10 As shown, the change in crack aperture occurs earlier and with a greater magnitude near the grouting hole, which is consistent with the spatial decay law of grout pressure. Figure 11 As shown, with the passage of grouting time, the aperture variation exhibits a trend of first increasing rapidly and then gradually slowing down. This indicates that in the initial stage of grouting, grout pressure dominates the elastic opening of the fracture; while in the later stage, with the intensification of particle deposition and pressure equilibrium, the aperture variation tends to stabilize. These simulation results fully demonstrate that the grouting control method of coupling fracture aperture variation and grout particle deposition in the embodiments of this application can accurately capture the microscopic evolution characteristics under the coupling of the "fluid-solid-particle" three fields.

[0109] To further verify the sensitivity of the proposed method in predicting the sealing effect under different operating conditions, Figures 12 to 15 The curve showing the cumulative decrease rate of permeability as a function of slurry diffusion distance, calculated based on step S5 of this application, is presented.

[0110] like Figure 12 As shown, the effect of fracture pore size on the sealing effect is illustrated. It can be seen that fracture networks with larger pore sizes (10.1-11.0 mm) can ultimately achieve a higher cumulative permeability reduction rate (approximately 90%), while fractures with smaller pore sizes (0.1-1.0 mm) have relatively lower final sealing efficiency due to high flow resistance and limited slurry diffusion.

[0111] like Figure 13 As shown, the effect of grouting pressure is illustrated. High grouting pressure (6.0-7.5 MPa) provides a stronger driving force, enabling the grout to overcome resistance and penetrate deeper into the fracture branches. Therefore, its permeability decreases the fastest and the final decrease rate is the highest (over 85%). In contrast, the sealing effect of low-pressure grouting (0.5-1.5 MPa) is significantly limited.

[0112] like Figure 14 As shown, the effect of water-cement ratio is illustrated. High water-cement ratio slurries (1.6-1.8) have good fluidity and can diffuse more widely to fill micro-cracks, thus exhibiting the highest permeability reduction rate; while low water-cement ratio slurries (0.8-1.0) have high viscosity and limited diffusion range, resulting in a lower overall sealing rate.

[0113] like Figure 15 The figure shows the effect of fracture network density. A high-density fracture network (6.6-7.0 m) −-1 It provides more connectivity paths, which is conducive to the uniform diffusion and deposition of slurry, thus achieving a better sealing effect.

[0114] The above Figures 12 to 15 The simulation results not only quantify the influence mechanism of each parameter on the grouting effect, but more importantly, they verify that the control method proposed in this application can accurately output the evolved permeability index through iterative calculation based on real-time working conditions, providing solid data support for determining whether the preset sealing threshold has been reached in step S6.

[0115] Finally, the reliability of the optimal grouting time determined by the grouting control method based on the coupling fracture aperture change and grout particle deposition in the embodiments of this application is verified. Figure 16 and Figure 17 It has been verified. Figure 16The optimal grouting time was statistically analyzed under different fracture geometry characteristics (fracture pore size, fracture network density). The results showed that when the fracture pore size was 5-15 mm and the network density was high, the optimal grouting time was relatively short (about 8-13 minutes). Figure 17 The optimal grouting time under different grout characteristics (grouting pressure, water-cement ratio) was statistically analyzed, showing that high pressure and high water-cement ratio are beneficial to shortening the grouting time. These statistical patterns are consistent with engineering experience, further confirming that the grouting control method proposed in this application can scientifically and quantitatively determine the optimal grouting time according to different geological conditions and construction parameters, thereby ensuring the sealing effect while avoiding material waste.

[0116] In summary, the grouting control method for coupled fracture aperture change and slurry particle deposition in this application has been clarified. This method abandons the traditional static fracture assumption and constructs a closed-loop dynamic iterative calculation system: First, based on a three-dimensional discrete fracture network model, a slurry spatiotemporal viscosity evolution model is introduced to finely characterize the rheological properties of the slurry at different transport paths and time nodes; then, the fluid-structure interaction mechanism is used to convert the slurry pressure across the entire field into a normal load on the fracture wall, and the dynamic opening deformation of the fracture under fluid pressure is solved in real time; on this basis, combined with the velocity distribution, real-time viscosity field and the geometric characteristics after deformation, the deposition behavior of suspended particles in the fracture is quantified, and the permeability and effective aperture of the fracture network are dynamically corrected accordingly; then, the corrected fracture aperture is fed back to the flow field calculation stage, forming an iterative cycle of flow field calculation—deformation solution—deposition analysis—parameter update, until the evolved permeability meets the preset sealing threshold. By implementing the above technical solutions, this application can accurately reflect the evolution mechanism of multi-physics field coupling in complex geological environments, precisely capture the dynamic competitive relationship between grout pressure-induced fracture expansion and particle deposition-induced channel blockage during grouting, and accurately simulate the nonlinear evolution trajectory of fracture network permeability. At the same time, by using a spatiotemporal viscosity model, it overcomes the shortcomings of traditional single time viscosity models in failing to reflect spatial flow field differences, significantly improving the fidelity of flow field simulation. Finally, based on the convergence state of the evolving permeability, it scientifically determines the critical point of grouting and sealing, accurately determines the optimal grouting time, and effectively avoids grout waste and formation disturbance that may be caused by excessive grouting while ensuring the sealing effect of the project, thereby greatly improving the level of fine control and construction efficiency of grouting projects.

[0117] Figure 5 This is a block diagram of a grouting control system for coupled fracture aperture variation and grout particle deposition according to an embodiment of this application. Figure 5As shown, the grouting control system 500 for coupling fracture aperture change and slurry particle deposition according to an embodiment of this application includes: a fracture aperture acquisition module 510 for acquiring the current fracture aperture; a slurry flow field estimation module 520 for estimating the slurry diffusion flow field and spatiotemporal viscosity distribution based on a three-dimensional discrete fracture network model and the current fracture aperture, combined with a slurry spatiotemporal viscosity evolution model, to obtain the overall slurry pressure distribution, slurry velocity distribution, and real-time viscosity field; and a fracture deformation calculation module 530 for applying the overall slurry pressure distribution as a normal load to the fracture wall boundary and calculating the dynamic deformation of the fracture under fluid-structure interaction to obtain the fracture aperture deformation. The particle deposition analysis module 540 is used to determine the volume fraction of deposited particles based on the slurry flow velocity distribution, real-time viscosity field, and fracture aperture deformation. The permeability update module 550 is used to update the fracture network permeability and effective aperture according to the volume fraction of deposited particles to obtain the evolved fracture network permeability and the corrected fracture aperture, wherein the corrected fracture aperture is used as the updated current fracture aperture. The plugging judgment module 560 is used to determine whether the evolved fracture network permeability meets the preset plugging threshold. If it does not meet the threshold, the corrected fracture aperture is fed back to the slurry flow field estimation module to trigger each module to perform iterative loops until the plugging threshold is met to determine the optimal grouting time.

[0118] It should be noted that the grouting control system for coupled fracture opening change and slurry particle deposition in this application embodiment is similar in principle to the aforementioned grouting control method for coupled fracture opening change and slurry particle deposition. Therefore, the implementation process, implementation principle, and beneficial effects of the grouting control system for coupled fracture opening change and slurry particle deposition can all be found in the description of the implementation process, implementation principle, and beneficial effects of the aforementioned method, and will not be repeated.

[0119] The various embodiments of this disclosure have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical application, or improvement of the technology in the market, or to enable others skilled in the art to understand the embodiments disclosed herein.

Claims

1. A grouting control method that couples fracture aperture variation with grout particle deposition, characterized in that, include: S1: Get the current crack aperture; S2: Based on the three-dimensional discrete fracture network model and the current fracture opening, and combined with the slurry spatiotemporal viscosity evolution model, the slurry diffusion flow field and spatiotemporal viscosity distribution are estimated to obtain the full-field slurry pressure distribution, slurry velocity distribution and real-time viscosity field; S3: Apply the full-field slurry pressure distribution as a normal load to the fracture wall boundary and calculate the dynamic deformation of the fracture under fluid-structure interaction to obtain the fracture aperture deformation. This includes: extracting the nodal pressure values ​​in the full-field slurry pressure distribution and converting them into compressive stress acting on the fluid-structure interaction interface based on the fracture wall normal vector to obtain the fluid-structure interface stress load; substituting the fluid-structure interface stress load as a boundary condition into the linear elastic solid momentum conservation equation containing rock mechanics parameters for finite element solution to obtain the full-field rock mass displacement field; and determining the fracture aperture deformation based on the full-field rock mass displacement field. S4: Determine the volume fraction of deposited particles based on slurry flow velocity distribution, real-time viscosity field, and fracture aperture deformation. S5: The fracture network permeability and effective aperture are updated based on the volume fraction of deposited particles to obtain the evolved fracture network permeability and the corrected fracture aperture, which is used as the updated current fracture aperture. S6: Determine whether the permeability of the evolved fracture network meets the preset sealing threshold. If it does not meet the threshold, the fracture opening is corrected and fed back to steps S2 to S5 for iterative looping until the sealing threshold is met to determine the optimal grouting time. Among them, the deformation of fracture aperture is determined based on the displacement field of the entire rock mass, including: Based on discrete fracture network topology data and the full-field rock mass displacement field, a topological asymmetric correction factor is generated, and the calculation formula is as follows: in, As a correction factor; and These represent the equivalent volumes of local rock blocks on the positive and negative sides of the fracture, respectively. This is the sensitivity coefficient for volume difference; is the Euclidean distance from the current computation node to the nearest crack intersection line, and Preset rock mass damage characteristic length, For enhancement coefficient; Based on the overall slurry pressure distribution, fracture surface roughness parameter set, and initial in-situ stress parameters, the shear expansion aperture component is determined using the following formula: In the formula, This is the shear expansion component; This represents the local micro-shear displacement that occurs during the grouting process; For joint roughness coefficient, For the compressive strength of the joint wall, For the initial normal geostress, The current slurry fluid pressure; The non-uniform stiffness of the entire rock mass displacement field is corrected using a topological asymmetric correction factor, and the fracture aperture deformation is obtained by superimposing multiple deformation mechanisms in combination with the shear expansion aperture component. The calculation formula is as follows: in, This represents the deformation amount of the crack aperture. Based on the normal elastic displacement of the fracture wall, As a correction factor; This is the shear expansion opening component.

2. The grouting control method for coupling fracture aperture change and grout particle deposition according to claim 1, characterized in that, The construction process of the three-dimensional discrete fracture network model and the slurry spatiotemporal viscosity evolution model includes: Obtain geological statistics and initial rheological parameters of the slurry; The mean and variance parameters of fracture radius in geological statistics are processed using the Monte Carlo method and the log-normal distribution function, and the fracture azimuth parameter in geological statistics is processed using the Fisher distribution function module to obtain the set of fracture geometric attributes. Based on the fracture density parameter in the geological statistics, the fracture entities are instantiated and the connectivity is assembled using the Poisson distribution function and the fracture geometric attribute set to obtain the three-dimensional discrete fracture network model. Based on the initial viscosity and the viscosity growth rate constant related to the water-cement ratio in the initial rheological parameters of the slurry, a spatiotemporal viscosity evolution model of the slurry is constructed.

3. The grouting control method for coupling fracture aperture change and grout particle deposition according to claim 1, characterized in that, Step S2 includes: Using the initial three-dimensional discrete fracture network model as the solution domain, the current fracture aperture and slurry spatiotemporal viscosity evolution model are introduced into the Navier-Stokes equations and the continuity equation to obtain the fluid dynamics control equation set. Using a finite element solver and pressure inlet boundary conditions, the fluid dynamics control equations are solved iteratively with velocity and pressure to obtain the full field slurry pressure distribution and the full field slurry velocity distribution. The streamline tracing method is used to calculate the cumulative diffusion path and effective residence time of slurry particles by utilizing the full-field slurry velocity distribution, and the node mapping is performed by calling the slurry spatiotemporal viscosity evolution model to obtain the real-time viscosity field.

4. The grouting control method for coupling fracture aperture change and grout particle deposition according to claim 1, characterized in that, Step S4 includes: The instantaneous dynamic opening is obtained by adding the crack aperture deformation amount to the current crack aperture; Based on the full-field slurry velocity distribution, real-time viscosity field, and instantaneous dynamic opening, the dynamic deposition coefficient of the full-field distribution is determined; The dynamic deposition coefficient is substituted as the source term coefficient into the slurry mass conservation equation, which includes convection and deposition terms, to update the suspended particle volume fraction and calculate the cumulative sedimentary particle volume fraction.

5. The grouting control method for coupling fracture aperture change and grout particle deposition according to claim 1, characterized in that, Step S5 includes: Calculation of sediment-occupied volume based on porous media volume-filling model; The effective porosity is obtained by subtracting the normalized volume fraction of deposited particles from the initial porosity. Nonlinear mapping was performed on effective porosity, volume fraction of deposited particles, and volume fraction of suspended particles to obtain the evolutionary fracture network permeability that reflects the current clogging state. The corrected fracture aperture is obtained by inverse geometric derivation of the permeability of the evolved fracture network.

6. The grouting control method for coupling fracture aperture change and grout particle deposition according to claim 1, characterized in that, Step S6 includes: Calculate the rate of decrease in permeability of the evolved fracture network relative to the initial permeability, and compare this rate with a preset plugging threshold for convergence. If the reduction rate is lower than the preset blocking threshold, the corrected crack opening will be fed back to steps S2 to S5 for iterative looping until the blocking threshold is met. If the reduction rate is higher than or equal to the preset sealing threshold, the current grouting time is locked and marked as the optimal grouting duration.

Images