Complex fracture network rock mass permeation grouting model construction method based on COMSOL

Abstract

The invention discloses a method for constructing a complex fracture network rock mass permeation grouting model based on COMSOL, relates to the technical field of rock mass mechanics, and solves the problem that an established model is greatly different from actual engineering due to collapse of an obtained fracture sample caused by the defects of low sampling efficiency and large sample capacity. The rock mass permeation grouting model construction method comprises the steps that geological exploration is conducted on a rock mass, and fracture geometric parameters of the rock mass are extracted and processed; fitting correction is conducted on the geometric parameters of each fracture, and an optimal distribution probability function is obtained; according to the optimal distribution probability function, utilizing a Latin hypercube sampling algorithm to generate a discrete fracture network model; performing format conversion on the discrete fracture network model so as to import COMSOL simulation software; performing geometric repair and multi-physical field coupling on the discrete fracture network model through 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 invention relates to the technical field of rock mass mechanics, and in particular to a method for constructing a rock mass penetration grouting model containing a complex fracture network based on COMSOL. Background Art

[0002] In grouting engineering, the type of injected medium has a significant impact on the slurry flow mode and effect. Different from porous medium soft formations, fractured rock formations have unique geological characteristics, including uneven distribution of fractures, variability of directions, randomness of connectivity, and differences in opening, which have an important impact on the diffusion mechanism of slurry. As the main channel for slurry flow, the accurate characterization of the fracture network is of great significance for effectively evaluating the dynamic behavior of the grouting process, predicting the slurry diffusion range, and optimizing the grouting parameters.

[0003] The analysis of grouting behavior in fractured rock mass is usually based on the seepage model. The core of this process is to deeply understand how the slurry flows and distributes inside the fractured rock mass, so as to provide theoretical support for the optimization of grouting technology. The existing theoretical frameworks for grouting in fractured rock mass include: porous media theory, quasi-continuous media theory, fractured media theory and pore-fracture dual media theory. In natural fractured rock mass, 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, which limits its engineering applicability. The quasi-continuum theory is based on the equivalent continuum medium hypothesis. It proposes that the fracture system and the matrix system coexist continuously at all points in space, and realizes the macroscopic characterization of the seepage field through the equivalent permeability tensor. This theory has good applicability in rock masses with high fracture density and small scale. However, for rock masses containing macroscopic through fractures or dominant seepage channels, the equivalent continuum hypothesis 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 fracture 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 the existing experimental technology is difficult to accurately determine the dynamic exchange parameters of the two-phase medium. As a result, the 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, the prior art has attempted to improve model accuracy by combining digital image processing. Patent No. CN202010287865.6 discloses a broken rock modeling and seepage test method 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 the histogram of the RGB color components; analyzing the histogram to determine the main features; when the gray value difference is small, use pure colors to mark the same features, increase the feature contrast, separate the blocks and voids, use threshold separation, and determine the main features; redraw the RGB component histogram through Matlab software, and embed it in Comsol software to establish a numerical analysis model, call the physical and mechanical properties corresponding to the RGB characteristics, and establish a three-dimensional numerical analysis model; use the maximum particle size and the minimum particle size to determine the layered acquisition range of the image, and determine the data interpolation range to obtain the corresponding actual physical property range; establish a numerical analysis model for broken rock mass. This method accurately reproduces the spatial broken rock structure and solves the problem that the rock accumulation is too 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 the 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 images and generate geometric models for numerical simulation; the above invention can be applied in 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 acquired 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, which is specifically manifested in the inaccurate prediction of the slurry diffusion path and the 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 technical 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 purpose of the present invention is to provide a method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL, which can use a Latin hypercube sampling method to follow the rule of selecting only one sample point in each equal probability interval, and 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; and increase 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: A method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL, comprising the following steps; S1: Conduct geological exploration of rock mass, and extract and process fracture geometry parameters of rock mass; fracture geometry parameters include fracture length, fracture tendency, fracture inclination, fracture opening, fracture filling and fracture profile; S2: Fit and correct each fracture geometry 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 to be imported 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: Perform physical verification, inversion analysis and grouting optimization on the virtual discrete fracture model to obtain the optimal discrete fracture model.

[0011] Preferably, 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; perform borehole photography on 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 surface of the rock mass. S13: Digitize the crack sketch of the exposed surface of the rock mass to extract the crack length, crack tendency and crack dip angle; 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 surface point cloud data to obtain the crack contour; S14: Use the log-normal 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.

[0012] Preferably, 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: according to the crack length distribution information, the crack opening distribution information and the crack angular distribution information, each preliminary probability distribution function is optimized and corrected by using a data association algorithm to obtain a number of secondary probability distribution functions; S23: using a fusion algorithm to combine an exponential distribution function or a Weibull distribution function, to fuse the quadratic probability distribution functions, and 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.

[0013] Preferably, step S3 includes the following steps: S301: According to 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 fracture geometric parameters to determine the parameter range, sample number N and sampling interval ; S303: Determine the extraction probability of each fracture geometric parameter according to 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 in each sub-interval , according to the crack parameter sequence code Using the inverse probability distribution function: , calculate the corresponding crack parameter sampling value ; 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 of S307: Classify the fracture morphology into simple fractures and complex fractures according to the complexity of the fracture morphology; S308: For simple fractures: using fracture generation algorithm to randomly generate simple fracture surfaces according to 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 interfaces between fracture surfaces and the rock matrix, and between fracture surfaces, and thus completing the construction of a discrete fracture network model; S310: Verify, analyze and optimize the discrete fracture network model: compare the mean, variance and quantile of the 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.

[0014] Preferably, step S5 includes the following steps: S51: 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; S52: The Navier-Stokes two-phase flow module of COMSOL simulation software is used to simulate the infiltration and diffusion process of slurry in discrete fracture networks. At the same time, the fluid-solid coupling process is determined according to the fracture aperture, thereby completing 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 To set up, represents the original pressure, represents 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; 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.

[0015] Preferably, step S6 includes the following steps: S61: making a resin-based fracture model according to the virtual discrete fracture model, injecting fluorescent slurry, and calculating grouting pressure-flow data; 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; 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; 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 screen the dominant grouting parameters; the dominant grouting parameters include fracture aperture and slurry viscosity; S65: Dynamically correct the rock mass penetration grouting scheme using the monitoring correction model according to the dominant grouting parameters to obtain the optimal rock mass penetration grouting scheme; S66: Based on the optimal rock mass penetration grouting scheme combined with the virtual discrete fracture model, the optimal discrete fracture model is obtained.

[0016] Preferably, the monitoring and correction model first uses the feature extraction layer to extract feature information of the fracture geometry parameters, the optimal distribution probability function, the fracture parameter sampling values ​​and the grouting dominant parameters, respectively, 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 by 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.

[0017] Preferably, the iterative training layer includes 2 Convolutional layers, 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 connects 2 layers in series Upsampling layer, 2 Downsampling layer, 2 InceptionV4 blocks, 2 average pooling layers and 1 Mish activation function.

[0018] 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, 2 random dropout layers and 1 SoftMax activation function are located between 2 fully connected layers; the global loss function is: , 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 iterative training loss function for each layer, Represents the preset optimization weight coefficient, Indicates the optimization efficiency of each layer’s iterative training.

[0019] 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; , 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; 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; Finally, the judgment correction layer corrects the optimal distribution probability function according to the grouting sequence.

[0020] 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

[0021] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or related technologies, the drawings required for use in the embodiments or related technology descriptions are briefly introduced below. Obviously, the drawings in the following description are only embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on the provided drawings without creative work.

[0022] Figure 1 A block diagram of the method for modeling infiltration grouting; Figure 2 is a diagram of a discrete fracture network model; Figure 3The discrete fracture network model diagram of the coupled tunnel and grouting hole; Figure 4 This is a cloud diagram showing the change of slurry diffusion range over time. DETAILED DESCRIPTION

[0023] The present invention is described in detail below with reference to the accompanying drawings and embodiments: like Figure 1-Figure 4 As shown, the method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL of the present invention comprises the following steps: S1: Conduct geological exploration of rock mass, and extract and process fracture geometry parameters of rock mass; fracture geometry parameters include fracture length, fracture tendency, fracture inclination, fracture opening, fracture filling and fracture profile; S2: Fit and correct each fracture geometry 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 to be imported 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: Perform physical verification, inversion analysis and grouting optimization on the virtual discrete fracture model to obtain the optimal discrete fracture model.

[0024] In the present invention, 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; perform borehole photography on 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 surface of the rock mass. S13: Digitize the crack sketch of the exposed surface of the rock mass to extract the crack length, crack tendency and crack dip angle; 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 surface point cloud data to obtain the crack contour; S14: Use the log-normal 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.

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

[0026] In the present invention, 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: according to the crack length distribution information, the crack opening distribution information and the crack angular distribution information, each preliminary probability distribution function is optimized and corrected by using a data association algorithm to obtain a number of secondary probability distribution functions; S23: using a fusion algorithm to combine an exponential distribution function or a Weibull distribution function, to fuse the quadratic probability distribution functions, and 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.

[0027] In the present invention, step S3 includes the following steps: S301: According to 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 fracture geometric parameters to determine the parameter range, sample number N and sampling interval ; S303: Determine the extraction probability of each fracture geometric parameter according to 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 in each sub-interval , according to the crack parameter sequence code Using the inverse probability distribution function: , calculate the corresponding crack parameter sampling value ; 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 of S307: Classify the fracture morphology into simple fractures and complex fractures according to the complexity of the fracture morphology.

[0028] In this embodiment, simple cracks generally refer to cracks with a relatively simple shape and structure, and mainly include the following types: 1. Linear cracks: The cracks are straight, with a relatively consistent extension direction and a simple shape;

[0029] 2. Planar fissure: The fissure surface is relatively flat, without obvious bends or bifurcations;

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

[0031] Complex fissures refer to fissures with diverse shapes and complex structures, mainly including the following types: 1. Curved cracks: The cracks are curved, with variable extension directions and irregular shapes;

[0032] 2. Bifurcated fissure: The fissure bifurcates during the extension process, forming multiple branches;

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

[0034] 4. Irregular cracks: The cracks have complex shapes, no obvious rules, and may contain a variety of morphological features;

[0035] S308: For simple fractures: using fracture generation algorithm to randomly generate simple fracture surfaces according to 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 interfaces between fracture surfaces and the rock matrix, and between fracture surfaces, and thus completing the construction of a discrete fracture network model; S310: Verify, analyze and optimize the discrete fracture network model: compare the mean and variance of the 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.

[0036] In the present invention, step S5 includes the following steps: S51: Through the virtual operation of COMSOL simulation software, the tiny gaps (<0.1mm) 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; S52: The Navier-Stokes two-phase flow module of COMSOL simulation software is used to simulate the infiltration and diffusion process of slurry in discrete fracture networks. At the same time, the fluid-solid coupling process is determined according to the fracture aperture, thereby completing 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 To set up, represents the original pressure, represents 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; 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.

[0037] In the present invention, step S6 includes the following steps: S61: making a resin-based fracture model according to the virtual discrete fracture model, injecting fluorescent slurry, and calculating grouting pressure-flow data; 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; 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; 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 screen the dominant grouting parameters; the dominant grouting parameters include fracture aperture and slurry viscosity; S65: Dynamically correct the rock mass penetration grouting scheme using the monitoring correction model according to the dominant grouting parameters to obtain the optimal rock mass penetration grouting scheme; S66: Based on the optimal rock mass penetration grouting scheme combined with the virtual discrete fracture model, the optimal discrete fracture model is obtained.

[0038] In the present invention, the monitoring correction model firstly 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 respectively to obtain the data feature matrix and the probability feature matrix.

[0039] In this embodiment, the feature extraction layer includes 4 Convolutional layer, 3 Max pooling layer and 1 normalization layer; every 2 The convolutional layers are connected in parallel with 1 The maximum pooling layers are connected in series and then connected in parallel one by one A large pooling layer and a normalization layer are used to fuse the characteristic information of 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; Then, the data feature matrix is ​​iteratively trained several times using the iterative training layer to generate a data feature weight matrix; In the present invention, the iterative training layer includes 2 Convolutional layers, 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 connects 2 layers in series Upsampling layer, 2 Downsampling layer, 2 InceptionV4 blocks, 2 average pooling layers and 1 Mish activation function; Then, the data feature weight matrix is ​​updated and optimized through the identification output layer according to the probability feature matrix and the global loss function to obtain the optimal weight matrix.

[0040] 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 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: , 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 iterative training loss function for each layer, Represents the preset optimization weight coefficient, Indicates the optimization efficiency of each layer’s iterative training; 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.

[0041] In the present invention, 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; , 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; 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; Finally, the judgment correction layer corrects the optimal distribution probability function according to the grouting sequence.

[0042] In this embodiment, the penetration and diffusion process of slurry in the fracture network is simulated using the Navier-Stokes two-phase flow module in the COMSOL simulation software. The slurry diffusion range is the most concerned issue in grouting engineering. Grouting is regarded as a process in which the slurry phase drives to replace the water phase or air phase originally existing in the fracture. The slurry front, that is, the interface between the slurry and the water phase or air phase, is tracked through the phase function.

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

[0044] Embodiment: non-contact geological exploration of rock mass; photographing cracks on the exposed surface of rock mass by drone to obtain a sketch of cracks on the exposed surface of rock mass; performing borehole photography on rock mass to obtain three-dimensional spatial distribution data of cracks in the borehole of rock mass; obtaining point cloud data of rock mass surface by high-resolution camera or laser scanner; digitizing the sketch of cracks on the exposed surface of rock mass to extract crack length, crack tendency and crack inclination; identifying the three-dimensional spatial distribution data of cracks to obtain crack opening and crack filling; using edge detection algorithm to extract contours of point cloud data on the surface of rock mass to obtain crack contours; using uniform distribution function to classify and analyze crack length to obtain crack length distribution information; using 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, the crack opening distribution information and the crack angular distribution information, use the data association algorithm to optimize and correct each preliminary probability distribution function respectively to obtain several secondary probability distribution functions; use the construction fusion algorithm to combine the exponential distribution function or the Weibull distribution function to fuse the secondary 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.

[0045] 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 parameters to determine the parameter range, sample number N and sampling interval ; Determine the extraction probability of each fracture geometry parameter based on the optimal distribution probability function combined with the joint probability function; Divide the parameter range of each fracture geometry parameter into N sub-intervals according to the extraction probability; Randomly extract a representative value in each sub-interval , according to the crack parameter sequence code Using the inverse probability distribution function: , calculate the corresponding crack parameter sampling value ; Randomly arrange and combine all the crack parameter sampling values ​​to eliminate the unintentional correlation between the parameters and generate a parameter dimension of 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: complex crack surfaces are generated according to the vertex coordinates combined with the crack geometric parameters.

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

[0047] The discrete fracture network model is converted into a format for import into the COMSOL simulation software. Through the virtual operation of the 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 the COMSOL simulation software is used to simulate the infiltration and diffusion process of the 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, the segregated solver combined with algebraic multigrid is used to simulate the time-varying pressure boundary of the grouting hole. To set up, represents the original pressure, represents 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; in the multi-physical 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.

[0048] A resin-based fracture model was made according to 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 front position of the fluorescent slurry was recorded by 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 using the Levenberg-Marquardt algorithm based on the grouting pressure-flow data, 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 according to the dominant grouting parameters to dynamically correct the rock penetration grouting scheme and obtain the optimal rock penetration grouting scheme. The optimal discrete fracture model was obtained based on the optimal rock penetration grouting scheme combined with the virtual discrete fracture model.

Claims

1. A method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL, characterized in that: The following steps are involved: S1: Conduct geological exploration of rock mass, and extract and process fracture geometry parameters of rock mass; fracture geometry parameters include fracture length, fracture tendency, fracture inclination, fracture opening, fracture filling and fracture profile; S2: Fit and correct each fracture geometry 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 to be imported 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: Perform physical verification, inversion analysis and grouting optimization on the virtual discrete fracture model to obtain the optimal discrete fracture model.

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; perform borehole photography on 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 surface of the rock mass. S13: Digitize the crack sketch of the exposed surface of the rock mass to extract the crack length, crack tendency and crack dip angle; 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 surface point cloud data to obtain the crack contour; S14: Use the log-normal 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: according to the crack length distribution information, the crack opening distribution information and the crack angular distribution information, each preliminary probability distribution function is optimized and corrected by using a data association algorithm to obtain a number of secondary probability distribution functions; S23: using a fusion algorithm to combine an exponential distribution function or a Weibull distribution function, to fuse the quadratic probability distribution functions, and 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 permeation 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: According to 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 fracture geometric parameters to determine the parameter range, sample number N and sampling interval ; S303: Determine the extraction probability of each fracture geometric parameter according to 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 in each sub-interval , according to the crack parameter sequence code Using the inverse probability distribution function: , calculate the corresponding crack parameter sampling value ; 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 of S307: Classify the fracture morphology into simple fractures and complex fractures according to the complexity of the fracture morphology; S308: For simple fractures: using fracture generation algorithm to randomly generate simple fracture surfaces according to 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 interfaces between fracture surfaces and the rock matrix, and between fracture surfaces, and thus completing the construction of a discrete fracture network model; S310: Verify, analyze and optimize the discrete fracture network model: compare the mean, variance and quantile of the 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.

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 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; S52: The Navier-Stokes two-phase flow module of COMSOL simulation software is used to simulate the infiltration and diffusion process of slurry in discrete fracture networks. At the same time, the fluid-solid coupling process is determined according to the fracture aperture, thereby completing 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 To set up, represents the original pressure, represents 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; 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 permeation grouting model containing a complex fracture network based on COMSOL according to claim 1, characterized in that: The step S6 includes the following steps: S61: making a resin-based fracture model according to the virtual discrete fracture model, injecting fluorescent slurry, and calculating grouting pressure-flow data; 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; 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; 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 screen the dominant grouting parameters; the dominant grouting parameters include fracture aperture and slurry viscosity; S65: Dynamically correct the rock mass penetration grouting scheme using the monitoring correction model according to the dominant grouting parameters to obtain the optimal rock mass penetration grouting scheme; S66: Based on the optimal rock mass penetration grouting scheme combined with the virtual discrete fracture model, the optimal discrete fracture model is obtained.

7. The method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL according to claim 6, characterized in that: The monitoring and 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; then uses the iterative training layer to perform several iterative training on the data feature matrix to generate a data feature weight matrix; then, through the identification output layer, the data feature weight matrix is ​​updated and optimized according to the probability feature matrix and the global loss function to obtain the optimal weight matrix; finally, through the judgment and correction layer, the data feature matrix is ​​decomposed and reconstructed using the reconstruction decoder, and the optimal distribution probability function is corrected according to the real-time monitoring results of the optimal weight matrix.

8. The method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL according to claim 7, characterized in that: The iterative training layer includes 2 Convolutional layers, 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 concatenated 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 Upsampling layer, 2 Downsampling layer, 2 average pooling layers and 1 Mish activation function; the second training branch connects 2 layers in series Upsampling layer, 2 Downsampling layer, 2 InceptionV4 blocks, 2 average pooling layers and 1 Mish activation function.

9. The method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL according to claim 7, characterized in that: 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: ， 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 iterative training loss function for each layer, Represents the preset optimization weight coefficient, Indicates the optimization efficiency of each layer’s iterative training.

10. The method for constructing a rock mass permeation grouting model containing a complex fracture network based on COMSOL according to claim 7, characterized in that: The judgment correction layer first uses a reconstruction decoder combined with a continuity equation and a 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; , 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; 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; Finally, the judgment correction layer corrects the optimal distribution probability function according to the grouting sequence.

Images