Grouting Process Simulation Based on 3D Fracture Network Considering Fluid–Structure Interaction

Abstract

Grouting has always been the main engineering measure of ground improvement and foundation remediation of hydraulic structures. Due to complex geological conditions and the interactions between the grout and the fractured rock mass, which poses a serious challenge to the grouting diffusion mechanism analysis, fracture grouting has been a research hotspot for a long time. In order to throw light on the grout diffusion process in the fractured rock mass and the influence of grout on the fracture network, and to achieve more realistic grouting numerical simulation, in this paper a grouting process simulation approach considering fluid–structure interaction is developed based on the 3D fractured network model. Firstly, the relationship between fracture apertures and trace lengths is used to obtain a more realistic value of fracture aperture; then a more reliable model is established; subsequently, based on the 3D fracture network model, different numerical models are established to calculate fluid dynamics (grout) and structure deformation (fractured rock mass), and the results are exchanged at the fluid–structure interface to realize the grouting process simulation using two-way fluid-structure interaction method. Finally, the approach is applied to analyze the grouting performance of a hydropower station X, and the results show that the grouting simulation considering fluid–structure interaction are more realistic and can simultaneously reveal the diffusion of grout and the deformation of fracture, which indicates that it is necessary to consider the effect of fluid–structure interaction in grouting simulation. The results can provide more valuable information for grouting construction.

1. Introduction

Many high dams are built in areas with complicated geological conditions, and the numerous fractures and voids will increase the permeability and decrease the strength of the rock mass, which affect the stability and safety of the dam’s foundations. As a common and effective measure to improve the geological conditions of dam foundations, grouting is used to fill up the joints and fractures in the rock mass so as to prevent seepage and improve the bearing capacity and deformation resistance [1,2,3,4]. However, fracture grouting is still a difficult issue due to complicated fracture distribution, complex fluid–structure interaction effects, and incomplete information on grout diffusion behavior and corresponding fracture deformation. In order to reveal the grouting mechanisms in the fractured rock mass, it is necessary to select an effective tool for studying the grouting process, especially for considering the fluid–structure interaction.

Computational fluid dynamics (CFD) is an effective tool used to simulate fracture grouting, which can partly overcome some limitations of experiments. In recent years, various researchers have carried out an abundance of research. Saeidi et al. [5] established a numerical model to study the effect of fracture properties on grout flow and penetration length in fractured rock mass using Universal Distinct Element Code (UDEC). Yang et al. [6] simulated the cement grout diffusion process in a single rough fracture by the finite element method. Fu et al. [7] performed numerical simulation on the diffusion process of cement grouting in the fractures of the rock mass to determine reasonable hole spacing and other parameters. Hao et al. [8] developed a numerical simulation of polymer grout diffusion in a single fracture to analyze the pressure distribution. Deng et al. [4] proposed a CFD simulation approach based on 3D fracture network model to study the grouting process of a dam’s foundation. Kim et al. [9] used UDEC to simulate the flow of Bingham grout in a single joint with smooth parallel surfaces and considered the hydromechanical coupling to study its effect on grouting performance. Ao et al. [10] simulated the grouting process in underground goaf and analyzed the stability by applying one-way fluid–structure interaction. Liu et al. [11] combined a finite-discrete element method (FDEM) and a grouting flow simulator to consider the hydromechanical coupling effect in the parallel-plate model. The review of previous studies reveals that most of the aforementioned approaches simulated the grouting process in a single fracture, 2D fracture network or simplified rock mass, which cannot reflect the actual diffusion of grout flow in complex fractured rock mass completely. Furthermore, these studies only studied the diffusion of grout in fractures or the effect of grout on the rock mass, without taking the complex fluid–structure interaction between grout flow and rock mass into account, which were different from the actual conditions.

In order to obtain a more realistic simulation of grouting process, the establishment of a precise and reliable three-dimensional fracture network model is an important prerequisite [4]. Since deterministic data of the fracture aperture are not available, in current studies of fracture network modeling and its application the fracture aperture is usually ignored [12,13,14], reduced to a given value [15], or randomly generated from a given range or geologically conditioned statistical distributions [4,16]. The fracture aperture significant influence the permeability of the fractured rock masses, so it is necessary to establish a fracture network model with a more authentic and accurate fracture aperture to simulate the grouting process. For fractured rock mass, fractures are usually random and complex which are fractal structure with self-similarity; according to this characteristic of a fracture, some scholars have put forward the formula for the relationship between fracture apertures and trace lengths, and this relationship has been widely investigated [17,18,19,20]. In this study, this relationship will be introduced into the fracture network modeling to make up for the deficiency of the existing research.

As the pressurized grout penetrates the fractures inside the rock mass, the grout will separate the fracture surfaces from each other, causing an interaction between the grout and the rock mass [21]. Tsang et al. [22] indicated the coupling of processes implies that the both interact in the initiation and progress of each other. On the one hand, the grout pressure induces stresses on the surfaces of the fracture, this will lead to the deformation of the fractures. On the other hand, the change of the fracture aperture will affect the fracture permeability, and then results in the variation of grout performance. So the grouting performance cannot be determined by considering each process independently. Some theoretical studies on the interaction between grout and rock mass have been carried out. GothäLl and Stille [21] analyzed the interaction of two parallel fracture during high pressure grouting and discussed the effect of fracture dilation on the penetrability of fine fractures. Rafi and Stille [23] proposed a procedure for optimizing grouting pressure based on the estimation of grout spread and the identification of jacking of the fracture. Rafi and Stille [24] described the basic mechanism of elastic deformation during grouting and discussed its impact on the spread of grout. In literature [9], the authors strongly recommended that the interaction between the grout flow and fractured rock mass should be included in the grouting analysis in order to have a precise prediction of grout performance. With this in mind, the fluid–structure interaction (FSI) method as a numerical technique is used to solve problems that involve the mutual interaction of fluid and structure. In recent years, with the development of computer performance and increasing interest in more realistic modeling, FSI has attracted extensive attention in the computational field [25]. For instance, this technique has been applied in the problem of internal slip during the operation of progressive cavity pump (PCP)s in oilfield production [26], biomedical problems where blood flow interacts with blood vessel walls [27], and the fluid–structure interaction problem of fracturing structures under impulsive loads [28]. Of interest from a fractured grouting perspective, there is still much work to do in the grouting simulation considering fluid–structure interaction to capture the interaction between grout flow characteristics and deformation of the fractures.

In summary, most of the existing studies simulated the grouting process in a single fracture, 2D fracture network or simplified rock mass, which is inconsistent with actual fractures under complex geological conditions. Moreover, the value of fracture aperture is usually ignored or inaccurate in the fracture network modeling, which will affect the authenticity of grouting simulation. Additionally, rich theoretical research achievements have proposed on the interaction between grout and fracture, and some of the grouting simulations considering fluid-structure interaction are based on single-fracture or just using one-way fluid–solid coupling. Therefore, based on the 3D fractured network model the numerical simulation of the grouting process considering two-way fluid-structure interaction still needs further study.

In this study, a grouting simulation approach considering fluid–structure interaction is developed based on the 3D fractured network model. Firstly, fracture parameters are randomly simulated by the Latin hypercube sampling (LHS) method based on the statistical information from fracture survey and borehole imaging of the exposed surface, the relationship between fracture apertures and trace lengths is used to obtain the value of fracture aperture, then a more reliable 3D fracture network model for dam foundation rock mass is established with VisualGeo software [29]. Next, the CFD simulation model of grout (fluid) and the finite element model of fractured rock mass (structure) are established respectively, and their governing equations are solved in different ways, with the results exchanged through the fluid–structure interface to realize the two-way fluid-structure interaction simulation of the grouting process. Finally, the approach is used in a case study to analyze the dam foundation grouting to investigate the effects of fluid-structure interaction on grouting processes; the results show that the grouting simulation considering fluid–structure interaction are more realistic and can simultaneously reveal the grout diffusion and fracture deformation under the interaction of grout and rock mass, which can provide more valuable information for optimizing the grouting process.

The remaining parts of this paper are organized as follows: the methodology of 3D fracture network modeling and fluid–structure interaction simulation are introduced in Section 2. In Section 3, the approach is applied to analysis of the grouting performance of hydropower station X and the studies on the grouting characteristics are given in this section. Finally, the conclusions are provided in Section 4.

2. Methodology

2.1. Modeling of 3D Fracture Network

2.1.1. Modeling Process of 3D Fracture Network

Due to the large amount of complex fractures in the rock mass of a dam’s foundations, it is difficult to determine the exact position and occurrence of each fracture by using a deterministic model. A large number of engineering practices and studies have shown that the fractures have obvious statistical distribution rules and characteristics, so we established the 3D fracture network model which is close to the real fracture conditions in a statistical sense.

The modeling process (Figure 1) mainly includes the following steps: (1) the statistical homogeneous zone is divided firstly, then the fractures in the statistical homogeneous zone are divided into dominant sets and the cracks with similar properties are clustered; (2) the fracture space density and the distributions of the geometry parameters could be obtained based on the statistical analysis; (3) the fracturing parameters are randomly simulated by the LHS method; and (4) a 3D fracture network model is constructed in VisualGeo software.

2.1.2. Statistical Analysis of Fracture Geometric Characteristic Parameters

In this study, the fracture was simulated by the Baecher disc model [30] which assumed every fracture as a thin disc (Figure 2). The fracture disc model can be defined by the following formula:

$$C=c\left(O,V,R,A\right)$$

(1)

the formula defines a fracture disc with center point O, occurrence V, radius R and aperture A. Where O = (x₀, y₀, z₀), V = (α, β), α and β are the dip direction and dip angle of the fracture disc respectively, and n is the normal vector of the fracture disc.

The relationship between the various parameters can be expressed as follows:

$$\{\begin{array}{l}a\left(x-x_0\right)+b\left(y-y_0\right)+c\left(z-z_0\right)=0\\ {\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2\le R^2\\ a=sin\beta sin\alpha ,b=sin\beta cos\alpha ,c=cos\beta \\ A=f\left(R\right)\end{array}$$

(2)

According to the exposed surface fracture catalog data and digital borehole data, the distributions of fracture geometric characteristic parameters can be determined.

(1)

Fracture space location

The Poisson process [31] is widely used to describe fracture location. The fractures are mutually independent and the uniform distribution function are adopted to obtain the coordinates (x₀, y₀, z₀) of the fracture center point.

(2)

Fracture density

The Mauldon method [32] is adopted to estimate the fracture volume density. The following equation is used to estimate the trace area density:

$${\lambda }_a=\frac{n_1+2n_2}{2WH}$$

(3)

where λ_a is the trace area density, n₁ is the number of traces which one end can be observed, n₂ is the number of traces which both end can be observed, W is the width of rectangular window, H is the height of rectangular window. Then, Equation (4) is adopted to obtain the fracture volume density:

$$E\left({\lambda }_v\right)=\frac{E\left({\lambda }_a\right)}{E\left(D\right)E\left|sin\nu \right|}$$

(4)

where λ_v is the fracture volume density, D is the fracture diameter, sinν is the sine value of the dip.

(3)

Fracture size

To simulate the size of the fracture surface, statistical analysis of the fracture trace length is needed first. Huang et al. [33] put forward the estimation formula of trace length:

$$l=\frac{n_1+2n_0}{2N}\frac{\pi WH}{W+H}$$

(5)

where l is the fracture trace length in the window, n₀ is the number of traces which neither end can be observed, n₁ is the number of traces which one end can be observed, N is the total number of fracture traces in the window, W is the width of rectangular window, H is the height of rectangular window.

When the disc model is used to simulate the fracture, the fracture size is expressed by its diameter. the fracture diameter distribution can be confirmed based on the distribution of trace length.

(4)

Fracture occurrence

According to Kemeny and Post [34], the fisher distribution can be used to fit fracture occurrence and obtained relatively better results.

(5)

Fracture aperture

Schultz et al. [17,18] conducted a lot of statistical studies and obtained the expression of the relation between fracture aperture and fracture trace length:

$$A=\beta l^n$$

(6)

where A is the fracture aperture, l is the fracture trace length, β and n are constants related to the properties of fractured rock mass.

In this study, the value of n = 1 is chosen to reflect the self-similarity and fractal of the fracture network [35,36,37], and then the relationship between fracture aperture and fracture trace length is linear; based on this relationship, the existing survey data are linearly fitted to get the value of γ. Finally, the fracture aperture can be obtained by Equation (6) based on the data of trace length.

2.1.3. Latin Hypercube Sampling (LHS) Random Sampling

After obtaining the determined probability distribution model of each fracture parameter, random sampling of parameters is needed. The essence of the LHS method is to divide the sampling interval according to the sampling times, and then random sampling is carried out from each subinterval. This method avoids the collapse of the sample data and the simulation results are more stable. Therefore, the LHS method is adopted to randomly simulate the fracture parameters in this study.

Then, taking the center point coordinates, diameters, aperture, dip direction and dip angle of fractures simulated by LHS method as input parameters, a three-dimensional model of rock mass fracture network is established by using VisualGeo software.

2.2. Fluid–Structure Interaction Model

The fluid–structure interaction model mainly consists of two parts: computational fluid dynamics (CFD) and computational structure dynamics (CSD). The grouting process is simulated by CFD and the deformation of fractures is calculated by CSD.

The solution of fluid-structure interaction includes a directly coupled solution and partitioned solution. The first one solves the governing equations of fluid and structure simultaneously in the same solver by coupling the governing equations of fluid and structure to the same equation matrix, so its advantage is that there is no time lag problem. However, a direct coupling solution may result in poor convergence and huge computational cost. As a consequence, it is difficult to realize in fact [38]. On the contrary, the second one solves the fluid governing equations and the structure governing equations in different solvers, and the results are exchanged and transmitted through the fluid-structure interface. In this study, we choose the partitioned solution to solve the fluid–structure interaction between grout and fractures.

2.2.1. Computational Fluid Dynamics (CFD) Grouting Numerical Model

The governing equations of the grouting can be described by the continuity equation, momentum equation, two-phase volume of fluid (VOF) equation and Papanastasiou regularized equation.

(1)

The two-phase VOF equation

In the process of grouting, the grout drives out air or groundwater, which should be treated as a two-phase flow [4]. The accurate description of the interface between two kinds of incompatible and incompressible fluids is one of the most important issues in multi-fluid flow computations [39], this can be solved by the VOF method which is proposed by Hirt and Nichols [40] to track free fluid surfaces under fixed grid condition. Therefore, the VOF method is used to keep track of the grout-air interface in this paper. In this method, a volume fractional variable F = F (x, y, z, t) for each phase of the model in the computational domain is introduced. F_g = 1 indicates that the volume is occupied by grout while F_g = 0 indicates that the volume contains no grout and is in the air phase, and 0 < F_g < 1 stands for the volume that contains both grout and air. Equation (7) is used to describe the motion of the grout-air interface:

$$\frac{\partial F_g}{\partial t}+\rho \nu \nabla F_g=0$$

(7)

where F_g is the volume fraction of grout; ρ is the density of fluid in kg/m³; ν is the kinematic viscosity of fluid in m²/s.

(2)

The continuity equation:

$$\frac{\partial \rho }{\partial t}+\nabla \left(\rho u\right)=0$$

(8)

where ρ is the density of fluid in kg/m³; t is the time in s; u is the velocity of the unit section in m/s.

(3)

The momentum equation:

$$\rho \frac{du}{dt}=-\nabla p+\rho g+\nabla \left(\eta \dot{\gamma }\right)+F_{\mathrm{st}}+S$$

(9)

where p is the pressure on the fluid micro-unit in Pa, g is the acceleration of gravity in m/s²; η is the apparent viscosity of fluid in Pa·s; the relationship between ν and η is ν = η/ρ; $\dot{\gamma }$ is the shear rate in 1/s, F_st is the surface tension force in N/m³ and is presented as Equation (10) [41]; and S is the momentum resistance source term in N/m³, including inertia loss term S_i and viscosity loss term Sν; in this paper, S_i can be neglected because of the low velocity of the grout and S = S_ν. Equation (11) is the expression of the viscosity loss term S_ν:

$$F_{\mathrm{st}}=-\sigma \nabla ·\left(\frac{\nabla F_i}{\left|\nabla F_i\right|}\right)\nabla F_i$$

(10)

where F_i is the volume fraction of phases; σ is the surface tension coefficient in N/m.

$$S_{\nu }=-\frac{\rho \nu }{\alpha }u$$

(11)

where $\frac{1}{\alpha }$ is the viscous drag coefficient and its expression is as follows:

$$\frac{1}{\alpha }=\frac{g}{K\nu }$$

(12)

where K is the permeability coefficient.

In order to obtain single set of equations, ρ and ν in Equations (7)–(12) are no longer constants but are variables weighted by the volume fraction of fluid [42]:

$$\rho =F_g{\rho }_g+\left(1-F_g\right){\rho }_a$$

(13)

$$\nu =F_g{\nu }_g+\left(1-F_g\right){\nu }_a$$

(14)

where ρ_g, ρ_a, ν_g, ν_a are the density of grout, the density of air, the kinematic viscosity of grout and the kinematic viscosity of air, respectively.

(4)

The Papanastasiou regularized equation

The cement grout with a w/c ratio of less than 1 is usually described by the Bingham model. However, in the Bingham constitutive equation, when the shear rate is close to zero, the apparent viscosity will become infinite, which causes problems in numerical simulation. In order to solve this problem, the Papanastasiou regularized model is used to describe the rheological properties of cement grout [4], as shown in Equation (15):

$$\eta =\{\begin{array}{ll}m{\tau }_0& \dot{\gamma }=0\\ {\eta }_0+\frac{{\tau }_0}{\dot{\gamma }}\left[1-e^{-m\dot{\gamma }}\right]& \dot{\gamma }\ne 0\end{array}$$

(15)

where η is the apparent viscosity of fluid in Pa·s; m is the stress growth parameter in s; τ₀ is the yield stress in Pa; $\dot{\gamma }$ is the shear rate in 1/s; η₀ represents the plastic viscosity in Pa·s; and e is a natural constant. From our previous research [4], it is considered that m = 100 can meet the research needs and can make the numerical model effectively express the rheological properties of cement grout.

2.2.2. Computational Structure Dynamics (CSD) Model

The rock mass is elastoplastic materials, since the objective has been to establish a model applicable for grouting problems, the Mohr–Coulomb (M-C) shear yield criterion which is comparable to actual rock was applied to the rock mass. The expressions of the criterion are shown as follows:

$${\tau }_n=C+{\sigma }_{n}tan\varphi$$

(16)

if $\tau <{\tau }_n$, the rock mass is linear elastic material, if $\tau \ge {\tau }_n$, the rock mass become yield.

The linear elastic model is based on the generalized Hooke law, and the constitutive equation is as follows:

$$\left\{\epsilon \right\}=\left[D\right]\left\{\sigma \right\}$$

(17)

where ε is the strain, D is the elastic matrix, and σ is the stress component.

In the stress space, the form of the yield function is as follows:

$$f=\frac{1}{2}\left({\sigma }_1-{\sigma }_3\right)-\frac{1}{2}\left({\sigma }_1+{\sigma }_3\right)sin\varphi -Ccos\varphi$$

(18)

where σ₁ and σ₂ are the maximum and minimum principal stresses of material damage, respectively; C is the cohesion; and φ is the internal friction angle. When f ≥ 0, shear failure will occur in the material.

2.2.3. Fluid–Structure Interaction Analysis Solution

In this paper, Finite Element Analysis (FEA) software solves the structural domain (fractured rock mass) and CFD software solves the fluid domain (grout). The two domains are interconnected on the fluid–structure interface using the SIMULIA Co-Simulation Engine (CSE) [41] which is widely used to couple CFD software and FEA software for fluid–structure interaction simulation [43,44,45].

As show in Figure 3, in the fluid-structure interaction solving process, CFD software initializes the fluid field and passes loads to the FEA software (pressure + wall shear stress), and the deformations in the structure can be computed by FEA software and passes displacements back to CFD software; this can provide a new deformed geometry for the CFD software to solve the fluid field. This iteration can be repeated until the end of the coupling process.

On the fluid–structure interface, the stress (τ) and the displacement (d) of the fluid and structure should be equal:

$${\tau }_g{n}_g={\tau }_f{n}_f$$

(19)

$$d_g=d_f$$

(20)

where the subscript g and f represent grout and fractured rock mass.

2.2.4. Boundary Conditions

(1)

Inlet boundary conditions: according to the data of grouting pressure measured by grouting recorder and taking the mean value of grouting pressure during grouting period, pressure inlet is set at the boundary of grouting borehole interval. The corresponding grout VOF at the inlet is set to 1.

(2)

Outlet boundary conditions: the pressure outlet is set at the end boundary of the fracture, and the pressure satisfies the second boundary condition.

(3)

Initial conditions: assuming that there is no groundwater during grouting, the fractures are filled with air before grouting, and the initial air VOF in the fracture is set to 1.

(4)

Displacement boundary conditions: the bottom boundary of the computational domain is the z-axis constraint, the lateral boundaries are the x- and y-axis constraints.

3. Case Study

Hydropower station X is located in the upper reaches of Lancang River. It is a large-scale hydropower project which mainly generates electricity and takes into account the comprehensive utilization benefits of irrigation and water supply. The hydropower project is mainly composed of a gravelly soil core rockfill dam, left bank slope spillway, water diversion and power generation system, ground workshop, and so on. The installed capacity of the hydropower station is 1400 MW, the maximum dam height is 139.80 m, and total dam crest length is 576.68 m. The dam area exposed strata are mainly from the Middle Jurassic flower group (J2h) and Quaternary strata (Q). The layout of the dam’s foundation grouting curtain and the project profile are shown in Figure 4. The dam foundation curtain grouting project is divided into several continuous grouting units. In this study, a typical grouting unit is taken as a case study to simulate a three-dimensional random fracture network. The location of the study area is shown in the red wireframe in Figure 4b. This area is the foundation curtain grouting unit of the river bed dam section and is also the main area of the dam foundation seepage control of this project. In the study area, there are 41 grouting boreholes including 21 in the upstream row and 20 in the downstream row with the hole spacing of 1.5 m, the diameter of the grouting borehole is 75 mm, and the total depth is 60 m.

3.1. Simulation of 3D Fracture Network

According to the size of the exposed surface and the depth of the grouting hole, the study area is 30 m × 16 m × 60 m (length × width × depth). In the study area, a total of 83 fractures were recorded on the exposure surface, the fracture sketch picture and the fracture dominant sets diagram are shown in Figure 5. The fractures were divided into 3 sets based on the occurrence.

The distribution of fracture parameters in each set is fitted according to mathematical statistics, the results shows that the trace lengths of the three sets of fractures obey the logarithmic normal distribution and the occurrences obey the Fisher distribution. The statistical results of the fracture parameters of the dam foundation are shown in

Table 1. Statistical results of fracture parameters.

Set	Fracture Number	Regional Volume (m3)	Parameter	Mean	Minimum	Maximum	Distribution
1	2587	28,800	Coordinate X/m	15	0	30	Uniform
			Coordinate Y/m	8	0	16
			Coordinate Z/m	30	0	60
			Diameter/m	1.93	0.31	2.89	Lognormal
			Dip direction/degree	345.00	340.00	349.98	Fisher
			Dip angle/degree	15	10.00	19.98

(Taking first set fractures as an example).

According to the statistical results of fracture parameters of the dam foundation rock mass, the fracturing parameters are randomly simulated 10 times by the LHS method and the optimal results were obtained. Then, the 3D fracture network models are constructed in VisualGeo software. The effect drawing of the final 3D fracture network model for the study area is shown in Figure 6.

3.2. Grouting Simulation Considering Fluid–Structure Interaction

In this study, the second stage (3–6 m) of the grouting borehole 3-LR1-33 is selected from the established 3D fracture network model as the simulation object, the selected model is a cylindrical area centered on the grouting borehole with the radius of 1.5 m, which contains 20 fractures, including 7 fractures in the first set, 4 fractures in the second set and 9 fractures in the third set (Figure 7). There are 9 independent fractures which are not connected to the existing hydraulic fracture network and removed from the model. Hence, as shown in Figure 8a, the fracture model for fluid computation contains 11 fractures, of which 5 fractures are the primary fractures intersecting with the grouting borehole and the parameters of each fractures are arranged in

Table 2. Geometric parameters of the computation fracture mode.

Fracture	Coordinate/m			Radius/m	Dip Direction/Deg	Dip Angle/Deg	Aperture/m
	X	Y	Z
1-1724	19.874	9.182	56.403	0.835	344.432	18.027	0.0089
1-1647	19.053	9.902	55.445	1.213	347.699	13.503	0.0129
1-872	19.157	8.701	55.255	0.726	347.079	14.107	0.0077
1-952	17.818	10.036	52.685	0.543	349.080	15.092	0.0058
2-220	20.193	7.555	57.400	1.246	41.687	80.996	0.0132
2-562	18.444	14.515	54.789	0.755	41.635	86.409	0.0080
2-121	17.862	13.714	51.170	1.162	39.079	85.291	0.0124
3-1755	18.708	7.543	56.899	1.150	9.167	14.944	0.0122
3-2007	16.190	8.332	56.259	1.099	13.901	11.962	0.0117
3-2881	19.736	7.530	55.015	0.784	7.736	10.360	0.0083
3-699	19.515	9.665	54.138	1.035	9.766	19.581	0.0110

¹ The primary fractures are marked in bold in the table.

. Figure 8b shows the fracture grid model, and the number of cell meshes is 2,319,281. Figure 8c is the rock mass geometrical model and Figure 8d is its grid model with 34,346 elements. The conditions of the simulated calculation are as follows: the water–cement ratio of cement grout is 1:1, the grouting pressure is 1.0 MPa, and the permeability rate is 7.23 Lu. and the physical mechanical parameters of rock masses is shown in

Table 3. Physical mechanical parameters of rock masses.

Medium	Density (kg/m3)	Elastic Modulus (MPa)	Poisson Ratio μ	Cohesion C (MPa)	Internal Friction Angle φ (°)
Slate	2690	4000	0.25	9.5	37.7

.

3.2.1. Analysis of Fluid Calculation

1. Analysis of Grout Diffusion Process

Driven by grouting pressure, the injected grout penetrates into the primary fractures intersecting with grouting borehole and then migrates into the fracture network. As shown in Figure 9, the grout diffusion length increases gradually over time. After 60 s, most of the primary fractures have been filled with grout except primary fracture 1-1647; this is because the radius of fracture 1-1647 is larger and there are two secondary fractures intersecting with it, which hinders the further diffusion of grout in it. After grouting, primary fractures are completely filled with grout, secondary fractures are mostly filled, and other fractures are filled with a little grout.

For details, five primary fractures (1-1724, 1-1647, 1-872, 3-2881, 1-952) intersecting directly with the grouting borehole are the main passages of the grout, and then the grout diffuses from the primary fracture to the secondary fractures. From Figure 9, it can be seen that there are secondary fractures intersecting with primary fractures, and the grout diffusion at the intersection presents a non-uniform spread shape. Thus, the filling of the primary fractures is influenced by the number of intersecting secondary fractures. As shown in Figure 9, the filling rate of secondary fracture 2-220 is 100%, this is due to the aperture of fracture 2-220 being 13.2 mm, which is bigger than other fractures, and it intersects with two primary fractures, of which the intersecting length with primary fractures 1-1647 is relatively large, reaching 1.6447 m. Secondary fracture 2-562 is basically filled with grout because it intersects with three primary fractures (1-1647, 1-872, 3-2881), but the length of intersection is short, and its aperture is 0.8 mm which is relatively small resulting in incomplete grout filling. Secondary fracture 2-121 is not fully filled with grout due to the hindrance of fracture 3-699 and the small intersections with fracture 3-2881 and 1-952. In summary, the filling of the secondary fractures is related to the fracture aperture and the intersection length between secondary facture and primary fracture. For other fractures (3-1755, 3-699), the grout comes from secondary fractures and they are far from the grouting borehole, so the filling rate is low.

In the actual project, the grout borehole spacing is 1.5 m, and thus the fractures that are not fully filled will be grouted by adjacent grouting boreholes. In addition, the real cement consumption is 244.74 kg, and the simulated cement consumption is 278.41 kg with an error of 13.76% compared to the real cement consumption. Therefore, considering the complexity of geological conditions, there is a good agreement between the simulated and actual values, and this approach can be used effectively to determine the grouting effect.

The maximum radial diffusion length variation of the primary fractures is analyzed and compared in Figure 10. At the beginning of grouting, there is no significant difference in the variation of grout diffusion length with time in each fracture, after 10 s, the difference of diffusion rate in each fracture gradually appeared. The diffusion rate in fracture 1-952 is the lowest because its aperture is the smallest with 5.8 mm and the diffusion rate in fracture 1-1647 is the fastest because its aperture is the biggest with 12.9 mm. As for the other three fractures (1-1724, 1-872, 3-2881), they have similar fracture apertures, but their grout diffusion rate is different, this is due to their different dip angles, and the grout diffusion rate increases with the dip angle. In summary, fracture aperture and dip angle will affect grout diffusion rate and the grout diffusion rate is proportional to them. Besides, fracture 1-1647 with the biggest aperture but smaller dip angle, and the maximum grout diffusion rate indicates that the fracture aperture is the main effect factor.

2. Comparison Analysis

In order to further analyze the simulation results, the effects of fluid–structure interaction in the grouting process are investigated through a comparison with the conventional CFD simulation neglecting these effects.

It can be seen from Figure 10 that the variation trend of the grout diffusion length in each fracture obtained by conventional CFD simulation is same as that obtained by fluid–structure interaction simulation. However, at the same grouting time, the fluid–structure interaction simulation results display a larger grouting diffusion length than the conventional CFD simulation results, that is, the grouting diffusion rate is faster. This can be explained thus: when considering the fluid–structure interaction effect, the fracture deformation occurs due to the grout pressure induces stresses on the surface of the fracture during the grouting process, which will improve the grout penetrability because the grout will diffuse more easily in a fracture with larger aperture. Therefore, ignoring the fluid–structure interaction effects in the simulation of grouting will lead to an underestimation of the grout diffusion ability.

3.2.2. Analysis of Structural Calculation

In this section, the fractured rock mass in the study area is analyzed to explore the effects of fluid–structure interaction on the rock mass. The maximum and minimum principal stresses of rock mass are shown in Figure 11, which indicates that during the grouting process, most of the rock mass is subjected to compressive stress and the maximum compressive stress is 4.936 MPa. The area near the grouting borehole, the local boundary of the rock mass and fractures are partially subjected to tensile stress, and the maximum tensile stress is 3.361 MPa. The maximum tensile stress and the maximum compressive stress of the rock mass during grouting are less than the tensile strength (45.56 MPa) and the compressive strength (6.32 MPa) of the rock mass. Therefore, the rock mass of the grouting area will not generate new fractures or lifting deformation.

In order to understand the influence of grout on the fracture aperture during grouting process, the aperture variation of fracture 1-1724 is analyzed at 50 s and 100 s. As shown in Figure 12, the maximum aperture increment happens at the intersection of grouting holes and fracture, and when at 50 s, the minimum aperture increment occurs at about 0.6 m where the grout front diffuses to. Then at 100 s, the grout diffusion length increases, the minimum aperture increment occurs at about 1.0 m, and the increment of the aperture along the grout diffusion distance is larger than that at 50 s. In summary, the aperture increment decreases non-linearly along the grout diffusion length and as the grouting time increases, the increment of aperture increases and the fracture deformation range expands.

Figure 13 shows the displacement vectorgraph of fractured rock mass, it can be seen that the aperture variation trend of primary fractures is approximately the same, for each fracture, the largest displacement occurs at the intersection of the fracture and grouting borehole, and the smallest at the end of the fracture. For secondary fractures, the displacement happens at the intersection of primary and secondary fractures. As shown in Figure 13, the maximum displacement of the whole area occurs at the lower fracture surface of fracture 1-1647, i.e., the region ②. This is because that the distance between fracture 1-1647 and fracture 1-872 is 0.43 m which is relatively close, fracture 2-220 intersects with fracture 1-1647 and inserts into the rock mass between fracture 1-1647 and 1-872, so the rock mass in this region is the most vulnerable and prone to deformation. Similarly, the region ① is also the case where the fracture 2-220 intersects the fracture 1-1724, so that the rock mass in the region ① becomes less rigid, but the fracture 1-1724 is farther from the surrounding fractures, so the deformation of the region ① is smaller than the region ②. For fracture 1-872, we can see that the displacement of its upper surface is small. This is because the grouting process is conducted from top to bottom, so the fracture 1-1647 is grouted and deformation occurred before fracture 1-872, which offsets most deformation of the upper surface of the fracture 1-872. As for region ③, the distance between the fracture 3-2881 and the fracture 1-872 is small with 0.293 m, but the rock mass is still relatively rigid due to only a small portion of fracture 2-562 being inserted in the rock mass of this region, so that region 3 does not have large displacement. Overall, the maximum displacement of the entire area is 0.1181 mm, which is very small, so there will be no generation of new fractures and other adverse conditions in this area under the action of grouting.

In summary, the effects of fluid–structure interaction between the grout and the rock mass will affect the grout diffusion process and the rock mass deformation. Therefore, grouting process simulation considering the fluid–structure interaction can better reproduce the grout diffusion and rock deformation process and explore the grouting mechanism under real conditions.

3.3. Parameter Analysis

The parameters affecting the grouting process include grouting pressure, water–cement ratio and elastic modulus of rock mass. In order to discuss the influence of different parameters on the grouting process, the diffusion length of fracture 1-1724 at 20 s is analyzed at different grouting pressures (0.5 MPa, 1.0 MPa, 1.5 MPa) and different water-cement ratio grouts (0.7, 1.0, 2.0), and the maximum displacement of the rock mass is analyzed at different grouting pressure (0.5 MPa, 1.0 MPa, 1.5 MPa) and different elastic modulus (4 GPa, 10 GPa, 20 GPa, 40 GPa). The results are shown in Figure 14 and Figure 15, respectively.

As shown in Figure 14, under the same grouting pressure, the grout diffusion length increases with the water–cement ratio due to the larger the water-cement ratio, the smaller the density, and the yield stress and plastic viscosity of the grout increasing with the decrease of density, which results in a faster grout flow rate. In addition, in the case of the same grout water–cement ratio, the grout diffusion length increases as the grouting pressure increases. Furthermore, as shown in Figure 15, under the same grouting pressure, with the increase in elastic module, the maximum displacement of the rock mass decreases, and with the increase in grouting pressure, the maximum displacement of the rock mass increases accordingly. The maximum displacement under 1.5 MPa is small with 0.1391 mm, indicating that the grouting pressure can be appropriately increased to 1.5 MPa to improve the grouting efficiency without causing adverse conditions.

4. Conclusions

The purpose of this study is to investigate the grouting process in fractured rock mass and the influence of grout on the fracture network; thus, a grouting simulation approach considering fluid–structure interaction based on the 3D fractured network model is developed, and the hydropower station X is taken as a case study to do some in-depth analysis using the proposed approach. The results show that the grouting simulation considering fluid–structure interaction is more realistic and can simultaneously reveal the grout diffusion and fracture deformation, which indicated that it is necessary to consider the effect of fluid–structure interaction in grouting simulation. The main conclusions of this study are as follows:

(1)

In fracture network modeling studies, fracture aperture values are often ignored or inaccurate, which will affect the authenticity of grout simulation. Combined with the exposed surface fracture catalog data, we use the relationship between fracture apertures and trace lengths to obtain a more realistic value of fracture aperture and to establish a more reliable model for numerical simulation of grouting.

(2)

During the grouting process, the filling of the primary fractures is influenced by the number of intersecting secondary fractures, whilst the filling of the secondary fractures is related to the fracture aperture, and the length of the intersection between the secondary facture and the primary fracture. Fracture aperture and dip angle have a significant effect on the grout diffusion rate, while the fracture aperture is the major influencing factor. Moreover, the effect of fluid–structure interaction between the grout flow and the rock mass has a certain influence on the grout diffusion length and neglecting this effect will cause an underestimation of the grouting performance.

(3)

When the fractures in a certain region intersect with each other and are close to other fractures in the surrounding area, the rock mass between such fractures will be the least rigid and prone to deformation during the grouting process.

(4)

Grouting pressure, grout water–cement ratio and rock mass elastic module all have effects on the grouting process. Therefore, in the grouting construction process, the appropriate grouting pressure and grout water–cement ratio should be selected according to different geological conditions.

(5)

The effects of fluid–structure interaction between the grout and the rock mass will affect the grout diffusion process and the rock mass deformation. Therefore, the grouting process simulation considering the fluid–structure interaction can better analyze grout diffusion and rock deformation, and hence explore the grouting mechanism under real conditions.

These results can provide an important theoretical basis and valuable information for grouting construction, and in future research the grouting performance considering the effect of underground water can be further studied.