Three-dimensional numerical simulation method for grouting in fractured rock mass considering ground stress

Abstract

This invention presents a numerical simulation method for grouting in three-dimensional fractured rock masses considering the influence of geostress. First, a geometric model of the fractured rock mass is established using a three-dimensional DFN (Digital Flow Network), and this model is meshed to generate a three-dimensional fracture network. Next, the fluid flow within the three-dimensional fracture network is equivalently represented using the unified pipe network method, and the fluid transport control equations are discretized. Then, considering the influence of geostress, the fracture aperture is calculated. Finally, the pipes in the three-dimensional fracture network are divided into water flow pipes, grout flow pipes, and mixed flow pipes, and fluid flow models for each type of pipe are established. The fluid flow models for the corresponding pipes are substituted into the fluid transport control equations, and the equations are solved to obtain the head at each node in the three-dimensional fracture network, enabling the tracking of grout transport location. This method fully combines the two-phase flow of grout in saturated fractures with the unified pipe network method, while also considering the influence of magma-rock coupling in the geostress field on the fracture aperture, thus more realistically simulating the grout transport process in underground rock masses.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation technology for grouting in fractured rock masses, specifically a three-dimensional numerical simulation method for grouting in fractured rock masses that considers the influence of geostress. Background Technology

[0002] Underground tunnel, air defense, mining, and seepage prevention projects all face challenges such as extensive rock fissures, insufficient rock strength, and high permeability. Grouting is an effective method to solve these problems. In actual engineering, because the grout is injected into the fractured rock mass, it is difficult to visually observe the migration and distribution of the grout, resulting in the inherent concealment and uncertainty of grouting projects. Numerical simulation, however, can visually display the grout migration through simulation, making it an important tool for optimizing the design of grouting schemes.

[0003] Scholars have conducted a series of studies on numerical simulation of grouting in fractured rock mass. For example, the paper "Study on the diffusion mechanism of fracture grouting based on transfer matrix algorithm" uses the software COMSOL to establish a single fracture and multiple sets of parallel fracture models, introduces the transfer matrix algorithm into the Bingham body constitutive equation, and studies the grout diffusion under the condition of small changes in fracture aperture caused by grouting pressure. See Liu Ruokun et al. Study on the diffusion mechanism of fracture grouting based on transfer matrix algorithm [J]. Journal of Hydraulic and Architectural Engineering, 2021, 19(1):207-213. The paper "A two-dimensional grouting model considering hydromechanical coupling and fracturing for fractured rock mass" establishes a two-dimensional grouting model considering magma-rock coupling and fracture initiation based on the FDEM method to study the grout diffusion law. See Chengzeng Yan, et al. A two-dimensional grouting model considering hydromechanical coupling and fracturing for fractured rock mass [J]. Engineering Analysis with Boundary Elements, 2021, 133:385-397. The literature “Numerical simulation of grouting of three-dimensional random fracture rock mass in bedrock of hydropower dam” uses 3DEC to simulate the diffusion process of multi-level water-cement ratio grout in a three-dimensional random fracture network model where the fracture aperture follows a random distribution. For details, see Wang Xiaoling et al. Numerical simulation of grouting of three-dimensional random fracture rock mass in bedrock of hydropower dam [J]. Engineering Mechanics, 2018, 35(1):148-159.

[0004] As can be seen from the above, although existing technologies have studied numerical simulation methods for grouting in single fractures, two-dimensional fracture networks, and three-dimensional fracture networks, research on grouting problems based on single fractures mostly employs the finite element method, using software such as COMSOL for simulation, to study the influence of different grouting parameters, fracture aperture, fracture roughness, or fracture dip angle on grout transport within the fracture. Grouting problems in two-dimensional fracture networks mostly employ the discrete element method, using commercial software UDEC for simulation, but this method cannot be used to analyze the grout flow characteristics within fracture networks. While a few commercial software programs, such as 3DEC, can achieve numerical simulation of grouting under three-dimensional conditions, they cannot take into account the influence of in-situ stress at different locations within the fracture on grout diffusion in real time, nor can they analyze the fluid flow state at different locations within the fracture. In practical engineering, there is a geostress field in underground rock mass, which has a significant impact on the migration and diffusion of grout. However, existing numerical simulation methods for grouting in three-dimensional fractured rock mass fail to realize the influence of geostress on grout diffusion. Therefore, this application proposes a numerical simulation method for grouting in three-dimensional fractured rock mass that considers the influence of geostress. Summary of the Invention

[0005] To address the shortcomings of existing technologies, the technical problem this invention aims to solve is to provide a three-dimensional numerical simulation method for grouting in fractured rock masses that takes into account the influence of geostress.

[0006] The technical solution adopted by the present invention to solve the aforementioned technical problem is as follows:

[0007] A numerical simulation method for grouting in three-dimensional fractured rock mass considering the influence of in-situ stress, characterized by comprising the following steps:

[0008] The first step is to use a 3D DFN to build a geometric model of the fractured rock mass in Matlab software, and then mesh the geometric model to generate a 3D fracture network.

[0009] The second step is to perform equivalent analysis of the fluid flow in the three-dimensional fracture network based on the unified pipeline network method, and to discretize the fluid transport control equation of equation (1).

[0010]

[0011] In the formula, [K] is the conduction coefficient matrix, and [Q] is the pipe flow matrix. Let n be the flow potential matrix, and n be the number of nodes in the three-dimensional fracture network.

[0012] The third step is to calculate the fracture aperture, taking into account the influence of geostress.

[0013] Calculate the stress distribution around the grouting hole in the cylindrical coordinate system based on the stress distribution state within the wellbore coordinate system. Where σ represents the normal stress component, τ represents the shear stress component, the first subscript of each stress component indicates the direction of the outward normal of the surface, the second subscript indicates the direction of the force, c represents the cylindrical coordinate system, r is the radial distance from the point in space to the center of the grouting hole, θ is the angle, and z represents the height;

[0014] By rotating the stress distribution around the grouting hole using equation (3), the stress distribution σ on the crack surface is obtained. wcs ;

[0015]

[0016] In the formula, W represents the rotation matrix from the global coordinate system to the weak bedding plane coordinate system, B represents the rotation matrix from the global coordinate system to the wellbore coordinate system, C represents the rotation matrix from the wellbore coordinate system to the cylindrical coordinate system, T represents the matrix transpose, and w represents the weak bedding plane coordinate system.

[0017] e h =e n -v s (4)

[0018] e n =e r -(e0-e r )exp(-ησ' n (5)

[0019] η = 1 / [K n0 (e0-e r (6)

[0020]

[0021]

[0022] In the formula, e h e represents the fracture aperture. n v represents the deformation of the crack surface under normal stress. s e is the deformation of the fracture surface under tangential stress; e0 is the initial fracture aperture, e r Where η is the residual crack aperture, and η is the stress-aperture correlation coefficient. Where P is the effective normal stress, and K is the fluid pressure. n0 For the initial normal stiffness, u s This is the shear displacement. For the shear expansion angle, u p =τ p / K s This represents the peak shear displacement. This represents the peak shear stress. Let K be the friction angle. s For crack shear stiffness, This represents the shear stress on the crack surface;

[0023] The fourth step is to divide the pipes in the three-dimensional fracture network into water flow pipes, slurry flow pipes and mixed flow pipes, and establish fluid flow models for each type of pipe.

[0024] The fluid flow model for a mixed-flow pipe is as follows:

[0025]

[0026] In the formula, V=U(1 / β-1)(δ1-δ(x f )), ν w ν g γ represents the kinematic viscosity of water and slurry, respectively. w γ g Here, ρ represents the unit weight of water and the unit weight of slurry, respectively; τ0 is the yield stress of the slurry; and ρ is the specific gravity of the slurry. g Let g be the density of the slurry, l be the acceleration due to gravity, b be the pipe length, and J be the crack width. g =(h g,j -h g,i ) / l represents the pressure gradient in the slurry pipeline, h g,i h g,j Let ξ be the head at both ends of the slurry flow pipe, ξ be the ratio of the length of the slurry flowing into the pipe to the total length of the pipe, and δ(x) be the head at both ends of the slurry flow pipe. f () indicates the position of the leading edge of the slurry, x f The elevation at h U,i h U,j For the head at both ends of the mixed flow pipe, q U,i q U,j Let δ1 be the flow rate at both ends of the mixed flow pipeline, and δ2 and δ1 be the elevations at both ends of the pipeline.

[0027] If ξ is taken as 0 in the water pipe, then δ(x) f )=δ1, and from equation (9), the fluid flow model of the water pipe is obtained as follows:

[0028]

[0029] Where: h w,i h w,j Let q be the head at both ends of the water pipe. w,i q w,j The flow rate at both ends of the water pipe;

[0030] If ξ is taken as 1 in the slurry pipe, then δ(x) f )=δ2, and from equation (9), the fluid flow model of the slurry flow pipe is obtained as follows:

[0031]

[0032] In the formula: N=M(1-1 / β)(δ2-δ1), q g,i q g,j The flow rate at both ends of the slurry flow pipe;

[0033] The coefficient matrix in equation (9) The conduction coefficient matrix [K] n×n submatrix, The flow potential matrix submatrix, This is a submatrix of the pipeline flow matrix [Q].

[0034] Substitute the fluid flow model of the corresponding pipeline into the fluid transport control equation of equation (1), solve the equation, and obtain the water head at each node in the three-dimensional fracture network to realize the tracking of the slurry transport position.

[0035] Furthermore, in the third step, the expressions for each stress component in the stress distribution around the grouting hole in cylindrical coordinates are as follows:

[0036]

[0037] In the formula, b represents the coordinate system of the wellbore, v is Poisson's ratio, and P... m The weight of the slurry is given by ψ, where ψ is the Biot coefficient and P is the weight of the slurry. p This refers to the formation pore pressure.

[0038] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0039] 1. This invention fully combines the two-phase flow of slurry in saturated fractures with the unified pipe network method. First, based on the unified pipe network method, the fluid flow in the three-dimensional fracture network is equivalently represented, and the fluid transport control equation is discretized. The fracture pipes are divided into three types: water flow pipes, slurry pipes, and mixed flow pipes. The fluid flow models for different pipes are derived separately. The fluid flow models are equivalent to the fluid transport control equations. Therefore, the coefficient matrices of the fluid flow models for each pipe are substituted into the fluid transport control equations to solve the control equations, thereby tracking the slurry transport position and realizing a visual simulation of the slurry diffusion process.

[0040] 2. This invention considers the influence of magma-rock coupling on fracture aperture in the geostress field, transforms geostress onto the fracture surface, and comprehensively considers the influence of changes in fracture aperture at different nodes of the pipeline network on grout diffusion under the action of the normal component, shear component of geostress, and fluid pressure. This can more realistically simulate the diffusion of grout in underground rock mass and provide more valuable theoretical guidance for grouting projects. Attached Figure Description

[0041] Figure 1 This is an overall flowchart of the present invention;

[0042] Figure 2 This is a schematic diagram of a three-dimensional fracture network;

[0043] Figure 3 This is an equivalent schematic diagram of the unified pipe network method for triangular elements.

[0044] Figure 4 This is an equivalent schematic diagram of the unified pipeline method for tetrahedral elements.

[0045] Figure 5 A schematic diagram of fluid flow in a fracture network;

[0046] Figure 6 A schematic diagram showing the division of different pipelines;

[0047] Figure 7 This is a three-dimensional distribution map of grout diffusion in fractured rock masses. Detailed Implementation

[0048] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments, but this does not limit the scope of protection of this application.

[0049] This invention relates to a three-dimensional numerical simulation method for grouting in fractured rock masses that considers the influence of in-situ stress (hereinafter referred to as the method, see [link]). Figures 1-6 This includes the following steps:

[0050] Step 1: Generate a three-dimensional fracture network;

[0051] A geometric model of the fractured rock mass was established using a 3D Discrete Facture Network (DFN) in the commercial software Matlab. The geometric model was then meshed to obtain a 3D fracture network. The geometric model includes a rock matrix and random fractures generated within a defined region using the Monte Carlo method. The rock matrix was meshed using tetrahedral elements, and the fracture surfaces were meshed using triangular elements. The upper and lower limits of the mesh elements were 0.3 m and 0.15 m, respectively. The 3D fracture network can be found in [reference needed]. Figure 2 ;

[0052] The second step involves applying the unified pipe network method to the equivalent flow of fluid in a three-dimensional fractured network, treating the edges of the mesh elements as pipes, and discretizing the following fluid transport control equation:

[0053]

[0054] In the formula, [K] is the conduction coefficient matrix, n is the number of nodes in the three-dimensional fracture network, and [Q] is the pipe flow matrix. Here is the flow potential matrix. Here, p represents the fluid pressure, ρ represents the slurry density, g represents the gravitational acceleration, and h represents the slurry head height.

[0055] The unified pipeline method is a numerical method that treats fractures and rock matrix in fractured rock masses as equivalent one-dimensional pipelines, assuming fluid flows within these pipelines, and then uses a unified approach to treat the fractures and rock matrix. In other words, it treats the edges of triangular and tetrahedral elements as equivalent pipelines. Figure 3 The diagram illustrates the equivalent principle of fluid flow within a triangular element ijm. Each vertex of the triangular element represents a node in the mesh, and the three sides mi, ij, and jm are the equivalent pipes. Point o is the circumcenter of the triangular element, and oe, of, and og are the perpendicular bisectors from the circumcenter o to the three sides mi, ij, and jm, respectively. Using the unified pipe network method, the fluid flow between regions ifoe and jgof is equivalent to flow within pipe ij, with the flow rate denoted as Q. of The flow of fluid between regions ifoe and mgoe is equivalent to the flow in pipe mi, and the flow rate is denoted as Q. oe Similarly, the flow of fluid between regions mgoe and jgof is equivalent to the flow in pipe jm, with the flow rate denoted as Q. og ; Conductivity within pipe ij l of The length of the perpendicular bisector of 'of' is given by 'k', where 'k' represents permeability, 'μ' is the hydrodynamic viscosity, and 'l' is the length of the perpendicular bisector of 'of'. ij Let be the length of pipe ij; Figure 4 The diagram illustrates the equivalent principle of fluid flow within a tetrahedral element uijm. The edges of the tetrahedral element represent the equivalent pipes. Point o is the circumcenter of the tetrahedral element, f is the midpoint of edge ij, and c1, c2, c3, and c4 are the circumcenters of the four faces. Connecting the circumcenters of the tetrahedral element, the midpoints of each edge, and the circumcenters of each face sequentially yields six dividing faces. These six dividing faces divide the tetrahedral element into four parts. Similar to the triangular element, the unified pipe network method equates the fluid flow between two regions corresponding to the same edge to the flow within the corresponding pipe. Thus, the fluid flow within the tetrahedral element is equivalent to the flow within the six pipes of the tetrahedral element. The conductivity within pipe ij is... Where A is the contact area between two adjacent regions, that is, the area of ​​the corresponding dividing surface.

[0056] Step 3: Calculate the crack aperture;

[0057] Calculate the stress distribution around the grouting hole in the cylindrical coordinate system based on the stress distribution state within the wellbore coordinate system. The expressions for each stress component are as follows:

[0058]

[0059] In the formula, σ represents the normal stress component, τ represents the shear stress component, the first subscript of each stress component indicates the direction of the outward normal to the surface, the second subscript indicates the direction of the force, the superscript of the stress component indicates the coordinate system, i.e., c represents the cylindrical coordinate system, b represents the wellbore coordinate system, a is the radius of the grouting hole, r is the radial distance from the point in space to the center of the grouting hole, θ is the angle, v is Poisson's ratio, and P... m The weight of the slurry is given by ψ, where ψ is the Biot coefficient and P is the weight of the slurry. p The pore pressure is the formation pressure, and z represents the height.

[0060] The stress distribution around the grouting hole is transformed by the rotation of equation (3) to transfer the influence of the in-situ stress to the fracture surface, thus obtaining the stress distribution σ on the fracture surface. wcs ;

[0061]

[0062] In the formula, W represents the rotation matrix from the global coordinate system to the weak bedding plane coordinate system, B represents the rotation matrix from the global coordinate system to the wellbore coordinate system, C represents the rotation matrix from the wellbore coordinate system to the cylindrical coordinate system, T represents the matrix transpose, and w represents the weak bedding plane coordinate system.

[0063] The fracture aperture is calculated according to the following formulas (4)-(8):

[0064] e h =e n -v s (4)

[0065] e n =e r -(e0-e r )exp(-ησ' n (5)

[0066] η = 1 / [K n0 (e0-e r (6)

[0067]

[0068]

[0069] In the formula, e h e represents the fracture aperture. n v represents the deformation of the crack surface under normal stress. s e is the deformation of the fracture surface under tangential stress; e0 is the initial fracture aperture, e rWhere η is the residual crack aperture, and η is the stress-aperture correlation coefficient. Where P is the effective normal stress, and K is the fluid pressure. n0 For the initial normal stiffness, u s This is the shear displacement. For the shear expansion angle, u p =τ p / K s This represents the peak shear displacement. This represents the peak shear stress. Let K be the friction angle. s For crack shear stiffness, This represents the shear stress on the crack surface.

[0070] Step 4: Divide the channels in the three-dimensional fracture network into water flow channels, slurry flow channels, and mixed flow channels. Channels filled with slurry are slurry channels, channels not filled with slurry are mixed channels, and channels where slurry has not entered are water flow channels; for example... Figure 5 , 6 As shown, the slurry flows in a triangular element. The leading edge of the slurry passes through point m. The leading edge of the slurry divides the pipes of the triangular element into different types of pipes: im is the slurry pipe, mj is the water flow pipe, and ij is the mixed flow pipe.

[0071] Establish fluid flow models for various pipelines;

[0072] The fluid flow model for a mixed-flow pipe is as follows:

[0073]

[0074] In the formula, V=U(1 / β-1)(δ1-δ(x f )), ν w ν g γ represents the kinematic viscosity of water and slurry, respectively. w γ g Here, ρ represents the unit weight of water and the unit weight of slurry, respectively; τ0 is the yield stress of the slurry; and ρ is the specific gravity of the slurry. g Let g be the density of the slurry, l be the acceleration due to gravity, b be the pipe length, and J be the crack width. g =(h g,j -h g,i ) / l represents the pressure gradient in the slurry pipeline, h g,i h g,j Let ξ be the head at both ends of the slurry flow pipe, ξ be the ratio of the length of the slurry flowing into the pipe to the total length of the pipe, and δ(x) be the head at both ends of the slurry flow pipe. f () indicates the position of the leading edge of the slurry, x f The elevation at h U,i hU,j For the head at both ends of the mixed flow pipe, q U,i q U,j Let ξ be the flow rate at both ends of the mixed flow pipe; if ξ is taken as 0 in the water flow pipe, then δ(x) f )=δ1, therefore equation (9) can be written as equation (10); in the slurry pipe, ξ is taken as 1, then δ(x) f )=δ2, so equation (9) can be written as equation (11), where δ2 and δ1 are the elevations at both ends of the pipe;

[0075] The fluid flow model for a water pipe is as follows:

[0076]

[0077] Where: h w,i h w,j Let q be the head at both ends of the water pipe. w,i q w,j The flow rate at both ends of the water pipe;

[0078] The fluid flow model for the slurry flow pipe is as follows:

[0079]

[0080] In the formula: N=M(1-1 / β)(δ2-δ1), q g,i q g,j The flow rate at both ends of the slurry flow pipe;

[0081] Equation (9) is equivalent to Equation (1), and the coefficient matrix in Equation (9) is... The conduction coefficient matrix [K] n×n submatrix, The flow potential matrix submatrix, It is a submatrix of the pipeline flow matrix [Q]; therefore, according to the pipeline type, the fluid flow model of the corresponding pipeline is substituted into the fluid transport control equation of equation (1), the equation is solved, and the head at each node in the three-dimensional fracture network is obtained, so as to realize the tracking of the slurry transport position.

[0082] Fluid migration within a fracture network involves changes in fluid pressure due to diffusion. These changes in pressure affect fracture aperture, which in turn influences fluid diffusion, resulting in magma-rock coupling during fluid flow. Iterative analysis of fracture aperture and fluid flow can provide real-time tracking of magma migration location. Figure 7 This is a diagram showing the final grout diffusion distribution within a three-dimensional fractured rock mass.

[0083] Any aspects not covered in this invention are applicable to existing technologies.

Claims

1. A numerical simulation method for grouting in three-dimensional fractured rock mass considering the influence of geostress, characterized in that, Includes the following steps: The first step is to use a 3D DFN to build a geometric model of the fractured rock mass in Matlab software, and then mesh the geometric model to generate a 3D fracture network. The second step is to perform equivalent analysis of the fluid flow in the three-dimensional fracture network based on the unified pipeline network method, and to discretize the fluid transport control equation of equation (1). In the formula, [K] is the conduction coefficient matrix, and [Q] is the pipe flow matrix. Let n be the flow potential matrix, and n be the number of nodes in the three-dimensional fracture network. The third step is to calculate the fracture aperture, taking into account the influence of geostress. Calculate the stress distribution around the grouting hole in the cylindrical coordinate system based on the stress distribution state within the wellbore coordinate system. Where σ represents the normal stress component, The shear stress component is represented by the first subscript of each stress component, which indicates the direction of the outward normal of the surface on which it acts, and the second subscript indicates the direction of the force. c represents the cylindrical coordinate system, r is the radial distance from the point in space to the center of the grouting hole, θ is the angle, and z represents the height. By rotating the stress distribution around the grouting hole using equation (3), the stress distribution σ on the crack surface is obtained. wcs ; In the formula, W represents the rotation matrix from the global coordinate system to the weak bedding plane coordinate system, B represents the rotation matrix from the global coordinate system to the wellbore coordinate system, C represents the rotation matrix from the wellbore coordinate system to the cylindrical coordinate system, T represents the matrix transpose, and w represents the weak bedding plane coordinate system. yes h =e n -v s (4) e n =e r -(e0-e r )exp(-ησ' n ) (5) η=1 / [K n0 (e0–e r )] (6) In the formula, e h e represents the fracture aperture. n v represents the deformation of the crack surface under normal stress. s e is the deformation of the fracture surface under tangential stress; e0 is the initial fracture aperture, e r Where η is the residual crack aperture, and η is the stress-aperture correlation coefficient. Where P is the effective normal stress, and K is the fluid pressure. n0 For the initial normal stiffness, u s This is the shear displacement. For the shear expansion angle, u p =τ p / K s This represents the peak shear displacement. This represents the peak shear stress. Let K be the friction angle. s For crack shear stiffness, This represents the shear stress on the crack surface; The fourth step is to divide the pipes in the three-dimensional fracture network into water flow pipes, slurry flow pipes and mixed flow pipes, and establish fluid flow models for each type of pipe. The fluid flow model for a mixed-flow pipe is as follows: In the formula, V=U(1 / β-1)(δ1-δ(x f )), ν w ν g γ represents the kinematic viscosity of water and slurry, respectively. w γ g Here, ρ represents the unit weight of water and the unit weight of slurry, respectively; τ0 is the yield stress of the slurry; and ρ is the specific gravity of the slurry. g Let g be the density of the slurry, and g be the acceleration due to gravity. where b is the pipe length, b is the crack width, and J is the crack width. g =(h g,j -h g,i ) / h represents the pressure gradient in the slurry pipeline. g,i h g,j Let ξ be the head at both ends of the slurry flow pipe, ξ be the ratio of the length of the slurry flowing into the pipe to the total length of the pipe, and δ(x) be the head at both ends of the slurry flow pipe. f () indicates the position of the leading edge of the slurry, x f The elevation at h U,i h U,j For the head at both ends of the mixed flow pipe, q U,i q U,j Let δ1 be the flow rate at both ends of the mixed flow pipeline, and δ2 and δ1 be the elevations at both ends of the pipeline. If ξ is taken as 0 in the water pipe, then δ(x) f )=δ1, and from equation (9), the fluid flow model of the water pipe is obtained as follows: Where: h w,i h w,j Let q be the head at both ends of the water pipe. w,i q w,j The flow rate at both ends of the water pipe; If ξ is taken as 1 in the slurry pipe, then δ(x) f )=δ2, and from equation (9), the fluid flow model of the slurry flow pipe is obtained as follows: In the formula: N=M(1-1 / β)(δ2-δ1), q g,i q g,j The flow rate at both ends of the slurry flow pipe; The coefficient matrix in equation (9) The conduction coefficient matrix [K] n×n submatrix, The flow potential matrix submatrix, This is a submatrix of the pipeline flow matrix [Q]. Substitute the fluid flow model of the corresponding pipeline into the fluid transport control equation of equation (1), solve the equation, and obtain the water head at each node in the three-dimensional fracture network to realize the tracking of the slurry transport position.

2. The numerical simulation method for grouting in three-dimensional fractured rock mass considering the influence of geostress as described in claim 1, characterized in that, In the third step, the expressions for each stress component in the stress distribution around the grouting hole in cylindrical coordinates are as follows: In the formula, b represents the coordinate system of the wellbore, v is Poisson's ratio, and P... m The weight of the slurry is given by ψ, where ψ is the Biot coefficient and P is the weight of the slurry. p This refers to the formation pore pressure.

Images