Research Article  Open Access
A Quasi3D Numerical Model for Grout Injection in a Parallel Fracture Based on Finite Volume Method
Abstract
The purpose of the present work is to develop a quasi3D numerical method that can be used to study the diffusion mechanism of grout injection in a rock fracture based on the collocated structured grid of the finite volume method (FVM). Considering the characteristics of fracture in geometry that the aperture is much less than its length and width, the HeleShaw model is introduced to deduce the zderivatives of velocities u and v at walls, which is a function of the relevant average velocity and the fracture aperture. The traditional difference scheme for the diffusive term is partly substituted with the derived analytical expressions; hence a threedimensional problem of grout flow in the parallel fracture can be transformed into a twodimensional one that concerns fracture aperture. The new model is validated by the analytical solution and experimental data on three cases of grouting in the parallelplate fracture. Compared with the results from ANSYSFluent software, the present model shows better agreement with the analytical solution for the distribution of pressure and velocity. Furthermore, the new model needs less grid unit, spends less time, but achieves greater accuracy. The complexity of the grout flow field in the rock fracture is reduced; thus the computational efficiency can be improved significantly.
1. Introduction
Grouting is widely used as the main means of antiseepage reinforcement of fractured rock mass in the fields of civil engineering, water conservancy and mining. The knowledge of grout diffusion mechanism in the rock fracture can significantly guide the design and construction of the grouting project. The rock mass is usually considered as a discontinuity separated by fractures, in which the grout flows, spreads, and finally seals the channel. For decades, many researchers suggested that the grout propagation in a fracture may have a close relationship with the fracture aperture, the grout properties, and the grouting procedure such as the grouting pressure and time [1]. The theoretical analysis, experimental investigation, and numerical simulation are all effective ways to carry on the correlational study.
In terms of the analytical method, the first discussion of Bingham fluid flow in the fracture channel emerged with the work of Wallner [2]. The grout transport was described as linear flow in a channel, and the penetration length was calculated by the onedimensional channel flow model. Dai and Bird [3] proceeded to study the radial flow rate of a Bingham fluid between the parallel plates. Taking the rheological properties of the grouts into account, Hässler [4] used similar flow rate equations to estimate the flow of cementitious grout in fractures. Subjected to pressure gradients, the transient flow of Bingham viscoplastic grouts in parallel fractures was solved analytically by Amadei and Savage [5]. The results show that the extent of a rigid central layer in fractures depends on the pressure gradient, the grout properties, the fracture aperture, and surface roughness. After that, the analytic relations of grout spread to time period were presented systematically by El Tani [6], including channel flow between parallel plates and radial flow between parallel discs, nonlinear Newtonian fluids, and Bingham material. Considering the influence of inertia forces in the governing equation, a simplified model for predicting grout flow was proposed by Xiao et al. [7]. Zhang et al. [8] derived the analytical expressions of evaluating the propagation of cement and sodium silicate grout in a planar fracture, revealing that the rheology of the grout has a most significant influence on the penetration length. All the abovementioned models share a similar basic assumption that the grout propagation in the fracture is onedimensional or twodimensional flow.
Meanwhile, some field tests are conducted to verify the analytical model [9–11]. However, the concealment of rock fractures causes great difficulty for site test, while lab experiment is more feasible to observe the grout flow in the fracture. In this regard, the parallelplates fracture, formed by two stiff rectangular plates or discs, is the most widely used model to investigate the grout flow and propagation [12–15]. On this basis, the effects of the grout properties, the grouting procedure, the geometry of the fracture, and the underground environmental conditions such as underground water on the diffusion behavior of grout in the fracture were studied [16–18].
Compared to the analytical or experimental method, the numerical simulation has great advantages in predicting the grout propagation with wide applicability, more efficiency, and low cost. It is flexible to simulate the grouting under different conditions and analyze the influences of various factors on the grouting effect. Based on the work of Wallner [2], Hässler et al. [19, 20] made the first effort to simulate the grout flow. The fracture plane in a rock mass was described as an orthogonal network of onedimensional channels. The study was improved by Eriksson et al. [21]. In that work, a fracture network with varying aperture, constructed from the onedimensional rectangular channel elements, was modeled to calculate the grout spread and predict the grouting result. Based on MonteCarlo method, the stochastic fracture networks were generated to predict the grout penetration [22, 23]. Similarly, a stochastic discrete feature network model incorporating boreholes and the fractures intersecting them was presented to evaluate the foundation grouting performance in fractured rock [24]. However, this approach is limited by the requirement of sitespecific fracture information. Mohajerani et al. [25, 26] developed an efficient algorithm of explicit grout forehead pressure to simulate the grout propagation within the fractures. In brief, whether applying a structured or stochastic network, the aforementioned works mostly take onedimensional channel model as the basic unit. This base model always assumes that the grout flows in a straight line and occupies the entire crosssection along the channel. As a matter of fact, at the initial stage of grouting into a largesize fracture, the grout diffuses in two even three dimensions before the flow front of the grout reaches the fracture boundary. Only when the grout is constrained by the boundary and develops fully, it can be approximated as onedimensional flow.
In recent years, the quicksetting grout, such as cement and sodium silicate grout, has been widely applied to prevent the seepage in fractured rock mass [27–29]. The curing time of such grouts is short, generally within a few minutes even in tens of seconds in the water disaster prevention engineering. Therefore, it is very difficult for the grout to reach a onedimensional flow state during a short period of time before curing in the largesize fracture. In this case, it may be more appropriate to conduct twodimensional or threedimensional numerical analysis for the grout propagation process based on the computational fluid dynamics theory.
The aim of this paper is to provide an available numerical method for the study of such grouting engineering. Based on the conventional FVM, a quasi3D numerical model of grout injection in the rectangular parallel fracture is presented considering the fracture aperture. Because the fracture is length and width far greater than the aperture, the HeleShaw approach is employed to deduce the partial derivative of velocity with respect to z, which is related to the average velocity and the fracture aperture. In this way, the threedimensional flow of grout in the fracture is converted to a twodimensional one. The proposed model is solved by the SemiImplicit Method for PressureLinked Equations (SIMPLE) algorithm and the moving interface of the grout tracked by the volume of fluid (VOF) method. The new method solves faster with higher precision than ANSYSFluent. The numerical results, including the diffusion radius, the pressure and velocity distribution, are compared with the analytical solution and experimental data to verify the validity of the new method.
2. Quasi3D Numerical Model
2.1. Governing Equations
For analyzing the grout flow in a parallel fracture illustrated in Figure 1, the following assumptions are made:(i)The grout is an incompressible Newtonian fluid.(ii)Flow is laminar.(iii)Fracture walls are stiff.(iv)No slip occurs at the solid interface.
It is also assumed that the velocity in the zdirection is zero; accordingly, there is no pressure gradient. It can be expressed as
The geometric feature of the rock fracture is that the aperture value is much smaller than the length and width; therefore the HeleShaw model [30] can be used for analyzing the grout flow in the fracture, with the following flow equations:
Based on Stokes’ hypothesis, the Newtonian flow of the grout in the fracture can be expressed by the mass and momentum conservation equations. The mass conservation equation is given bywhere is the density and u is the velocity vector.
The momentum equations can be written as follows:where p is the pressure and is the viscosity. u, v, and are the flow velocities in the x, y, and z directions, respectively. Similarly, , , and are three components of the constant source term.
If introducing a general variable ϕ, (4) and (5) can be expressed as the following form:where is the diffusion coefficient and ϕ stands for velocity u, v, or in momentum equations. Equation (6) is called the transport equation for variable ϕ. It clearly highlights the various transport processes: the rate of change term and the convective term on the lefthand side and the diffusive term and the source term, respectively, on the righthand side.
2.2. ThreeDimensional Control Volume Element
The simulation domain is divided into a number of discrete control volumes (CVs) based on a Cartesian grid of equally spaced horizontal and vertical grid lines. For each control volume, the continuity and momentum equations are approximated using algebraic expressions involving the values of the unknown u, v, p at the center of that CV and at the centers of neighbouring CVs [31]. A typical threedimensional control volume element is shown in Figure 2. Positioned midway between adjacent nodes, the six faces of the control volume are labelled n, , n, s, t, and b, which represent East, West, North, South, Top, and Bottom. The centroid of the element is identified by the nodal point P and its neighbours in a threedimensional geometry, the nodes to the west, east, south, north, top, and bottom, are denoted by W, E, S, N, T, and B, respectively. The positive directions along the coordinate are also given.
2.3. Analytical Expressions for Derivatives of Velocity
The two sides of (2) and (3) are integrated twice with respect to z, which derives the formulation as follows:with the boundary condition of u = v = 0 at the fracture wall. The zderivatives of flow velocities at the boundary are expressed as
According to the cubic law for Newtonian flow in the fracture [32], the averaged velocities in the x and ydirections are separately given as
Respectively substituting (11) and (12) into (9) and (10) yields the following equations:
If the general variable is introduced, (13)(14) can be usefully written in the following form:
2.4. Discretization of the Transport Equations
The key step of FVM is the integration of governing equations over a control volume and over a time interval to yield discretized equations at its nodal point P. For the threedimensional control volume defined in Figure 2, from t to t+Δt this gives
Taking the average velocities and as unknown, the volume integral of the rate of change term can be written as
The resulting volume integral in the convective term is converted into an integral over the entire bounding surfaces of the control volume by using Gauss’s divergence theorem. Due to the assumed condition , the convective term can be transformed into the following equations:Similarly, applying Gauss’s divergence theorem and the central differencing method as well as substituting (15)(16) into the diffusive term, the volume integral in the diffusive term can be written as follows:The source term is approximated by means of a linear form:The fully implicit scheme is applied to conduct the time integral of Eq. (17). Then two variables F and D are defined to represent the convective mass flux per unit area and diffusion conductance at cell faces, which are calculated with the following formulae:
By substituting (18)(22) into (17), the threedimensional numerical calculation is transformed into twodimensional calculation, and the corresponding discretized momentum equation can be rearranged aswherewithThe neighbour coefficients of (23) for the exponential scheme are shown as follows:
The discretized form of continuity equation for the control volume gives
Based on SIMPLE algorithm, the discretized continuity equation will be turned into a pressure correction equation. A variable, pressure correction , is firstly defined as the difference between the correct pressure and the guessed pressure. The pressure correction values of nodal point P and its neighbours are denoted as , , , , and , the subscript of which indicates the position of the node in the threedimensional geometry.
The pressure correction equation is treated in a similar manner. Its discretized equation is written aswherewith the coefficients as follows:
From the above calculation process, it can be seen that the threedimensional flow has been transformed into a twodimensional problem by introducing the hypothesis of HeleShaw model and deducing the analytical expression for partial derivatives of velocity with respect to z. The present model taking the fracture aperture into account only needs to make the twodimensional numerical calculation; therefore it is called the quasi3D model.
2.5. VOF Model
The grout diffusion in the fracture is a dynamic process. The flow problem involves two immiscible phases separated by a sharp interface. One phase is the grout, and the other is usually either air or water. As the grout propagates, the interface between grout and air in a fracture gradually changes. The grout flow front can be positioned by the interface tracking technology such as MAC method, LevelSet method, and VOF method [33]. In this study, the VOF method is selected to calculate the multiphase flow. A fractional volume function, which is defined as the fractional volume of one fluid in each cell, is used to reconstruct and track the interface at each moment. The evolution of the grout position and shape depends on how the function changes with time.
The volume of one phase in each cell is calculated firstly by building a scalar function as follows:where and represent the area occupied by air and grout, respectively. The scalar function satisfies the following equation:where t is time. u and v are the velocities in x and ydirection, respectively.
The computational domain is divided into multiple noncoincident cells. The volume of a cell is denoted as . The integral of λ(x,y,t) over each cell I_{i,j} is F_{i,j}, which can be written asThe function F_{i,j} is called the fractional volume function, which obeys the following convection equation:
Equation (34) is the VOF equation and its solving results can be divided into three cases: if F = 1, only the grout is present in the cell; if F = 0, the cell is full of air; if , the phase interface can be in the cell. According to the value of F, the corresponding interface can be constructed.
It is obvious that the key to tracking the moving interface is solving the VOF equation [34]. There are many methods at present, such as the HirtNichols method, the Flair method, the Youngs method, and the FCT method. The Youngs method based on geometric principles has higher computational accuracy compared with other methods [33, 35]. The basic implementation steps of the Youngs method include (1) constructing the fluid interface according to the value of F in the cell; (2) advancing the moving interface (updating the VOF function) according to the current flow velocity. Repeatedly execute the above two steps in each time step to track the movement of the interface between grout and air.
2.6. Numerical Details
The overall numerical procedure for predicting the grout diffusion in the fracture is as follows:(1)Fractional volume function F, velocity vector u, and pressure p are initialized at all cells and the boundary conditions are defined.(2)Material properties at all cells are initialized according to the time t and fractional volume function F.(3)Momentum equations are solved to obtain the current velocity.(4)Pressure correction equation is solved to obtain pressure correction value.(5)Update current velocity and pressure.(6)Repeat steps (3)(5) until the flow field converges.(7)Update fractional volume function using Youngs method.(8)Return to step (2). Time step is advanced and the numerical process proceeds to the next iteration.
3. Results and Discussion
3.1. Validation of Numerical Model
3.1.1. Case 1: h = 2 mm
In order to validate the proposed model, an example of grouting a parallel fracture with a constant rate was taken. Its analytical solution is listed in Appendix A.
The fracture is 1 m long and 1 m wide with an aperture of 2 mm, as shown in Figure 3. A noslip boundary condition is used for upper and lower wall surfaces and the pressure outlet boundary for four sides. A grouting hole, 10 mm in radius, is located at the center of upper wall. It is assumed that the fracture is initially filled with air under ambient pressure of zero. The grout with a density of 1500 kg/m^{3} and viscosity of 0.04 Pa·s is pumped through the grouting hole at a constant flow rate of 37.7 mL/s. Air viscosity at normal temperature is 1.8×10^{−5} Pa·s and its density is 1.205kg/m^{3}. Considering the geometric symmetry of the fracture model, the simulation was simplified by modeling only a quarter of the flow field, as shown in Figure 4. Case 1 was computed separately by the present quasi3D numerical model and the commercial code ANSYSFluent. In the new model the element length in the zdirection equaled the aperture, while the fracture was decomposed into 2 layers and 6 layers in the direction perpendicular to the wall before the calculation is performed using Fluent. All the calculations were finished by the same computer with the following basic configuration:(i)CPU: Intel Core i7 7700(ii)Mainboard: ASUS B150PLUS(iii)Memory: 16 GB (DDR4 2401MHz), 400
In the quasi3D numerical model, the computational domain was meshed into 2500 structured hexahedral elements at the interval of 10 mm along both length and width. The edge length in the zdirection was 2 mm, and the time step was set to 0.1 s. The calculation kept going until the front flow of grout reached the boundary. The total computing time was about 0.5 hours.
Simultaneously, the grout diffusion behavior in the fracture was simulated by ANSYSFluent 17.0 under the same boundary condition. The simulated domain was divided into 500064 hexahedral cells with an edge length of 1 mm in the zdirection. Most of the meshes had an edge length between 1 mm and 3.5 mm. A variable time step size limited by global Courant number was adopted, the maximum time step was less than 0.006 s, and the calculation finished after 48 hours. Correspondingly, when the fracture aperture was divided into 6 layers, the total number of cells was 1500264. The minimum edge length is 0.3 mm, the time step was less than 0.003 s, and the total computing time was more than 432 hours.
Figures 5–7 display the advancement of the grout flow front as a function of time in the fracture. The diffusion front positions of grout at t = 10 s, 20 s, 30 s, and 40 s are represented based on three different models. The results show that the filling range of grout gradually increases in the fracture with the continuously grouting. Besides, the diffusion shape is always circular, which is consistent with the diffusion pattern described by the analytical model.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The numerical results and analytical solution of the grout diffusion radius are shown in Table 1. According to these data, the variation curves of grout diffusion radius with time are plotted in Figure 8. The average relative errors of the Quasi3D model and 2 layers and 6 layers model are 0.41%, 0.37%, and 0.38%, respectively. The numerical predictions of the three models are all in good agreements with the analytical solution for diffusion radius variation with time.

The grout pressure at different positions calculated by the three models is presented in Table 2 at t = 20 s and Table 3 at t = 40 s. The corresponding average flow velocity is summarized in Tables 4 and 5. The contrastive pressure distribution over radius is plotted in Figure 9. The comparative study for the distribution of average flow velocity at different times is illustrated by Figure 10. It shows that the presented model agrees with the analytical solution better than other models in the distributions of the pressure and velocity. For the new model, the maximum relative errors of the pressure and velocity are 4.44% and 1.91% at t = 20 s, respectively, while 4.25% and 1.90% at t = 40 s, respectively. By contrast, 6 layers model calculated by ANSYSFluent has a lower accuracy. At t = 20 s, the relative errors of the pressure and velocity are up to 7.30% and 5.63% separately. At t = 40 s, they are 7.52% and 5.77%, respectively. Then the solution accuracy of 2 layers model calculated by ANSYSFluent is further reduced. The maximum relative error is more than 33%.




(a)
(b)
(a)
(b)
In general, when grouting at a constant rate, both the grout pressure and the velocity decrease with the increase of the radius, and the rate of change declines gently. Their maximum values are located in the grouting hole. Obviously, as the diffusion radius increases, the grouting pressure rises up. The trend is in accord with the analytical model.
Table 6 shows an overview of the quasi3D model and 2 layers model and 6 layers model about the implementation details and the average relative errors of solving results for the diffusion radius, the pressure distribution, and the average velocity. It is obvious that the quasi3D model takes less time but obtains a higher solving precision. The new model has more advantage in simulating the grout diffusion in the fracture than others.
 
Note: △x is the edge length of a grid cell. 
3.1.2. Case 2: h = 0.1 mm
The second verification example is to calculate the diffusion characteristics of grout in a parallelplate fracture under constant grouting pressure, and its analytical solution is listed in Appendix B.
In Case 2, the grout was pumped into a fracture of size 1 m × 1 m × 0.1 mm with a constant grouting pressure of 200 kPa. As Case 1 has mentioned, the fracture was still initially full of air under ambient pressure of zero. The density and viscosity of the grout are 1500kg/m^{3} and 0.02 Pa·s, respectively. Other conditions were identical with Case 1. Still only a quarter of the fluid region was calculated. For the quasi3D model, the number of the grid unit is 2500.
Figure 11 shows the change of grout diffusion front position with time under the grouting pressure of 200 kPa. The interface between the grout and air gradually goes forward away from the grouting hole in a circular way as the injection keeps going, which is the same as the diffusion behavior obtained by the analytical model.
(a)
(b)
(c)
(d)
(e)
Table 7 presents the comparison of the pumping rate and grout diffusion radius between the simulation results and the analytical solution, and their contrast curves are plotted in Figures 12 and 13, respectively. It can be seen that the diffusion radius gradually increases over time while the diffusion rate descends slowly under constant grouting pressure. The results of the quasi3D model agree well with the analytical solution. The relative error of the diffusion radius is less than 1.25%, and of the pumping rate is 1.6%.

Tables 8 and 9 show the contrast of the radial pressure and average velocity calculated by the quasi3D model and the analytical model at t = 20 s and t = 50 s, respectively. The distributions of the pressure and average velocity at different moments are plotted in Figures 14 and 15, respectively. It demonstrates that the developed numerical model matches well with the analytical result. The maximum average relative error of the pressure is about 1.91%, and of the average velocity is less than 1.95%. Figure 15 indicates that the average flow velocity in sections declines over the diffusion radius, which is due to the decrease in the pressure gradient of the flow field and pumping rate. To summarize, at different times the results of the pressure and average velocity obtained by the numerical model are in good agreement with the theoretical solution.


(a)
(b)
(a)
(b)
3.2. Comparison with Experimental Results
To further verify the applicability of the quasi3D model, an experiment of grout injection in a fracture was carried out. The main device consists of a transparent 12 mm thick plexiglass plate, a 20 mm thick steel plate and 1mm thick rubber gasket in the middle of both to form the aperture, as shown in Figure 16. The grout with a viscosity of 0.2 Pa·s was injected through a 1 cm diameter grouting pipe mounted on the bottom of steel plate at a pumping rate 5 mL/s. The grout diffusion process was taken photos by a digital camera.
(a)
(b)
The quasi3D model, set the same condition as experimental, was used to simulate the grout diffusion process. Due to the symmetric structure of the fracture, a quarter of simulated domain in the upper right was selected to implement the numerical work, as shown in Figure 17. The domain was divided into 50 cells in the xdirection and 30 cells in the ydirection.
Figure 18 shows the grout diffusion region and morphology at t = 5 s, 10 s, 15 s, 20 s, 25 s, and 30 s in the experiment. The grout diffuses in a circular shape in the fracture, which is consistent with the theoretical model. The grout front advances forward, then the filling area becomes larger over the time. Figure 19 displays the comparison of the grout diffusion region and morphology between numerical and experimental results. The solid line represents the experimentally measured front position of the grout. It indicates that the results of both methods are basically in agreement. The grout diffusion radius was obtained by averaging the measured values in four orthogonal directions. Table 10 shows the grout diffusion radius at different times obtained from the analytical, experimental, and numerical methods. The average relative error of the numerical results versus the analytical solution is 0.64%. The average relative error between the experimental results and the analytical solution is 0.86%. To summarize, the results of the grout diffusion distance at different times calculated by the quasi3D model are in good accordance with the experimental and analytical results, demonstrating that the presented model can accurately simulate the grouting diffusion process

(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
4. Conclusion
Numerical simulation is always an important way to study the diffusion mechanism of grout in the rock fracture. Based on the traditional threedimensional FVM, the aim of this paper is to develop a quasi3D model to simulate the grout injection in the parallel fracture considering its geometric characteristics. The present model is validated by the analytical model and the experimental measurement. The main work and conclusion are summarized as follows:(i)The geometric feature of fracture is that the length and width are much larger than the aperture. In the presented model, the average velocity parallel to the fracture surface was used as the unknown variable to discretize the governing equations. The partial derivatives of the velocity with respect to z at fracture walls were deduced by introducing the HeleShaw model. The derived equations were proven to be a function of the aperture and the average velocity. By substituting the expressions into the discretized momentum equations, the threedimensional diffusion of the grout in the parallel fracture can be transformed into a twodimensional problem. Combined with the VOF method of tracking the grout front, a quasi3D simulation method was established to predict the flow and diffusion behavior of grout in the fracture.(ii)The quasi3D model was validated by two cases under different conditions. It was also compared with the experimental results. The calculation results, including the diffusion radius, shape, pressure, and average velocity, show a good agreement with the analytical and experimental results. The contrast of the present model with the simulation by ANSYSFluent indicates that the new model needs less grid unit and spends much less time when they both satisfy the same accuracy requirement.(iii)The proposed model provides an efficient way to simulate the grout injection in the fracture. The process of grout diffusion under different environmental conditions and the influences of different factors on the grout, especially the quicksetting grout, will be focused on in the future study.(iv)This method succeeds to simulate the grout diffusion behavior in a fracture. On this basis, the quasi3D model may be further extended to a threedimensional intersected fracture network in the next phase; thus it is possible to simulate the grout injection in a more complex environment [36, 37].
Appendix
A. Constant Flow
The analytical expressions for grout pressure distribution, diffusion radius, and average flow velocity in different sections are given as follows, respectively [38]:where t is time and R is the diffusion radius at time t. is the kinematic viscosity of grout. q is the pumping rate. h is the fracture aperture. p is grout pressure at radius r, p_{0} is the ambient pressure, and u is the radial average velocity.
B. Constant Pressure
When the grouting pressure is constant, the analytical solutions of the grout diffusion radius, pumping rate, pressure distribution, and the average velocity are shown below, respectively [38]:where t is time and R is the diffusion radius at time t. is the kinematic viscosity of grout. q is the pumping rate. h is the fracture aperture. p is grout pressure at radius r, P_{0} is the grouting pressure, p_{0} is the ambient pressure, and u is the radial average velocity. r_{0} is the radius of the grouting hole.
Symbols
:  Face area 
:  Coefficient of general variable 
:  Constant term in the momentum discretized equation 
:  Source term in the pressure correction equation 
:  Diffusion conductance at cell faces 
:  Coefficient of momentum interpolation equation 
:  Convective mass flux per unit area 
:  Fractional volume function 
:  Fracture aperture 
:  Grouting pressure 
:  Grout pressure 
:  Ambient pressure 
:  Pressure correction value 
:  Pumping rate 
:  Diffusion radius 
:  Radius 
:  Radius of the grouting hole 
:  Source term 
:  Constant term in the linearised form of the source term 
:  Source term of momentum 
:  Source term of momentum 
:  Source term of momentum 
:  Slope of the function of source term with respect to at node P 
:  Time 
u:  Velocity vector 
:  Velocity in the direction 
:  Average velocity in the direction 
:  Velocity in the direction 
:  Average velocity in the direction 
:  Velocity in the direction 
:  Coordinate value 
:  Diffusion coefficient 
:  Time interval 
:  Volume of a grid cell 
:  Volume of a cell 
:  Edge length of a grid cell 
:  Node spacing in the direction 
:  Node spacing in the direction 
:  Node spacing in the direction 
:  Scalar function 
:  Viscosity 
:  Density 
:  General variable in the traditional model 
:  General variable at node P at new time level 
:  General variable at node P at (old) time level 
:  Area occupied by air 
:  Area occupied by air. 
:  Node to the bottom 
:  Bottom face 
:  Node to the east 
:  East face 
:  Node to the north 
:  North face 
:  Nodal point P 
:  Between node P and node E 
:  Between node P and node N 
:  Node to the south 
:  Between node S and node P 
:  South face 
:  Node to the top 
:  Top face 
:  Node to the west 
:  Between node W and node P 
:  West face. 
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest to this work.
Acknowledgments
This research was supported by National Key Research and Development Plan [Grant no. 2016YFC0802207], the Program for Science and Technology Innovation Talents in Universities of Henan Province [Grant no. 19HASTIT041], the National Natural Science Foundation of China [Grants nos. 51878624 and 51878622], Key Research Projects of Higher Education in Henan Province [Grant no. 18A580001], the Natural Science Foundation of Henan province [Grant no. 182300410116], the Program for Innovative Research Team (in Science and Technology) at the University of Henan Province [Grant no. 18IRTSTHN007], and Major Scientific, Technological Special Project in Henan [Grant no. 181100310400], and the Postdoctoral Research Sponsorship in Henan Province [Grant no. 19030026].
References
 G. Gustafson and H. Stille, “Stop criteria for cement grouting,” Zeitschrift für Geomechanik und Ingenieurgeologie im Bauwesen und Bergbau, vol. 25, no. 3, pp. 62–68, 2005. View at: Google Scholar
 M. Wallner, “Propagation of sedimentation stable cement pastes in jointed rock,” in Rock Mechanics and Waterways Construction, pp. 449–461, University of Aachen, BRD, 1976. View at: Google Scholar
 G. Dai and R. Byron Bird, “Radial flow of a Bingham fluid between two fixed circular disks,” Journal of NonNewtonian Fluid Mechanics, vol. 8, no. 34, pp. 349–355, 1981. View at: Publisher Site  Google Scholar
 L. Hässler, Grouting of Rock–Simulation and Classification, Department of Soil and Rock Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden, 1991.
 B. Amadei and W. Z. Savage, “An analytical solution for transient flow of Bingham viscoplastic materials in rock fractures,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 2, pp. 285–296, 2001. View at: Publisher Site  Google Scholar
 M. El Tani and H. Stille, “Grout spread and injection period of silica solution and cement mix in rock fractures,” Rock Mechanics and Rock Engineering, vol. 50, no. 9, pp. 2365–2380, 2017. View at: Publisher Site  Google Scholar
 F. Xiao, Z. Zhao, and H. Chen, “A simplified model for predicting grout flow in fracture channels,” Tunnelling and Underground Space Technology, vol. 70, pp. 11–18, 2017. View at: Publisher Site  Google Scholar
 W. Zhang, S. Li, J. Wei et al., “Grouting rock fractures with cement and sodium silicate grout,” Carbonates and Evaporites, vol. 33, no. 2, pp. 211–222, 2018. View at: Publisher Site  Google Scholar
 M. Eriksson, “Grouting field experiment at the Äspö Hard Rock Laboratory,” Tunnelling and Underground Space Technology, vol. 17, no. 3, pp. 287–293, 2002. View at: Publisher Site  Google Scholar
 J. Funehag and Å. Fransson, “Sealing narrow fractures with a Newtonian fluid: Model prediction for grouting verified by field study,” Tunnelling and Underground Space Technology, vol. 21, no. 5, pp. 492–498, 2006. View at: Publisher Site  Google Scholar
 J. Funehag and G. Gustafson, “Design of grouting with silica sol in hard rock – New design criteria tested in the field, Part II,” Tunnelling and Underground Space Technology, vol. 23, no. 1, pp. 9–17, 2008. View at: Publisher Site  Google Scholar
 H. Alehossein, “Viscous, cohesive, nonNewtonian, depositing, radial slurry flow,” International Journal of Mineral Processing, vol. 93, no. 1, pp. 11–19, 2009. View at: Publisher Site  Google Scholar
 A. Draganović and H. Stille, “Filtration of cementbased grouts measured using a long slot,” Tunnelling and Underground Space Technology, vol. 43, pp. 101–112, 2014. View at: Publisher Site  Google Scholar
 J. Funehag and J. Thörn, “Radial penetration of cementitious grout – Laboratory verification of grout spread in a fracture model,” Tunnelling and Underground Space Technology, vol. 72, pp. 228–232, 2018. View at: Publisher Site  Google Scholar
 M. H. Mohammed, R. Pusch, and S. Knutsson, “Study of cementgrout penetration into fractures under static and oscillatory conditions,” Tunnelling and Underground Space Technology, vol. 45, pp. 10–19, 2015. View at: Publisher Site  Google Scholar
 W. Sui, J. Liu, W. Hu, J. Qi, and K. Zhan, “Experimental investigation on sealing efficiency of chemical grouting in rock fracture with flowing water,” Tunnelling and Underground Space Technology, vol. 50, pp. 239–249, 2015. View at: Publisher Site  Google Scholar
 H. B. Lee, T. Oh, E. Park, J. Lee, and H. Kim, “Factors affecting waterproof efficiency of grouting in single rock fracture,” Geomechanics and Engineering, vol. 12, no. 5, pp. 771–783, 2017. View at: Publisher Site  Google Scholar
 Z. Gailing, Z. Kaiyu, G. Yue, and W. Wenxue, “Comparative experimental investigation of chemical grouting into a fracture with flowing and static water,” Mining Science and Technology, vol. 21, no. 2, pp. 201–205, 2011. View at: Google Scholar
 L. Hässler, H. Stille, and U. Håkansson, “Simulation of grouting in jointed rock,” in Proceedings of the 6th International Society for Rock Mechanics and Rock Engineering Congress, ISRM 1987, pp. 829–832, Montreal, Canada, September 1987. View at: Google Scholar
 L. Hässler, U. Håkansson, and H. Stille, “Computersimulated flow of grouts in jointed rock,” Tunnelling and Underground Space Technology, vol. 7, no. 4, pp. 441–446, 1992. View at: Publisher Site  Google Scholar
 M. Eriksson, H. Stille, and J. Andersson, “Numerical calculations for prediction of grout spread with account for filtration and varying aperture,” Tunnelling and Underground Space Technology, vol. 15, no. 4, pp. 353–364, 2000. View at: Publisher Site  Google Scholar
 M. J. Yang, Z. Q. Yue, P. K. Lee, B. Su, and L. G. Tham, “Prediction of grout penetration in fractured rocks by numerical simulation,” Canadian Geotechnical Journal, vol. 39, no. 6, pp. 1384–1394, 2002. View at: Publisher Site  Google Scholar
 M. Eriksson, Prediction of Grout Spread and Sealing Effect, Royal Institute of Technology, Stockholm, Sweden, 2002.
 D. A. Shuttle and E. Glynn, “Grout curtain effectiveness in fractured rock by the discrete feature network approach,” in Proceedings of the Third International Conference on Grouting and Ground Treatment, pp. 1405–1416, New Orleans, Louisiana, USA, 2003. View at: Publisher Site  Google Scholar
 S. Mohajerani, A. Baghbanan, R. Bagherpour, and H. Hashemolhosseini, “Grout penetration in fractured rock mass using a new developed explicit algorithm,” International Journal of Rock Mechanics and Mining Sciences, vol. 80, pp. 412–417, 2015. View at: Publisher Site  Google Scholar
 S. Mohajerani, A. Baghbanan, G. Wang, and S. Forouhandeh, “An efficient algorithm for simulating grout propagation in 2D discrete fracture networks,” International Journal of Rock Mechanics and Mining Sciences, vol. 98, pp. 67–77, 2017. View at: Publisher Site  Google Scholar
 M. Zhu, Q. Zhang, X. Zhang, and B. Hui, “Comparative study of soil grouting with cement slurry and cementsodium silicate slurry,” Advances in Materials Science and Engineering, vol. 2018, Article ID 1893195, 10 pages, 2018. View at: Google Scholar
 R. Valentino, E. Romeo, and D. Stevanoni, “An experimental study on the mechanical behaviour of two polyurethane resins used for geotechnical applications,” Mechanics of Materials, vol. 71, pp. 101–113, 2014. View at: Publisher Site  Google Scholar
 G. Molinda, “Reinforcing coal mine roof with polyurethane injection: 4 case studies,” Geotechnical and Geological Engineering, vol. 26, no. 5, pp. 553–566, 2008. View at: Publisher Site  Google Scholar
 T. Mitani and H. Hamada, “Prediction of flow patterns in the polyurethane foaming process by numerical simulation considering foam expansion,” Polymer Engineering & Science, vol. 43, no. 9, pp. 1603–1612, 2003. View at: Publisher Site  Google Scholar
 H. K. Versteeg, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Prentice Hall, New Jersey, NJ, USA, 1988.
 T. Chen, L. Zhang, and D. Zhang, “An FEM/VOF hybrid formulation for fracture grouting modelling,” Computers & Geosciences, vol. 58, pp. 14–27, 2014. View at: Publisher Site  Google Scholar
 M. Rudman, “Volumetracking methods for interfacial flow calculations,” International Journal for Numerical Methods in Fluids, vol. 24, no. 7, pp. 671–691, 1997. View at: Publisher Site  Google Scholar  MathSciNet
 C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,” Journal of Computational Physics, vol. 39, no. 1, pp. 201–225, 1981. View at: Publisher Site  Google Scholar
 Y. Liu, N. Ma, and X. Gu, “3D simulation of large amplitude liquid sloshing using YoungsVOF method,” Harbin Gongcheng Daxue Xuebao/Journal of Harbin Engineering University, vol. 33, no. 9, pp. 1075–1078, 2012. View at: Google Scholar
 S. Mohajerani, D. Huang, G. Wang, S.M. E. Jalali, and S. R. Torabi, “An efficient algorithm for generation of conforming mesh for threedimensional discrete fracture networks,” Engineering Computations (Swansea, Wales), vol. 35, no. 8, pp. 2860–2882, 2018. View at: Publisher Site  Google Scholar
 S. Mohajerani, G. Wang, D. Huang, S. M. Jalali, S. R. Torabi, and F. Jin, “An efficient computational model for simulating stressdependent flow in threedimensional discrete fracture networks,” KSCE Journal of Civil Engineering, vol. 23, no. 3, pp. 1384–1394, 2019. View at: Publisher Site  Google Scholar
 Z. Feng, Material Development and Research of Osmosis and Diffusion on Chemical Grouting for Extraordinary Cracked Coal and Rockmass, China Coal Research Institute, Beijing, China, 2007.
Copyright
Copyright © 2019 Xiaolong Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.