Transient Analysis of Grout Penetration With Time-Dependent Viscosity Inside 3 D Fractured Rock Mass by Unified Pipe-Network Method

Grouting is widely used for mitigating the seepage of underground water and enhancing the stability of fractured rock mass. After injection, the viscosity of the grout gradually increases until solidification. Conventional multifield analysis models ignoring such effects greatly overestimate the penetration region of the grout and the stability of the grouted rock structures. Based on the 3D unified pipe-network method (UPM), we propose a novel numerical model considering the time-dependent viscosity of the grout, therein being a quasi-implicit approach of high efficiency. The proposed model is verified by comparing with analytical results and a time-wise method. Several large-scale 3D examples of fractured rock mass are considered in the numerical studies, demonstrating the effectiveness and robustness of the proposed method. The influence of the time-dependent viscosity, fracture properties, and grouting operation methods are discussed for the grout penetration process.


Introduction
The existence of discontinuous fractures provides paths for underground water, which reduce the mechanical strength and increase the permeability of the rock mass.In engineering practices, sealing the fractured rock mass is a critical issue in large-scale rock engineering projects such as mining, civil, and hydrogeology engineering for avoiding potential water and mud inrush [1,2].As one of the most effective methods, grouting is a widely used strategy for improving the mechanical properties of rock mass and mitigating groundwater leakages [3][4][5][6][7].One of the biggest challenges in grouting flow in fractured rock masses is computing the penetration regions of the grout.Because of the complexity of the fracture networks, bare experimental and analytical investigations [8][9][10][11][12] are not sufficient for evaluating their grouting qualities [13][14][15][16].The results achieved by these approaches are mainly suitable for rock masses with a single fracture or regular fracture systems, which cannot reproduce the grouting flow in masses with multiple irregular fractures [17].
The 3D unified pipe-network method (3D UPM) [66][67][68][69] is an efficient numerical approach for capturing the energy/mass transport processes in 3D fracture networks.UPM transforms the 3D discontinuous fracture network into a 3D system with intersected artificial 1D pipe segments.From this point of view, UPM is similar to lattice elements approaches (LEM) [70][71][72], both of which simulate complicate 3D processes in an equivalent lower-dimensional system.On the other hand, UPM is different from LEM as

•
UPM uses pipes for capturing transport processes but ignores mechanical behaviors, while LEM uses trusses for capturing the mechanical behaviors; • UPM introduces pipes parallel to the fracture planes to capture the transport processes in fractures.
On the other hand, in LEM the breakages of trusses represent damage/fractures which originally intersect with the fracture planes, but are not parallel to; • Unlike LEM, UPM cannot simulate fracture propagations and nucleations, which is suitable for simulating transport processes in naturally highly fractured rock mass.
Most numerical models assumes Darcy's flow and uses heat equation [56][57][58] for simulating grouting penetration, which is different from models based on computational fluid dynamics and does not trace the grout flow.Hence, only the distributions of pressure are obtained as results.It is impossible to directly account for the time-dependent viscosity.Most methods assume constant viscosity and ignore the solidification processes of the grout after injection.However, the grouting materials for quickly sealing water-bearing fractures in rock masses possess time-dependent viscosities that increase rapidly after injection due to chemical reactions (solidification), such as with silica sol and cement-silicate sodium grout [73][74][75], which have been proven by experimental investigations [76,77].To the best of our knowledge, there is limited research on simulating the time-dependent viscosity of grout and its impact on the grouting process in the framework of DFN models.
In this work, we present a novel and simple strategy for accounting for the time-dependent viscosity in 3D UPM, which is a quasi-implicit method and exploits the advantages of UPM.A spatialand time-dependent parameter is introduced for indicating the elapsed time of the grouting flow to reach a specific position from the injection point.With this strategy, the viscosity of the grout will also be updated at every time step, providing results in good agreement with experimental investigations.The remaining parts of the paper are organized as follows: In Section 2, the control equations of the grout are briefly introduced, with corresponding discretized forms under the UPM framework.Then, the strategy for accounting for the time-dependent viscosity of the grouting is prescribed.In Section 3, our novel strategy is verified by comparing with analytical solutions and experimental results given in [77].In Section 4, the grouting process in the fracture networks is simulated, and the influential factors of the grouting are analyzed.The main influential factors are found to include (i) the material properties of the grout, (ii) the characteristics of the fractures, and (iii) the grout operation method.Finally, the concluding remarks are given in Section 5.

Rheological Models of Grout
In this paper, we assume single-phase flow in our analysis and consider only the grout flow.The grout penetration depth is estimated by the region with grout pressure higher than or equal to the targeted values.The rheological models of grout are divided into two main categories: Newtonian fluid and Bingham fluid.The linear constitutive model of the Newtonian fluid can be expressed as: where τ is the shear stress (Pa), µ is the dynamic viscosity (Pa • s), and γ is the shear rate (1/s).Fluids such as very-fine-grained cement-based grout (silica sol) and bacterial grout (MICP) can be regarded as being of this type [78].
A Bingham fluid is a type of viscoplastic fluid that can only flow at a higher shear stresses.The constitutive model of a Bingham fluid is still linear, and can be expressed as where τ 0 is the initial yield stress (Pa).When the water-cement ratio is smaller than 1, the cement-based grout can be considered as this type [79].
The grout flow in a one-dimensional fracture is assumed to be laminar, steady, and incompressible.Figure 1a is a schematic illustration of a Newtonian grout flow in a single fracture, where w is the hydraulic aperture (m), dx is the length of a small control volume in the x-direction (m), P g is the grouting pressure (Pa), and P w is the water pressure (Pa).The average flow velocity of a Newtonian fluid in a fracture can be described using the following equation [3]: where x is the global coordinate for the fracture (m), p is the grout pressure (Pa), µ(t) is the time-varying viscosity (Pa • s), and k is the fracture intrinsic permeability (m 2 ).The grout flow in fractures obeys the Darcy's law, thus the permeability for fractures can be written as k = w 2 /12. Figure 1b shows the Bingham grout flow between the parallel plates.Due to the effect of the initial yield stress, the mean velocity of the grout forehead is given:

Pg
where δ = −τ 0 /w is the starting pressure gradient.The grout will only flow in the fracture when the pressure gradient is higher than the starting pressure gradient.

Unified Pipe-Network Method (UPM) Discrete Model
In the UPM framework, the energy/mass transport inside the fracture is equivalently considered as energy/mass transport through the artificial pipes along the fracture boundaries.The equivalent hydraulic parameters of the pipes are derived based the unstructured triangular mesh (see Figure 2).Figure 3 shows a triangular fracture with node i, j, k. o is the center point of the triangular fracture and g, f , h are the corresponding three perpendicular feet of o to the edges.The flow in each pipe can be seen as the flow through the line linking the point o and points g, f , h.Points o, g, f , h separate the triangle into three polygons as i f og, jgoh, and kho f .The areas of these polygons provide information of the control volume of the nodes i, j, k as: Based on [66,67], for the pipe ik, the fluid flow rate Q ik (m 3 /s) shown in Figure 3 can be evaluated as where n of is a unit vector normal to the face o f , which can be calculated as: where l ik is the length of pipe ik.
In each fracture triangle element (composed of nodes i, j, and k), the pressure within the triangle mesh can be written as where N m is the linear shape function, expressed as The pressure gradient in the triangular domain is calculated as: where A ijk is the triangle element area.The coefficients b m and c m depend on the coordinates of the three nodes in each triangle element.These coefficients are represented as For the pipe ik, the flow rate of Newtonian fluid (Q N o f ) can be expressed as: where A o f is the cross-sectional area of the face o f regarding a 3D fracture with aperture, A o f = w l o f in this paper (no filling media inside the fracture), with l o f as the length of o f .l ik is the length of pipe ik.
Similarly, the flow rate of Bingham fluid (Q B o f ) in pipe ik can be expressed as: where the superscripts N and B represent the Newtonian fluid and Bingham fluid, respectively.For the fracture pipe ik, the equivalent conductance coefficient K f ik can be expressed as: With this procedure, the whole fracture network can be transformed into a 3D pipe network (see Figure 4) with equivalent mass transport properties.With the assumptions of the UPM, in this paper, the grout flow is simulated with a one-dimensional Dupuit-Forchheimer model [80,81] for saturated flow as: where S is the storage coefficient, h is the hydraulic head (m), and q s are the source terms (m/s).
For each node, the grout flow obeys the law of mass conservation.Therefore, the governing equation in the UPM can be expressed as: where B i = SV i /g, V i is the control volume (m 3 ) of node i (see Equation ( 5)), ρ is the density of the grout (kg/m 3 ), and g is the gravitational acceleration (m/s 2 ).

Considering the Time-Dependent Viscosity in UPM
In the framework of the UPM, the fractures are equivalently modeled as pipe networks.Points in the considered domain are connected with pipes.For each pipe, the velocity of the grout flow is determined by the pressures of the two nodes.When the pressure distribution is known, the elapsed time for the grout to transport from one point to another can be calculated by the length of the pipe and the velocity of the grout.The velocity of the grout depends on the pressure gradient as well as the viscosity, and the viscosity depends on the elapsed time for the grout to transport in the pipes.Hence, we introduce a new parameter at every node as "the elapsed time for the grout to reach one node from the injection point", denoted as "φ".Obviously, φ = 0 at the injection points of the grout.For a pipe connecting nodes i and j, a simple relationship exists: where L i,j is the length of the pipe.v i,j is the velocity of the grout flow, which is a function of pressure and viscosity at the nodes, as v(µ i , µ j , P i , P j ), built based on Equations ( 3) and (4).Furthermore, the viscosity of the grout is a function of φ (i.e., µ = µ(φ)).Li et al. [82] presented a power law based on experimental investigations: where k is the flow consistency index (Pa • s), n is the flow behavior index, and µ 0 is the initial viscosity (Pa • s).In summary, if P i , P j , and φ i are known, φ j is the only unknown in Equation ( 17), which can be solved by a standard Newton method.To solve for φ at time step t, a quasi-implicit strategy is used: the pressure distribution of time step t − 1 is used for implicitly solving µ at time step i.Compared to the analysis with constant viscosity, our method is very efficient, requiring negligible computing efforts.

Model Verification
In this section, a rectangular fracture (10 m × 10 m) in a cube with dimensions 10 m × 10 m × 10 m is used for verification.The grout was injected into a single fracture with an injection hole (see Figure 6).Two cases were considered: (i) comparing our results to the analytical solutions and (ii) comparing our results to the experimental results given in [77], considering time-dependent viscosities.With UPM, a self-developed mesh generator [83] is adopted to generate triangle elements for the fractures.The parameters employed in the simulation are listed in Table 1.

Verifying the Rheological Models
In this subsection, we first verify the rheological models, where the viscosity of the grout is considered to be a constant value µ 0 .The analytical solutions for Newtonian fluid and Bingham fluid have been proposed by Funehag [3] and Gustanafson [84], which are taken for verification.
Funehag [3] calculated the grout penetration for a Newtonian fluid with constant viscosity µ 0 as: where ∆p is the difference between the injection pressure and initial pore pressure, and t is the flow time.
Gustafson [84] proposed a solution for grout penetration distance (I) and penetration time (t): where t D and I D are relative penetration time and relative penetration length, respectively, t 0 is the characteristic grouting time, I max is the maximum grout penetration, and θ is a ratio influenced by t D .However, the ratio is different for one-dimensional flow and two-dimensional flow: As shown in Figure 7, the grout penetration length along the injection hole obtained by the UPM simulation was consistent with the analytical results at each time step, indicating the effectiveness of our rheological models of grout.

Transient Flow Considering the Time-Dependent Viscosity of Grout
Zhang et al. [77] experimentally investigated the grout penetration length of a type of grout with time-dependent viscosity.The experiments were reproduced with our model.The total simulation time was 60 s, and the viscosity of grout ultimately reached 29.4 Pa • s.Because our strategy is quasi-implicit and the pressure distributions at the last time step were used for calculating the distributions of φ and µ in this step.We attempted different time intervals ∆t for sensitivity analysis.As shown in Figure 8, comparisons between our results and the experimental results indicated that the results were not sensitive to the selection of ∆t when ∆t ≤ 0.05 s.Furthermore, we also tested the influences of UPM discretization.As shown in Figure 9, a total grid number larger than 10, 000 assured the reliability of the results.
Figure 10 shows the comparisons of the grouting pressure distribution between our results and the experimental results at different times, with ∆t = 0.05 s and a grid number equal to 15,514.Generally, our results were in agreement with the experimental results, indicating that the proposed model was capable of simulating grout flows with time-dependent viscosities.

Application in Simulating Grouting Process in Rock Mass with Fracture Networks
The UPM was coded with C++, and all the cases in this section were calculated within 2 h on a normal desktop.Parametric studies were performed to determine the influences of various parameters, including (i) the material properties of the grout (viscosity), (ii) the characteristics of the fractures (number, aperture, and roughness), and (iii) the grout operation method (the grout pressure, the number of boreholes, and the arrangement of the boreholes) on grout penetration in fractured rock masses.

The Influence of Viscosity
The commonly used quick-setting slurry has two components: cement and sodium silicate slurries.The viscosity function µ(φ) depends on the ratio of cement slurry to sodium silicate slurry (C:S).According to the experimental results provided in [85], when C:S = 1, the flow consistency index (k) and flow behavior index (n) are 0.003182 Pa • s and 2.23, respectively, while they are 0.0008427 Pa • s and 2.694 when C:S = 2 (see Figure 11).The other parameters were considered to be the same as listed in Table 1. Figure 12 shows the pressure distributions for a constant viscosity (µ = µ 0 ) and time-dependent viscosity with varying C:S.The results indicate that assuming a constant viscosity greatly overestimated the penetration length of the grout, and the grout flow lasted for a long time.In the cases with a time-dependent viscosity, on the other hand, the grout flow almost stopped after 30 s.

The Influence of Fractures and Grout Operation Method
We generated fracture networks to analyze the influences of the fractures on grout penetration with the same cube model size of 10 m × 10 m × 10 m.Two sets of elliptical fractures were generated by a random function with the parameters listed in Table 2.For each fracture, its center coordinate followed a uniform distribution, and the major and minor axes followed a logarithmic distribution.The Fisher coefficient was 22.The configuration of the random fractures and their mesh are illustrated in Figure 13.The red line represents the borehole through the fractured rock.In the model, there were 7088 nodes and 13, 210 triangle elements in total.The grout parameters were the same as those listed in Table 1.The grout was injected through the borehole with an injection pressure of 30 kPa.The pressure distributions are shown in Figure 14, therein indicating that the pressure distributions became stable after 90 s.Here, we define the grout filling rate as the ratio between the area of grout-filled fractures and the total fracture area.Figure 15 shows the influences of a uniform aperture and joint roughness coefficient (JRC) [86] on the grout filling rate, indicating that the grout filling rate increased with increasing joint aperture and decreased with increasing JRC.The influence of JRC on intrinsic permeability was considered by an empirical equation given in [47].Furthermore, we tested different grout operation methods.Figure 16 shows the influences of the grout pressure for a single borehole, and Figure 16 shows the influences of the borehole number at the same injection pressure of 30 kPa, indicating that the grout filling rate was increased by increasing either the injection pressure or the number of boreholes.

Conclusions
In this study, a three-dimensional numerical model based on the UPM is developed for modeling the transient grout penetration in fractured networks.Rheological models of grout, including Newtonian and Bingham fluids, are considered.A novel quasi-implicit method is presented to account for the time-dependent viscosity of grout, therein calculating the elapsed time φ of the grout to reach specific positions from the injection points.By comparison to analytical and experimental results, the new model is validated.Finally, parametric and case studies are presented regarding a complicated 3D fractured rock mass with boreholes.All the results evidence the effectiveness of our model, which estimates the penetration regions and time for the stabilization of grout with a time-dependent viscosity.

Figure 1 .
Figure 1.Schematic illustration of the grout flow in a fracture joint (a) Newtonian grout flow, (b) Bingham grout flow.

Figure 2 .
Figure 2. Transformation of the mass transport inside a fracture into the mass transport through artificial pipes on the fracture boundaries (regarding triangular mesh).

Figure 4 .
Figure 4. Transforming a 3D fracture network into an equivalent pipe network: (a) domain with fractures (b) domain with equivalent pipe networks.

Figure 5 .
Figure 5. Calculation of φ in the unified pipe-network method (UPM) model: (a) calculate φ from injection point, (b) determining the values of φ from one node to another.

Figure 6 .
Figure 6.Three-dimensional model with a fracture for verification.

Figure 7 .
Figure 7.Comparison of the UPM model with the analytical solutions: (a) Newtonian fluid, (b) Bingham fluid.

Figure 8 .Figure 9 .
Figure 8.The influence of the time step on the UPM results.

Figure 10 .
Figure 10.Comparison of the UPM model with results obtained by Zhang [77].

Figure 12 .
Figure 12.Distribution of pressure at different simulation times.

Figure 15 .
Figure 15.Sensitivity analyses with respect to the fracture aperture and joint roughness.(JRC: joint roughness coefficient.)

Figure 16 .
Figure 16.Sensitivity analyses with respect to the grout pressure and the number of boreholes.

Table 1 .
Parameters for verification simulation.

Table 2 .
Parameters adopted for the generation of the fracture models.

Fracture Number Mean Length (m) St Dev Dip Angle (Degree) Dip Direction (Degree)
Note: St dev is standard deviation.