A method for constructing a seepage grouting model of rock mass with complex fracture network based on COMSOL

Abstract

The invention discloses a method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL, which relates to the technical field of rock mass mechanics and solves the problem that fracture samples obtained collapse due to the shortcomings of low sampling efficiency and large sample capacity, resulting in a large difference between the established model and the actual project. The rock mass permeation grouting model construction method comprises conducting geological exploration on the rock mass, and extracting and processing fracture geometric parameters of the rock mass; performing fitting correction on each fracture geometric parameter to obtain an optimal distribution probability function; generating a discrete fracture network model based on the optimal distribution probability function by using a Latin hypercube sampling algorithm; converting the discrete fracture network model into a format so as to import it into COMSOL simulation software; performing geometric repair and multi-physical field coupling on the discrete fracture network model by using the COMSOL simulation software to obtain a virtual discrete fracture model, and performing entity verification, inversion analysis and grouting optimization to obtain an optimal discrete fracture model.

Description

Technical Field

[0001] The present invention relates to the technical field of rock mechanics, and in particular to a COMSOL-based method for constructing a rock mass penetration grouting model containing a complex fracture network. Background Art

[0002] In grouting projects, the type of injected medium significantly influences the slurry flow pattern and effectiveness. Unlike porous, weak formations, fractured rock formations, due to their unique geological characteristics, include uneven fracture distribution, variability in direction, random connectivity, and varying apertures, which significantly influence the slurry diffusion mechanism. As the primary pathway for slurry flow, accurate characterization of the fracture network is crucial for effectively evaluating the dynamic behavior of the grouting process, predicting the slurry diffusion range, and optimizing grouting parameters.

[0003] The analysis of grouting behavior in fractured rock masses is usually based on seepage models. The core of this process lies in a deep understanding of how slurry flows and distributes within fractured rock masses, thereby providing theoretical support for the optimization of grouting technology. The existing theoretical frameworks for grouting in fractured rock masses include: porous media theory, pseudo-continuum media theory, fracture media theory, and pore-fracture dual media theory. In natural fractured rock masses, the permeability of fractures is much stronger than that of the rock matrix. The porous media theory simulates the fracture system as a porous medium. Although this method simplifies the complex structural characteristics of the fracture network, it is unable to accurately capture the permeation path of the slurry in the rock mass, resulting in limited engineering applicability. Based on the equivalent continuum assumption, the quasi-continuum theory proposes that the fracture system and the matrix system coexist continuously at all points in space, and that the macroscopic representation of the seepage field is achieved through the equivalent permeability tensor. This theory is well applicable to rock masses with a high fracture density and small scale. However, for rock masses containing macroscopic through-fractures or dominant seepage channels, the equivalent continuum assumption will lead to significant deviations in the anisotropic characteristics of the permeability tensor, thereby affecting the reliability of the simulation results.

[0004] In contrast, the fractured medium theory is based on the physical fact that slurry migrates mainly in the fracture network, assuming that the permeability of the rock matrix is negligible. Considering that engineering grouting materials usually have high viscosity and small particle size characteristics, they are almost impermeable in the pores of the rock matrix. This theoretical assumption is physically reasonable in most engineering scenarios. The pore-fracture dual medium theory attempts to simultaneously describe the rapid migration of slurry in the fracture network and the slow diffusion process in the matrix pores, and introduces a mass exchange term to characterize the interaction between the two phases of the medium. Although this model can theoretically reflect the strong water conductivity of the fracture system and the high water storage characteristics of the matrix system, the mass transfer mechanism at the fracture-pore interface has not yet been clarified, and existing experimental techniques make it difficult to accurately determine the dynamic exchange parameters of the two phases of the medium. As a result, this theory is still in the stage of mathematical model construction and has not yet formed a mature engineering application system.

[0005] In terms of innovation in numerical modeling methods, existing technologies have attempted to improve model accuracy by combining digital image processing. Patent No. CN202010287865.6 discloses a method for modeling and seepage testing of broken rock masses based on digital image processing, which relates to the field of rock mechanics technology. The steps include: sampling the broken rock mass and taking images in layers; analyzing the image to obtain a histogram of the RGB color components; analyzing the histogram to determine the main features; when the grayscale value difference is small, using pure colors to mark the same features, increase the feature contrast, separate the blocks and voids, and use threshold separation to determine the main features; redrawing the RGB component histogram through Matlab software and embedding it in Comsol software to establish a numerical analysis model, calling the physical and mechanical properties corresponding to the RGB characteristics, and establishing a three-dimensional numerical analysis model; using the maximum and minimum particle sizes to determine the layered acquisition range of the image, and determining the data interpolation range to obtain the corresponding actual physical property range; and establishing a numerical analysis model for the broken rock mass. This method accurately reproduces the spatial structure of the broken rock mass, solving the problem that the rock accumulation is highly discrete and cannot be tested repeatedly.

[0006] Similarly, Patent No. CN202411209289.8 discloses a method for importing digital core CT slices into COMSOL. The method obtains a binary image of the digital core slice image obtained by CT scanning by performing steps such as slice image contrast enhancement and threshold division through MATLAB software; the binary image is then cut to obtain its representative region (REV) of the core, and the boundary coordinates of the particles and the isolated pores contained in the particles are obtained by processing the functions in MATLAB; the COMSOL WITH MATLAB interface is used to draw the image and generate a geometric model for numerical simulation; the above invention can be applied to porous media fields such as seepage, oil extraction, and geological reservoir analysis, and can accurately and efficiently import CT images into COMSOL for various simulations, and has wide application value.

[0007] However, the above-mentioned modeling method based on digital image processing still has significant limitations. Although patents CN202010287865.6 and CN202411209289.8 both use the Monte Carlo method for random sampling in fracture reconstruction, since this method follows the law of large numbers and requires a large number of samples to support it, it has inherent defects such as low sampling efficiency and large sample capacity requirements in practical applications, which can easily lead to statistical collapse of the obtained fracture samples. The fracture network model finally established has significant deviations from the actual engineering geological structure. More importantly, the model constructed by the existing method fails to effectively couple the rheological properties of the grouting material with the seepage law of the fracture network, resulting in the inability to accurately characterize the diffusion dynamics of the slurry in the discrete fracture network. This is specifically manifested in problems such as inaccurate prediction of the slurry diffusion path and distortion of the grouting pressure field distribution, which seriously restricts the guiding value of numerical simulation for engineering practice.

[0008] A comprehensive analysis of the existing technology system shows that the current modeling of fractured rock mass grouting faces three major technical bottlenecks: (1) The traditional random sampling method is inefficient and it is difficult to ensure the integrity of the statistical characteristics of the fracture network under limited sample capacity; (2) The dynamic coupling mechanism between the discrete fracture network model and the grouting process has not yet been established, and there is a lack of mathematical models that reflect the interaction between the time-varying characteristics of slurry viscosity and the spatial heterogeneity of fracture permeability; (3) The existing modeling process lacks a reliable monitoring and correction mechanism, and it is impossible to achieve adaptive optimization of the model through real-time data assimilation. These problems make it difficult for the existing grouting model to meet the requirements of complex fractured rock mass grouting projects for diffusion range prediction accuracy and construction parameter optimization reliability, which seriously restricts the actual application effect of numerical simulation technology in grouting projects. Summary of the Invention

[0009] The present invention aims to provide a COMSOL-based method for constructing a rock mass infiltration grouting model containing a complex fracture network. The method can adopt a Latin hypercube sampling method to follow the rule of selecting only one sample point in each equal probability interval, significantly reduce the variance of the sample points through a stratified method, thereby greatly improving the sampling efficiency and solving the problems of large sample capacity and low sampling efficiency in Monte Carlo simulation. The method also increases the reliability of the established discrete fracture network model through a monitoring correction model, thereby improving the accuracy and efficiency of the simulated slurry grouting diffusion process in the fracture network.

[0010] The present invention utilizes the following technical solutions:

[0011] A COMSOL-based method for constructing a rock mass permeation grouting model containing a complex fracture network comprises the following steps:

[0012] S1: Conduct geological exploration of the rock mass and extract and process the fracture geometry parameters of the rock mass; fracture geometry parameters include fracture length, fracture tendency, fracture inclination, fracture aperture, fracture filling and fracture profile;

[0013] S2: Fit and correct each crack geometric parameter to obtain the optimal distribution probability function;

[0014] S3: Generate a discrete fracture network model using the Latin hypercube sampling algorithm based on the optimal distribution probability function;

[0015] S4: Convert the discrete fracture network model into a new format for importing into COMSOL simulation software;

[0016] S5: Use COMSOL simulation software to perform geometric repair and multi-physics field coupling on the discrete fracture network model to obtain a virtual discrete fracture model;

[0017] S6: Perform physical verification, inversion analysis, and grouting optimization on the virtual discrete fracture model to obtain the optimal discrete fracture model.

[0018] Preferably, step S1 includes the following steps:

[0019] S11: Non-contact geological exploration of rock mass;

[0020] S12: Use drones to photograph cracks on the exposed surface of the rock mass to obtain a sketch of the cracks on the exposed surface of the rock mass; conduct borehole photography of the rock mass to obtain three-dimensional spatial distribution data of cracks in the rock mass borehole; and use high-resolution cameras or laser scanners to obtain point cloud data on the rock mass surface.

[0021] S13: Digitize the crack sketch of the exposed surface of the rock mass to extract the crack length, crack tendency and crack inclination; identify the three-dimensional spatial distribution data of the crack to obtain the crack opening and crack filling; use the edge detection algorithm to extract the contour of the rock mass surface point cloud data to obtain the crack contour;

[0022] S14: Use the lognormal distribution function or the uniform distribution function to classify and analyze the crack length to obtain the crack length distribution information; use the power law distribution function or the uniform distribution function to classify and analyze the crack aperture to obtain the crack aperture distribution information; use the Fisher distribution function or the Bingham distribution function to perform correlation analysis on the crack inclination and crack tendency to obtain the crack angular distribution information.

[0023] Preferably, step S2 includes the following steps:

[0024] S21: Use the maximum likelihood estimation function to perform data fitting on each fracture geometric parameter to obtain several preliminary probability distribution functions;

[0025] S22: Based on the crack length distribution information, the crack aperture distribution information, and the crack angular distribution information, each preliminary probability distribution function is optimized and corrected using a data association algorithm to obtain several secondary probability distribution functions;

[0026] S23: using a fusion algorithm to combine the exponential distribution function or the Weibull distribution function, fusing the quadratic probability distribution functions to obtain a plurality of distribution probability functions;

[0027] S24: Compare and evaluate several distribution probability functions according to the Akaike information criterion to obtain the optimal distribution probability function.

[0028] Preferably, step S3 includes the following steps:

[0029] S301: Based on the correlation between the fracture geometric parameters, a correlation analysis is performed on the optimal distribution probability function using a covariance matrix or a copula function to obtain a joint probability function;

[0030] S302: Perform statistical analysis on the geometric parameters of each type of crack to determine the parameter range, sample number N and sampling interval ;

[0031] S303: Determine the extraction probability of each fracture geometric parameter based on the optimal distribution probability function combined with the joint probability function;

[0032] S304: Divide the parameter range of each fracture geometric parameter into N sub-intervals according to the extraction probability;

[0033] S305: Randomly extract a representative value in each sub-interval , according to the crack parameter sequence code Use the inverse probability distribution function: , calculate the corresponding crack parameter sampling value ;

[0034] S306: Randomly arrange and combine all the crack parameter sampling values to eliminate the unintentional correlation between the parameters and generate a parameter dimension k. The crack sample matrix;

[0035] S307: Classify the crack morphology into simple cracks and complex cracks according to the complexity of the crack morphology;

[0036] S308: For simple fractures: using the fracture generation algorithm to randomly generate simple fracture surfaces according to the fracture sample matrix; for complex fractures: generating complex fracture surfaces according to vertex coordinates combined with fracture geometric parameters;

[0037] S309: embedding both complex fracture surfaces and simple fracture surfaces into the rock matrix, processing the interfaces between the fracture surfaces and the rock matrix, and between the fracture surfaces, and thus completing the construction of the discrete fracture network model;

[0038] S310: Validate, analyze, and optimize discrete fracture network models: Compare the mean, variance, and quantile of sampled fracture parameter values with those of fracture geometry parameters. Evaluate and optimize the grouting diffusion range or rock permeability by adjusting grouting pressure and fracture density.

[0039] Preferably, step S5 includes the following steps:

[0040] S51: Through virtual operations in COMSOL simulation software, tiny gaps in the discrete fracture network model are automatically stitched together. At the same time, cross fractures in the discrete fracture network model are connected using Boolean union, thereby completing the geometric repair of the discrete fracture network model.

[0041] S52: Use the Navier-Stokes two-phase flow module of COMSOL simulation software to simulate the slurry penetration and diffusion process in the discrete fracture network. At the same time, determine the fluid-solid coupling process according to the fracture aperture, and then complete the multi-physics field setting of the discrete fracture network model.

[0042] S53: Using a segregated solver combined with algebraic multigrid in a multiphysics environment to solve the time-varying pressure boundary of a grouting hole. Make settings, represents the original pressure, Indicates the increasing pressure, represents time; and the pore pressure at the far field boundary Fixation, represents the slurry density, represents the acceleration due to gravity, Indicates the slurry depth;

[0043] S54: In the multi-physics field, the discrete fracture network model that has completed geometric repair is fused according to the viscosity time-varying model, grouting holes and far-field boundaries to obtain a virtual discrete fracture model.

[0044] Preferably, step S6 includes the following steps:

[0045] S61: Create a resin-based fracture model based on the virtual discrete fracture model, inject fluorescent slurry, and calculate grouting pressure-flow data;

[0046] S62: Use COMSOL simulation software to perform finite element analysis on the virtual discrete fracture model and record the simulated diffusion radius of the simulated slurry;

[0047] S63: The position of the fluorescent slurry front was recorded by a high-speed camera and compared with the simulated diffusion radius of the simulated slurry;

[0048] S64: Based on the grouting pressure-flow data, the Levenberg-Marquardt algorithm is used to calibrate and perform sensitivity analysis on the virtual discrete fracture model, and the Morris method is combined to select the dominant grouting parameters; the dominant grouting parameters include fracture aperture and slurry viscosity;

[0049] S65: Dynamically modify the rock mass penetration grouting scheme using the monitoring correction model based on the dominant grouting parameters to obtain the optimal rock mass penetration grouting scheme;

[0050] S66: Based on the optimal rock mass penetration grouting scheme and the virtual discrete fracture model, the optimal discrete fracture model is obtained.

[0051] Preferably, the monitoring and correction model first uses the feature extraction layer to extract feature information of the crack geometric parameters, the optimal distribution probability function, the crack parameter sampling values and the grouting dominant parameters to obtain the data feature matrix and the probability feature matrix; then the data feature matrix is iteratively trained several times using the iterative training layer to generate a data feature weight matrix; then the data feature weight matrix is updated and optimized according to the probability feature matrix and the global loss function through the identification output layer to obtain the optimal weight matrix; finally, the data feature matrix is decomposed and reconstructed using the reconstruction decoder through the judgment and correction layer, and the optimal distribution probability function is corrected according to the real-time monitoring results of the optimal weight matrix.

[0052] Preferably, the iterative training layer includes 2 Convolutional layer, 4 Upsampling layer, 4 Downsampling layer, 2 InceptionV4 blocks, 4 InceptionV3 blocks, 4 residual blocks, 4 average pooling layers, 2 inverted residual blocks, 2 Mish activation functions and 2 RELU activation functions; for the data feature matrix, the iterative training layer learning framework is: 2 The convolutional layers are connected in parallel to form two training branches. After the two training branches are spliced using the channel shuffle layer, two InceptionV3 blocks, two residual blocks, one inverted residual block and one RELU activation function are connected in series. The first training branch is connected in series with two Upsampling layer, 2 Downsampling layer, 2 average pooling layers and 1 Mish activation function; the second training branch is connected in series with 2 layers Upsampling layer, 2 Downsampling layer, 2 InceptionV4 blocks, 2 average pooling layers and 1 Mish activation function.

[0053] Preferably, the recognition output layer includes 4 average pooling layers, 3 adaptive pooling layers, 2 fully connected layers, 2 random dropout layers and 1 SoftMax activation function; each adaptive pooling layer is located between 2 average pooling layers, and 2 random dropout layers and 1 SoftMax activation function are located between 2 fully connected layers; the global loss function is: ,

[0054] in, represents the feature extraction layer, represents the iterative training layer, represents the recognition output layer, represents the global loss function, represents the number of iterative training layers, represents the quantization function, represents the calculation accuracy coefficient, Represents the loss function of each layer iterative training, Indicates the preset optimization weight coefficient, Indicates the optimization efficiency of each layer's iterative training.

[0055] Preferably, the judgment correction layer first uses a reconstruction decoder in combination with the continuity equation and the percolation motion equation to perform a dual reconstruction of the structure and attributes of the data feature matrix and the probability feature matrix to generate a reconstructed decoding value; ,

[0056] in, Represents the reconstructed decoded value, represents the PRelu activation function, represents the data feature matrix, represents the dynamic learning parameter matrix, represents the weight parameter matrix, represents the SigMoid activation function, represents transpose, represents the probability feature matrix;

[0057] Then, the reconstructed decoded values are arranged in time sequence according to the real-time monitoring results of the optimal weight matrix to form a grouting time sequence;

[0058] Finally, the judgment correction layer corrects the optimal distribution probability function according to the grouting sequence.

[0059] The present invention adopts the Latin hypercube sampling method to follow the rule of selecting only one sample point in each equal probability interval, and significantly reduces the variance of the sample points through a stratified method, thereby greatly improving the sampling efficiency and solving the problems of large sample capacity and low sampling efficiency in Monte Carlo simulation; and increases the reliability of the established discrete fracture network model through a monitoring correction model, thereby improving the accuracy and efficiency of the simulated slurry grouting diffusion process in the fracture network. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or related technologies, the following briefly introduces the drawings required for use in the embodiments or related technologies. Obviously, the drawings described below are only embodiments of the present invention. Those skilled in the art can also derive other drawings based on the provided drawings without inventive effort.

[0061] Figure 1 A block diagram of the method for constructing a penetration grouting model;

[0062] Figure 2 is a diagram of a discrete fracture network model;

[0063] Figure 3 The discrete fracture network model diagram of the coupled tunnel and grouting hole;

[0064] Figure 4 This is a cloud diagram showing the change of slurry diffusion range over time. DETAILED DESCRIPTION

[0065] The present invention is described in detail below with reference to the accompanying drawings and embodiments:

[0066] like Figure 1-Figure 4 As shown, the COMSOL-based method for constructing a rock mass permeation grouting model containing a complex fracture network according to the present invention comprises the following steps:

[0067] S1: Conduct geological exploration of the rock mass and extract and process the fracture geometry parameters of the rock mass; fracture geometry parameters include fracture length, fracture tendency, fracture inclination, fracture aperture, fracture filling and fracture profile;

[0068] S2: Fit and correct each crack geometric parameter to obtain the optimal distribution probability function;

[0069] S3: Generate a discrete fracture network model using the Latin hypercube sampling algorithm based on the optimal distribution probability function;

[0070] S4: Convert the discrete fracture network model into a new format for importing into COMSOL simulation software;

[0071] S5: Use COMSOL simulation software to perform geometric repair and multi-physics field coupling on the discrete fracture network model to obtain a virtual discrete fracture model;

[0072] S6: Perform physical verification, inversion analysis, and grouting optimization on the virtual discrete fracture model to obtain the optimal discrete fracture model.

[0073] In the present invention, step S1 includes the following steps:

[0074] S11: Non-contact geological exploration of rock mass;

[0075] S12: Use drones to photograph cracks on the exposed surface of the rock mass to obtain a sketch of the cracks on the exposed surface of the rock mass; conduct borehole photography of the rock mass to obtain three-dimensional spatial distribution data of cracks in the rock mass borehole; and use high-resolution cameras or laser scanners to obtain point cloud data on the rock mass surface.

[0076] S13: Digitize the crack sketch of the exposed surface of the rock mass to extract the crack length, crack tendency and crack inclination; identify the three-dimensional spatial distribution data of the crack to obtain the crack opening and crack filling; use the edge detection algorithm to extract the contour of the rock mass surface point cloud data to obtain the crack contour;

[0077] S14: Use the lognormal distribution function or the uniform distribution function to classify and analyze the crack length to obtain the crack length distribution information; use the power law distribution function or the uniform distribution function to classify and analyze the crack aperture to obtain the crack aperture distribution information; use the Fisher distribution function or the Bingham distribution function to perform correlation analysis on the crack inclination and crack tendency to obtain the crack angular distribution information.

[0078] In this embodiment, the crack geometric parameters are shown in Table 1:

[0079]

[0080] In the present invention, step S2 includes the following steps:

[0081] S21: Use the maximum likelihood estimation function to perform data fitting on each fracture geometric parameter to obtain several preliminary probability distribution functions;

[0082] S22: Based on the crack length distribution information, the crack aperture distribution information, and the crack angular distribution information, each preliminary probability distribution function is optimized and corrected using a data association algorithm to obtain several secondary probability distribution functions;

[0083] S23: using a fusion algorithm to combine the exponential distribution function or the Weibull distribution function, fusing the quadratic probability distribution functions to obtain a plurality of distribution probability functions;

[0084] S24: Compare and evaluate several distribution probability functions according to the Akaike information criterion to obtain the optimal distribution probability function.

[0085] In the present invention, step S3 includes the following steps:

[0086] S301: Based on the correlation between the fracture geometric parameters, a correlation analysis is performed on the optimal distribution probability function using a covariance matrix or a copula function to obtain a joint probability function;

[0087] S302: Perform statistical analysis on the geometric parameters of each type of crack to determine the parameter range, sample number N and sampling interval ;

[0088] S303: Determine the extraction probability of each fracture geometric parameter based on the optimal distribution probability function combined with the joint probability function;

[0089] S304: Divide the parameter range of each fracture geometric parameter into N sub-intervals according to the extraction probability;

[0090] S305: Randomly extract a representative value in each sub-interval , according to the crack parameter sequence code Use the inverse probability distribution function: , calculate the corresponding crack parameter sampling value ;

[0091] S306: Randomly arrange and combine all the crack parameter sampling values to eliminate the unintentional correlation between the parameters and generate a parameter dimension k. The crack sample matrix;

[0092] S307: Classify the crack morphology into simple cracks and complex cracks according to the complexity of the crack morphology.

[0093] In this embodiment, simple cracks generally refer to cracks with a relatively simple shape and structure, and mainly include the following types:

[0094] 1. Linear cracks: The cracks are straight, extending in a relatively consistent direction and with a simple shape;

[0095] 2. Planar fissure: The fissure surface is relatively flat, without obvious curvature or bifurcation;

[0096] 3. Single fissure: The fissure has no obvious branches or intersections and exists independently.

[0097] Complex fissures refer to fissures with diverse shapes and complex structures, mainly including the following types:

[0098] 1. Curved cracks: The cracks are curved, with variable extension directions and irregular shapes;

[0099] 2. Bifurcation type crack: The crack bifurcates during the extension process, forming multiple branches;

[0100] 3. Network cracks: Multiple cracks intersect each other to form a network structure;

[0101] 4. Irregular cracks: The cracks have complex morphology, no obvious pattern, and may contain multiple morphological features;

[0102] S308: For simple fractures: using the fracture generation algorithm to randomly generate simple fracture surfaces according to the fracture sample matrix; for complex fractures: generating complex fracture surfaces according to vertex coordinates combined with fracture geometric parameters;

[0103] S309: embedding both complex fracture surfaces and simple fracture surfaces into the rock matrix, processing the interfaces between the fracture surfaces and the rock matrix, and between the fracture surfaces, thereby completing the construction of the discrete fracture network model;

[0104] S310: Validate, analyze, and optimize discrete fracture network models: Compare the mean and variance of fracture parameter sampling values with the fracture geometry parameters, and evaluate and optimize the grouting diffusion range or rock permeability by adjusting the grouting pressure and fracture density.

[0105] In the present invention, step S5 includes the following steps:

[0106] S51: Through virtual operations in COMSOL simulation software, small gaps (<0.1 mm) in the discrete fracture network model are automatically stitched together. At the same time, cross fractures in the discrete fracture network model are connected using Boolean union, thereby completing the geometric repair of the discrete fracture network model.

[0107] S52: Use the Navier-Stokes two-phase flow module of COMSOL simulation software to simulate the slurry penetration and diffusion process in the discrete fracture network. At the same time, determine the fluid-solid coupling process according to the fracture aperture, and then complete the multi-physics field setting of the discrete fracture network model.

[0108] S53: Using a segregated solver combined with algebraic multigrid in a multiphysics environment to solve the time-varying pressure boundary of a grouting hole. Make settings, represents the original pressure, Indicates the increasing pressure, represents time; and the pore pressure at the far field boundary Fixation, represents the slurry density, represents the acceleration due to gravity, Indicates the slurry depth;

[0109] S54: In the multi-physics field, the discrete fracture network model that has completed geometric repair is fused according to the viscosity time-varying model, grouting holes and far-field boundaries to obtain a virtual discrete fracture model.

[0110] In the present invention, step S6 includes the following steps:

[0111] S61: Create a resin-based fracture model based on the virtual discrete fracture model, inject fluorescent slurry, and calculate grouting pressure-flow data;

[0112] S62: Use COMSOL simulation software to perform finite element analysis on the virtual discrete fracture model and record the simulated diffusion radius of the simulated slurry;

[0113] S63: The position of the fluorescent slurry front was recorded by a high-speed camera and compared with the simulated diffusion radius of the simulated slurry;

[0114] S64: Based on the grouting pressure-flow data, the Levenberg-Marquardt algorithm is used to calibrate and perform sensitivity analysis on the virtual discrete fracture model, and the Morris method is combined to select the dominant grouting parameters; the dominant grouting parameters include fracture aperture and slurry viscosity;

[0115] S65: Dynamically modify the rock mass penetration grouting scheme using the monitoring correction model based on the dominant grouting parameters to obtain the optimal rock mass penetration grouting scheme;

[0116] S66: Based on the optimal rock mass penetration grouting scheme and the virtual discrete fracture model, the optimal discrete fracture model is obtained.

[0117] In the present invention, the monitoring correction model first uses the feature extraction layer to extract feature information of crack geometric parameters, optimal distribution probability function, crack parameter sampling values and grouting dominant parameters to obtain data feature matrix and probability feature matrix.

[0118] In this embodiment, the feature extraction layer includes 4 Convolutional layer, 3 Max pooling layer and 1 normalization layer; every 2 The convolutional layer is connected in parallel with 1 The maximum pooling layers are connected in series and then connected in parallel one at a time. A large pooling layer and a normalization layer are used to fuse the characteristic information of the fracture geometry parameters, fracture parameter sampling values and grouting dominant parameters into a data feature matrix; at the same time, the characteristic information of the optimal distribution probability function is fused into a probability feature matrix;

[0119] Then, the data feature matrix is iteratively trained several times using the iterative training layer to generate a data feature weight matrix;

[0120] In the present invention, the iterative training layer includes 2 Convolutional layer, 4 Upsampling layer, 4 Downsampling layer, 2 InceptionV4 blocks, 4 InceptionV3 blocks, 4 residual blocks, 4 average pooling layers, 2 inverted residual blocks, 2 Mish activation functions and 2 RELU activation functions; for the data feature matrix, the iterative training layer learning framework is: 2 The convolutional layers are connected in parallel to form two training branches. After the two training branches are spliced using the channel shuffle layer, two InceptionV3 blocks, two residual blocks, one inverted residual block and one RELU activation function are connected in series. The first training branch is connected in series with two Upsampling layer, 2 Downsampling layer, 2 average pooling layers and 1 Mish activation function; the second training branch is connected in series with 2 layers Upsampling layer, 2 Downsampling layer, 2 InceptionV4 blocks, 2 average pooling layers and 1 Mish activation function;

[0121] Then, the data feature weight matrix is updated and optimized according to the probability feature matrix and the global loss function through the identification output layer to obtain the optimal weight matrix.

[0122] In the present invention, the recognition output layer includes 4 average pooling layers, 3 adaptive pooling layers, 2 fully connected layers, 2 random dropout layers and 1 SoftMax activation function; each adaptive pooling layer is located between the two average pooling layers, and the two random dropout layers and 1 SoftMax activation function are located between the two fully connected layers; the global loss function is: ,

[0123] in, represents the feature extraction layer, represents the iterative training layer, represents the recognition output layer, represents the global loss function, represents the number of iterative training layers, represents the quantization function, represents the calculation accuracy coefficient, Represents the loss function of each layer iterative training, Indicates the preset optimization weight coefficient, Indicates the optimization efficiency of each layer’s iterative training;

[0124] Finally, the data feature matrix is decomposed and reconstructed by the reconstruction decoder through the judgment correction layer, and the optimal distribution probability function is corrected according to the real-time monitoring results of the optimal weight matrix.

[0125] In the present invention, the judgment correction layer first uses a reconstruction decoder combined with the continuity equation and the percolation motion equation to perform dual structural and attribute reconstruction on the data feature matrix and the probability feature matrix to generate a reconstructed decoding value;

[0126] ,

[0127] in, Represents the reconstructed decoded value, represents the PRelu activation function, represents the data feature matrix, represents the dynamic learning parameter matrix, represents the weight parameter matrix, represents the SigMoid activation function, represents transpose, represents the probability feature matrix;

[0128] Then, the reconstructed decoded values are arranged in time sequence according to the real-time monitoring results of the optimal weight matrix to form a grouting time sequence;

[0129] Finally, the judgment correction layer corrects the optimal distribution probability function according to the grouting sequence.

[0130] In this example, the slurry's diffusion process within the fracture network was simulated using the Navier-Stokes two-phase flow module in COMSOL simulation software. The slurry's diffusion range is a key concern in grouting projects. Grouting is considered a process whereby the slurry phase displaces the water or air phase originally present in the fractures. The slurry front, i.e., the interface between the slurry and the water or air phase, is tracked using a phase function.

[0131] By setting the boundary conditions of each phase and the material properties of the fluid at the initial moment, and through meshing and simulation, the diffusion process of the slurry in the fracture network can be solved.

[0132] Embodiment: Conduct non-contact geological exploration on rock mass; photograph cracks on rock surface exposure by drone to obtain sketch of cracks on rock surface exposure; conduct drilling photography on rock mass to obtain three-dimensional spatial distribution data of cracks in rock borehole; obtain point cloud data of rock surface by high-resolution camera or laser scanner; digitize sketch of cracks on rock surface exposure to extract crack length, crack tendency and crack inclination; identify three-dimensional spatial distribution data of cracks to obtain crack opening and crack filling; use edge detection algorithm to extract outline of point cloud data on rock surface to obtain crack outline; use uniform distribution function to classify and analyze crack length to obtain crack length distribution information; use uniform distribution function to classify and analyze crack opening to obtain Obtain the distribution information of crack opening; use the Fisher distribution function to perform correlation analysis on the crack inclination and crack tendency to obtain the crack angular distribution information; use the maximum likelihood estimation function to perform data fitting on each crack geometric parameter to obtain several preliminary probability distribution functions; according to the crack length distribution information, crack opening distribution information and crack angular distribution information, use the data association algorithm to optimize and correct each preliminary probability distribution function respectively to obtain several quadratic probability distribution functions; use the construction fusion algorithm combined with the exponential distribution function or the Weibull distribution function to fuse the quadratic probability distribution functions to obtain several distribution probability functions; compare and evaluate several distribution probability functions according to the Akaike information criterion to obtain the optimal distribution probability function.

[0133] According to the correlation between the fracture geometric parameters, the covariance matrix or Copula function is used to perform correlation analysis on the optimal distribution probability function to obtain the joint probability function; statistical analysis is performed on each type of fracture geometric parameter to determine the parameter range, sample number N and sampling interval. ; Determine the extraction probability of each fracture geometric parameter based on the optimal distribution probability function combined with the joint probability function; Divide the parameter range of each fracture geometric parameter into N subintervals according to the extraction probability; Randomly extract a representative value in each subinterval , according to the crack parameter sequence code Use the inverse probability distribution function: , calculate the corresponding crack parameter sampling value ; Randomly arrange and combine all crack parameter sampling values to eliminate the unintentional correlation between parameters and generate a parameter dimension k The crack sample matrix is obtained; according to the complexity of the crack morphology, the crack morphology is divided into simple cracks and complex cracks; for simple cracks: the crack generation algorithm is used to randomly generate simple crack surfaces according to the crack sample matrix; for complex cracks: the complex crack surfaces are generated according to the vertex coordinates combined with the crack geometric parameters.

[0134] Both complex and simple fracture surfaces are embedded in the rock matrix, and the interfaces between the fracture surfaces and the rock matrix, as well as between the fracture surfaces, are processed to complete the construction of a discrete fracture network model. The discrete fracture network model is verified, analyzed, and optimized by comparing the mean, variance, and quantile of the fracture parameter sampling values with the fracture geometric parameters. At the same time, by adjusting the grouting pressure and fracture density, the grouting diffusion range or rock permeability is evaluated and optimized.

[0135] The discrete fracture network model is converted into a format for import into COMSOL simulation software. Through the virtual operation of COMSOL simulation software, the tiny gaps in the discrete fracture network model are automatically stitched, and the cross fractures in the discrete fracture network model are connected using Boolean union, thereby completing the geometric repair of the discrete fracture network model. The Navier-Stokes two-phase flow module of COMSOL simulation software is used to simulate the infiltration and diffusion process of slurry in the discrete fracture network, and the fluid-solid coupling process is determined according to the fracture opening, thereby completing the multi-physics field setting of the discrete fracture network model. In the multi-physics field, a segregated solver combined with algebraic multigrid is used to solve the time-varying pressure boundary of the grouting hole. Make settings, represents the original pressure, Indicates the increasing pressure, represents time; and the pore pressure at the far field boundary Fixation, represents the slurry density, represents the acceleration due to gravity, represents the slurry depth; in the multi-physics field, the discrete fracture network model that has completed geometric repair is fused according to the viscosity time-varying model, grouting holes and far-field boundaries to obtain a virtual discrete fracture model.

[0136] A resin-based fracture model was fabricated based on the virtual discrete fracture model, and fluorescent slurry was injected to calculate the grouting pressure-flow data. The virtual discrete fracture model was subjected to finite element analysis using COMSOL simulation software, and the simulated diffusion radius of the simulated slurry was recorded. The position of the fluorescent slurry front was recorded with a high-speed camera and compared with the simulated diffusion radius of the simulated slurry. The virtual discrete fracture model was calibrated and sensitivity analyzed based on the grouting pressure-flow data using the Levenberg-Marquardt algorithm, and the dominant grouting parameters were screened using the Morris method. The dominant grouting parameters included fracture aperture and slurry viscosity. The monitoring correction model was used to dynamically correct the rock infiltration grouting scheme based on the dominant grouting parameters to obtain the optimal rock infiltration grouting scheme. The optimal discrete fracture model was obtained by combining the optimal rock infiltration grouting scheme with the virtual discrete fracture model.

Claims

1. A COMSOL-based method for constructing a rock mass permeation grouting model containing a complex fracture network, characterized by: The following steps are involved: S1: Conduct geological exploration of the rock mass and extract and process the fracture geometry parameters of the rock mass; fracture geometry parameters include fracture length, fracture tendency, fracture inclination, fracture aperture, fracture filling and fracture profile; S2: Fit and correct each crack geometric parameter to obtain the optimal distribution probability function; S3: Generate a discrete fracture network model using the Latin hypercube sampling algorithm based on the optimal distribution probability function; S4: Convert the discrete fracture network model into a new format for importing into COMSOL simulation software; S5: Use COMSOL simulation software to perform geometric repair and multi-physics field coupling on the discrete fracture network model to obtain a virtual discrete fracture model; S6: performing physical verification, inversion analysis, and grouting optimization on the virtual discrete fracture model to obtain an optimal discrete fracture model; said step S6 includes the following steps: making a resin-based fracture model based on the virtual discrete fracture model, injecting fluorescent slurry, and measuring grouting pressure-flow data; using the Levenberg-Marquardt algorithm based on the grouting pressure-flow data, calibrating and performing sensitivity analysis on the virtual discrete fracture model, and combining the Morris method to screen the dominant grouting parameters; the dominant grouting parameters include fracture aperture and slurry viscosity; using the monitoring correction model based on the dominant grouting parameters, dynamically correcting the rock penetration grouting scheme to obtain the optimal rock penetration grouting scheme; and combining the optimal rock penetration grouting scheme with the virtual discrete fracture model to obtain the optimal discrete fracture model; The monitoring correction model first uses the feature extraction layer to extract feature information of the fracture geometric parameters, the optimal distribution probability function, the fracture parameter sampling values, and the grouting dominant parameters to obtain the data feature matrix and the probability feature matrix. The iterative training layer then performs several iterative training on the data feature matrix to generate the data feature weight matrix. The identification output layer then updates and optimizes the data feature weight matrix based on the probability feature matrix and the global loss function to obtain the optimal weight matrix. Finally, the judgment correction layer uses the reconstruction decoder combined with the continuity equation and the seepage motion equation to perform dual structural and attribute reconstruction of the data feature matrix and the probability feature matrix to generate the reconstructed decoding value: Where S represents the reconstructed decoding value, P() represents the PRelu activation function, H represents the data feature matrix, A represents the dynamic learning parameter matrix, W represents the weight parameter matrix, SM() represents the SigMoid activation function, T represents the transpose, and M represents the probability feature matrix; The reconstructed decoded values are arranged in time sequence according to the real-time monitoring results of the optimal weight matrix to form a grouting time sequence; The judgment correction layer corrects the optimal distribution probability function according to the grouting sequence.

2. The method for constructing a rock mass infiltration grouting model containing a complex fracture network based on COMSOL according to claim 1, characterized in that: The step S1 includes the following steps: S11: Non-contact geological exploration of rock mass; S12: Use drones to photograph cracks on the exposed surface of the rock mass to obtain a sketch of the cracks on the exposed surface of the rock mass; conduct borehole photography of the rock mass to obtain three-dimensional spatial distribution data of cracks in the rock mass borehole; and use high-resolution cameras or laser scanners to obtain point cloud data on the rock mass surface. S13: Digitize the crack sketch of the exposed surface of the rock mass to extract the crack length, crack tendency and crack inclination; identify the three-dimensional spatial distribution data of the crack to obtain the crack opening and crack filling; use the edge detection algorithm to extract the contour of the rock mass surface point cloud data to obtain the crack contour; S14: Use the lognormal distribution function or the uniform distribution function to classify and analyze the crack length to obtain the crack length distribution information; use the power law distribution function or the uniform distribution function to classify and analyze the crack aperture to obtain the crack aperture distribution information; use the Fisher distribution function or the Bingham distribution function to perform correlation analysis on the crack inclination and crack tendency to obtain the crack angular distribution information.

3. The method for constructing a rock mass infiltration grouting model containing a complex fracture network based on COMSOL according to claim 1, characterized in that: The step S2 includes the following steps: S21: Use the maximum likelihood estimation function to perform data fitting on each fracture geometric parameter to obtain several preliminary probability distribution functions; S22: Based on the crack length distribution information, the crack aperture distribution information, and the crack angular distribution information, each preliminary probability distribution function is optimized and corrected using a data association algorithm to obtain several secondary probability distribution functions; S23: using a fusion algorithm to combine the exponential distribution function or the Weibull distribution function, fusing the quadratic probability distribution functions to obtain a plurality of distribution probability functions; S24: Compare and evaluate several distribution probability functions according to the Akaike information criterion to obtain the optimal distribution probability function.

4. The method for constructing a rock mass infiltration grouting model containing a complex fracture network based on COMSOL according to claim 1, characterized in that: The step S3 includes the following steps: S301: Based on the correlation between the fracture geometric parameters, a correlation analysis is performed on the optimal distribution probability function using a covariance matrix or a copula function to obtain a joint probability function; S302: Perform statistical analysis on each type of crack geometric parameters to determine the parameter range, sample number N and sampling interval [y min ,y max ]; S303: Determine the extraction probability of each fracture geometric parameter based on the optimal distribution probability function combined with the joint probability function; S304: Divide the parameter range of each fracture geometric parameter into N sub-intervals according to the extraction probability; S305: Randomly extract a representative value Y in each sub-interval i , using the inverse probability distribution function according to the crack parameter sequence code i: , calculate the corresponding crack parameter sampling value x i ; S306: Randomly arrange and combine all crack parameter sampling values to eliminate unintentional correlations between parameters, and generate an N×k crack sample matrix with parameter dimension k; S307: Classify the crack morphology into simple cracks and complex cracks according to the complexity of the crack morphology; S308: For simple fractures: using the fracture generation algorithm to randomly generate simple fracture surfaces according to the fracture sample matrix; for complex fractures: generating complex fracture surfaces according to vertex coordinates combined with fracture geometric parameters; S309: embedding both complex fracture surfaces and simple fracture surfaces into the rock matrix, processing the interfaces between the fracture surfaces and the rock matrix, and between the fracture surfaces, thereby completing the construction of the discrete fracture network model; S310: Validate, analyze, and optimize discrete fracture network models: Compare the mean, variance, and quantile of sampled fracture parameter values with those of fracture geometry parameters. Evaluate and optimize the grouting diffusion range or rock permeability by adjusting grouting pressure and fracture density.

5. The method for constructing a rock mass infiltration grouting model containing a complex fracture network based on COMSOL according to claim 1, characterized in that: The step S5 includes the following steps: S51: Through virtual operations in COMSOL simulation software, tiny gaps in the discrete fracture network model are automatically stitched together. At the same time, cross fractures in the discrete fracture network model are connected using Boolean union, thereby completing the geometric repair of the discrete fracture network model. S52: Use the Navier-Stokes two-phase flow module of COMSOL simulation software to simulate the slurry penetration and diffusion process in the discrete fracture network. At the same time, determine the fluid-solid coupling process according to the fracture aperture, and then complete the multi-physics field setting of the discrete fracture network model. S53: Using a segregated solver combined with algebraic multigrid in a multiphysics environment to solve the time-varying pressure boundary of a grouting hole. Set, p0 represents the original pressure, β represents the variable pressure, t represents the time; and the pore pressure at the far field boundary is fixed, ρ represents the slurry density, g represents the acceleration due to gravity, and h represents the slurry depth; S54: In the multi-physics field, the discrete fracture network model that has completed geometric repair is fused according to the viscosity time-varying model, grouting holes and far-field boundaries to obtain a virtual discrete fracture model.

6. The method for constructing a rock mass infiltration grouting model containing a complex fracture network based on COMSOL according to claim 1, characterized in that: The step S6 further includes the following steps: The virtual discrete fracture model was analyzed by finite element analysis using COMSOL simulation software, and the simulated diffusion radius of the simulated slurry was recorded; The position of the fluorescent slurry front was recorded by a high-speed camera and compared with the simulated diffusion radius of the simulated slurry.

7. The method for constructing a rock mass infiltration grouting model containing a complex fracture network based on COMSOL according to claim 1, characterized in that: The iterative training layer includes two 1×1 convolutional layers, four 3×1 upsampling layers, four 1×3 downsampling layers, two InceptionV4 blocks, four InceptionV3 blocks, four residual blocks, four average pooling layers, two inverted residual blocks, two Mish activation functions, and two RELU activation functions. For the data feature matrix, the iterative training layer learning framework is: two 1×1 convolutional layers are connected in parallel to form two training branches. After the two training branches are spliced using the channel shuffle layer, two InceptionV3 blocks, two residual blocks, one inverted residual block and one RELU activation function are connected in series. The first training branch connects two 3×1 upsampling layers, two 1×3 downsampling layers, two average pooling layers and one Mish activation function in series. The second training branch connects two 3×1 upsampling layers, two 1×3 downsampling layers, two InceptionV4 blocks, two average pooling layers and one Mish activation function in series.

Images