# A Quasi-3D Numerical Model for Grout Injection in a Parallel Fracture Based on Finite Volume Method

Authors: Xiaolong Li; Meimei Hao; Yanhui Zhong; Bei Zhang; Fuming Wang; Lianbang Wang
Journal: Complexity
 (2019)

---

## Abstract

The purpose of the present work is to develop a quasi-3D 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 Hele-Shaw model is introduced to deduce thez-derivatives of velocitiesu andv 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 three-dimensional problem of grout flow in the parallel fracture can be transformed into a two-dimensional one that concerns fracture aperture. The new model is validated by the analytical solution and experimental data on three cases of grouting in the parallel-plate fracture. Compared with the results from ANSYS-Fluent 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.


---

## Body

## 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 one-dimensional 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 above-mentioned models share a similar basic assumption that the grout propagation in the fracture is one-dimensional or two-dimensional 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 parallel-plates 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 one-dimensional channels. The study was improved by Eriksson et al. [21]. In that work, a fracture network with varying aperture, constructed from the one-dimensional rectangular channel elements, was modeled to calculate the grout spread and predict the grouting result. Based on Monte-Carlo 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 site-specific 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 one-dimensional channel model as the basic unit. This base model always assumes that the grout flows in a straight line and occupies the entire cross-section along the channel. As a matter of fact, at the initial stage of grouting into a large-size 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 one-dimensional flow.In recent years, the quick-setting 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 one-dimensional flow state during a short period of time before curing in the large-size fracture. In this case, it may be more appropriate to conduct two-dimensional or three-dimensional 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 quasi-3D 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 Hele-Shaw approach is employed to deduce the partial derivative of velocity with respect toz, which is related to the average velocity and the fracture aperture. In this way, the three-dimensional flow of grout in the fracture is converted to a two-dimensional one. The proposed model is solved by the Semi-Implicit Method for Pressure-Linked 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 ANSYS-Fluent. 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. Quasi-3D Numerical Model

### 2.1. Governing Equations

For analyzing the grout flow in a parallel fracture illustrated in Figure1, 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.Figure 1
Schematic diagram of parallel fracture model.It is also assumed that the velocity in thez-direction is zero; accordingly, there is no pressure gradient. It can be expressed as(1)
∂
p
∂
z
=
0The geometric feature of the rock fracture is that the aperture value is much smaller than the length and width; therefore the Hele-Shaw model [30] can be used for analyzing the grout flow in the fracture, with the following flow equations:(2)
∂
p
∂
x
=
∂
∂
z
μ
∂
u
∂
z
(3)
∂
p
∂
y
=
∂
∂
z
μ
∂
v
∂
zBased 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 by(4)
div
⁡
ρ
u
=
0where ρ is the density andu is the velocity vector.The momentum equations can be written as follows:(5)
∂
ρ
u
∂
t
+
d
i
v
ρ
u
u
=
-
∂
p
∂
x
+
d
i
v
μ
g
r
a
d
u
+
S
M
x
∂
ρ
v
∂
t
+
d
i
v
ρ
v
u
=
-
∂
p
∂
y
+
d
i
v
μ
g
r
a
d
v
+
S
M
y
∂
ρ
w
∂
t
+
d
i
v
ρ
w
u
=
-
∂
p
∂
z
+
d
i
v
μ
g
r
a
d
w
+
S
M
zwherep is the pressure and μ is the viscosity.u,v, and w are the flow velocities in thex-,y-, andz- directions, respectively. Similarly, S
M
x, S
M
y, and S
M
z are three components of the constant source term.If introducing a general variableϕ, (4) and (5) can be expressed as the following form:(6)
∂
ρ
ϕ
∂
t
+
d
i
v
ρ
ϕ
u
=
d
i
v
Γ
g
r
a
d
ϕ
+
Swhere Γ is the diffusion coefficient andϕ stands for velocityu,v, or w 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 left-hand side and the diffusive term and the source term, respectively, on the right-hand side.


### 2.2. Three-Dimensional 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 unknownu,v,p at the center of that CV and at the centers of neighbouring CVs [31]. A typical three-dimensional control volume element is shown in Figure 2. Positioned mid-way between adjacent nodes, the six faces of the control volume are labelledn, w,n,s,t, andb, which represent East, West, North, South, Top, and Bottom. The centroid of the element is identified by the nodal pointP and its neighbours in a three-dimensional 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.Figure 2
A cell in three dimensions and its neighbouring nodes.


### 2.3. Analytical Expressions forz-Derivatives of Velocity

The two sides of (2) and (3) are integrated twice with respect toz, which derives the formulation as follows:(7)
u
=
∂
p
∂
x
1
2
μ
z
2
-
h
2
4
(8)
v
=
∂
p
∂
y
1
2
μ
z
2
-
h
2
4with the boundary condition ofu =v = 0 at the fracture wall. Thez-derivatives of flow velocities at the boundary are expressed as(9)
∂
u
∂
z
z
=
h
/
2
=
-
∂
u
∂
z
z
=
-
h
/
2
=
∂
p
∂
x
h
2
μ
(10)
∂
v
∂
z
z
=
h
/
2
=
-
∂
v
∂
z
z
=
-
h
/
2
=
∂
p
∂
y
h
2
μAccording to the cubic law for Newtonian flow in the fracture [32], the averaged velocities in thex- andy-directions are separately given as(11)
u
-
=
-
∂
p
∂
x
h
2
12
μ
(12)
v
-
=
-
∂
p
∂
y
h
2
12
μRespectively substituting (11) and (12) into (9) and (10) yields the following equations:(13)
∂
u
∂
z
z
=
h
/
2
=
-
∂
u
∂
z
z
=
-
h
/
2
=
-
6
h
u
-
(14)
∂
v
∂
z
z
=
h
/
2
=
-
∂
v
∂
z
z
=
-
h
/
2
=
-
6
h
v
-If the general variableϕ
- is introduced, (13)-(14) can be usefully written in the following form:(15)
∂
ϕ
∂
z
z
=
h
/
2
=
-
6
h
ϕ
-
(16)
∂
ϕ
∂
z
z
=
-
h
/
2
=
6
h
ϕ
-


### 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 pointP. For the three-dimensional control volume defined in Figure 2, fromttot+Δtthis gives(17)
∫
t
t
+
Δ
t
∫
Δ
V
∂
ρ
ϕ
∂
t
d
V
d
t
+
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
ρ
u
ϕ
d
V
d
t
=
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
μ
g
r
a
d
ϕ
d
V
d
t
+
∫
t
t
+
Δ
t
∫
Δ
V
S
d
V
d
tTaking the average velocitiesu
- and v
- as unknown, the volume integral of the rate of change term can be written as(18)
∫
t
t
+
Δ
t
∫
Δ
V
∂
ρ
ϕ
∂
t
d
V
d
t
=
ρ
ϕ
-
P
-
ϕ
-
P
0
Δ
VThe 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 conditionw
t
=
w
b
=
0, the convective term can be transformed into the following equations:(19)
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
ρ
u
ϕ
d
V
d
t
=
∫
t
t
+
Δ
t
ρ
u
-
e
A
e
ϕ
-
e
-
ρ
u
-
w
A
w
ϕ
-
w
+
ρ
v
-
n
A
n
ϕ
-
n
-
ρ
v
-
s
A
s
ϕ
-
s
d
tSimilarly, 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:(20)
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
μ
g
r
a
d
ϕ
d
V
d
t
=
∫
t
t
+
Δ
t
μ
e
A
e
ϕ
-
E
-
ϕ
-
P
δ
x
P
E
-
μ
w
A
w
ϕ
-
P
-
ϕ
-
W
δ
x
W
P
+
μ
n
A
n
ϕ
-
N
-
ϕ
-
P
δ
y
P
N
-
μ
s
A
s
ϕ
-
P
-
ϕ
-
S
δ
y
S
P
+
-
μ
t
A
t
6
ϕ
-
P
h
-
μ
b
A
b
6
ϕ
-
P
h
d
tThe source term is approximated by means of a linear form:(21)
∫
t
t
+
Δ
t
∫
Δ
V
S
d
V
d
t
=
∫
t
t
+
Δ
t
S
c
+
S
P
ϕ
-
P
Δ
V
d
tThe fully implicit scheme is applied to conduct the time integral of Eq. (17). Then two variablesF andDare defined to represent the convective mass flux per unit area and diffusion conductance at cell faces, which are calculated with the following formulae:(22)
F
w
=
ρ
u
-
w
A
w
,
F
e
=
ρ
u
-
e
A
e
F
n
=
ρ
v
-
n
A
n
,
F
s
=
ρ
v
-
s
A
s
D
w
=
μ
w
A
w
δ
x
W
P
,
D
e
=
μ
e
A
e
δ
x
P
E
D
n
=
μ
n
A
n
δ
y
P
N
,
D
s
=
μ
s
A
s
δ
y
S
PBy substituting (18)-(22) into (17), the three-dimensional numerical calculation is transformed into two-dimensional calculation, and the corresponding discretized momentum equation can be rearranged as(23)
a
P
ϕ
-
P
=
a
E
ϕ
-
E
+
a
W
ϕ
-
W
+
a
N
ϕ
-
N
+
a
S
ϕ
-
S
+
bwhere(24)
a
P
=
a
W
+
a
E
+
a
S
+
a
N
+
μ
A
6
h
+
μ
A
6
h
+
F
e
-
F
w
+
F
n
-
F
s
+
a
P
0with(25)
a
P
0
=
ρ
Δ
V
Δ
t
,
b
=
S
Δ
V
+
a
P
0
ϕ
-
P
0The neighbour coefficients of (23) for the exponential scheme are shown as follows:(26)
a
W
=
D
w
max
⁡
0
,
1
-
0.1
F
w
D
w
5
+
max
⁡
0
,
F
w
a
E
=
D
e
max
⁡
0
,
1
-
0.1
F
e
D
e
5
+
max
⁡
0
,
-
F
e
a
S
=
D
s
max
⁡
0
,
1
-
0.1
F
s
D
s
5
+
max
⁡
0
,
F
s
a
B
=
D
b
max
⁡
0
,
1
-
0.1
F
b
D
b
5
+
max
⁡
0
,
F
bThe discretized form of continuity equation for the control volume gives(27)
ρ
u
A
e
-
ρ
u
A
w
+
ρ
v
A
n
-
ρ
v
A
s
=
0Based on SIMPLE algorithm, the discretized continuity equation will be turned into a pressure correction equation. A variable, pressure correctionp
′, is firstly defined as the difference between the correct pressure and the guessed pressure. The pressure correction values of nodal pointP and its neighbours are denoted as p
P
′, p
E
′, p
W
′, p
N
′, and p
S
′, the subscript of which indicates the position of the node in the three-dimensional geometry.The pressure correction equation is treated in a similar manner. Its discretized equation is written as(28)
a
P
′
p
P
′
=
a
E
′
p
E
′
+
a
W
′
p
W
′
+
a
N
′
p
N
′
+
a
S
′
p
S
′
+
b
P
′where(29)
b
P
′
=
ρ
u
-
A
w
-
ρ
u
-
A
e
+
ρ
v
-
A
s
-
ρ
v
-
A
n
,
a
P
′
=
a
E
′
+
a
W
′
+
a
N
′
+
a
S
′with the coefficients as follows:(30)
a
E
′
=
ρ
d
A
e
,
a
W
′
=
ρ
d
A
w
a
N
′
=
ρ
d
A
n
,
a
S
′
=
ρ
d
A
sFrom the above calculation process, it can be seen that the three-dimensional flow has been transformed into a two-dimensional problem by introducing the hypothesis of Hele-Shaw model and deducing the analytical expression for partial derivatives of velocity with respect toz. The present model taking the fracture aperture into account only needs to make the two-dimensional numerical calculation; therefore it is called the quasi-3D 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, Level-Set 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:(31)
λ
x
,
y
,
t
=
0
x
,
y
∈
Ω
1
1
x
,
y
∈
Ω
2where Ω
1 and Ω
2 represent the area occupied by air and grout, respectively. The scalar function satisfies the following equation:(32)
∂
λ
∂
t
+
u
∂
λ
∂
x
+
v
∂
λ
∂
y
=
0wheret is time.u andv are the velocities inx-andy-direction, respectively.The computational domain is divided into multiple noncoincident cells. The volume of a cell is denoted as△
V
i
,
j. The integral ofλ(x,y,t) over each cellI
i,j isF
i,j, which can be written as(33)
F
i
,
j
=
1
Δ
V
i
,
j
∫
I
i
,
j
λ
x
,
y
,
t
d
VThe functionF
i,j is called the fractional volume function, which obeys the following convection equation:(34)
∂
F
∂
t
+
u
∂
F
∂
x
+
v
∂
F
∂
y
=
0Equation (34) is the VOF equation and its solving results can be divided into three cases: ifF= 1, only the grout is present in the cell; ifF= 0, the cell is full of air; if 0
<
F
<
1, the phase interface can be in the cell. According to the value ofF, 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 Hirt-Nichols 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 ofF 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 functionF, velocity vectoru, and pressurep are initialized at all cells and the boundary conditions are defined.
(2)
Material properties at all cells are initialized according to the timet and fractional volume functionF.
(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.


## 2.1. Governing Equations

For analyzing the grout flow in a parallel fracture illustrated in Figure1, 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.Figure 1
Schematic diagram of parallel fracture model.It is also assumed that the velocity in thez-direction is zero; accordingly, there is no pressure gradient. It can be expressed as(1)
∂
p
∂
z
=
0The geometric feature of the rock fracture is that the aperture value is much smaller than the length and width; therefore the Hele-Shaw model [30] can be used for analyzing the grout flow in the fracture, with the following flow equations:(2)
∂
p
∂
x
=
∂
∂
z
μ
∂
u
∂
z
(3)
∂
p
∂
y
=
∂
∂
z
μ
∂
v
∂
zBased 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 by(4)
div
⁡
ρ
u
=
0where ρ is the density andu is the velocity vector.The momentum equations can be written as follows:(5)
∂
ρ
u
∂
t
+
d
i
v
ρ
u
u
=
-
∂
p
∂
x
+
d
i
v
μ
g
r
a
d
u
+
S
M
x
∂
ρ
v
∂
t
+
d
i
v
ρ
v
u
=
-
∂
p
∂
y
+
d
i
v
μ
g
r
a
d
v
+
S
M
y
∂
ρ
w
∂
t
+
d
i
v
ρ
w
u
=
-
∂
p
∂
z
+
d
i
v
μ
g
r
a
d
w
+
S
M
zwherep is the pressure and μ is the viscosity.u,v, and w are the flow velocities in thex-,y-, andz- directions, respectively. Similarly, S
M
x, S
M
y, and S
M
z are three components of the constant source term.If introducing a general variableϕ, (4) and (5) can be expressed as the following form:(6)
∂
ρ
ϕ
∂
t
+
d
i
v
ρ
ϕ
u
=
d
i
v
Γ
g
r
a
d
ϕ
+
Swhere Γ is the diffusion coefficient andϕ stands for velocityu,v, or w 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 left-hand side and the diffusive term and the source term, respectively, on the right-hand side.


## 2.2. Three-Dimensional 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 unknownu,v,p at the center of that CV and at the centers of neighbouring CVs [31]. A typical three-dimensional control volume element is shown in Figure 2. Positioned mid-way between adjacent nodes, the six faces of the control volume are labelledn, w,n,s,t, andb, which represent East, West, North, South, Top, and Bottom. The centroid of the element is identified by the nodal pointP and its neighbours in a three-dimensional 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.Figure 2
A cell in three dimensions and its neighbouring nodes.


## 2.3. Analytical Expressions forz-Derivatives of Velocity

The two sides of (2) and (3) are integrated twice with respect toz, which derives the formulation as follows:(7)
u
=
∂
p
∂
x
1
2
μ
z
2
-
h
2
4
(8)
v
=
∂
p
∂
y
1
2
μ
z
2
-
h
2
4with the boundary condition ofu =v = 0 at the fracture wall. Thez-derivatives of flow velocities at the boundary are expressed as(9)
∂
u
∂
z
z
=
h
/
2
=
-
∂
u
∂
z
z
=
-
h
/
2
=
∂
p
∂
x
h
2
μ
(10)
∂
v
∂
z
z
=
h
/
2
=
-
∂
v
∂
z
z
=
-
h
/
2
=
∂
p
∂
y
h
2
μAccording to the cubic law for Newtonian flow in the fracture [32], the averaged velocities in thex- andy-directions are separately given as(11)
u
-
=
-
∂
p
∂
x
h
2
12
μ
(12)
v
-
=
-
∂
p
∂
y
h
2
12
μRespectively substituting (11) and (12) into (9) and (10) yields the following equations:(13)
∂
u
∂
z
z
=
h
/
2
=
-
∂
u
∂
z
z
=
-
h
/
2
=
-
6
h
u
-
(14)
∂
v
∂
z
z
=
h
/
2
=
-
∂
v
∂
z
z
=
-
h
/
2
=
-
6
h
v
-If the general variableϕ
- is introduced, (13)-(14) can be usefully written in the following form:(15)
∂
ϕ
∂
z
z
=
h
/
2
=
-
6
h
ϕ
-
(16)
∂
ϕ
∂
z
z
=
-
h
/
2
=
6
h
ϕ
-


## 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 pointP. For the three-dimensional control volume defined in Figure 2, fromttot+Δtthis gives(17)
∫
t
t
+
Δ
t
∫
Δ
V
∂
ρ
ϕ
∂
t
d
V
d
t
+
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
ρ
u
ϕ
d
V
d
t
=
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
μ
g
r
a
d
ϕ
d
V
d
t
+
∫
t
t
+
Δ
t
∫
Δ
V
S
d
V
d
tTaking the average velocitiesu
- and v
- as unknown, the volume integral of the rate of change term can be written as(18)
∫
t
t
+
Δ
t
∫
Δ
V
∂
ρ
ϕ
∂
t
d
V
d
t
=
ρ
ϕ
-
P
-
ϕ
-
P
0
Δ
VThe 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 conditionw
t
=
w
b
=
0, the convective term can be transformed into the following equations:(19)
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
ρ
u
ϕ
d
V
d
t
=
∫
t
t
+
Δ
t
ρ
u
-
e
A
e
ϕ
-
e
-
ρ
u
-
w
A
w
ϕ
-
w
+
ρ
v
-
n
A
n
ϕ
-
n
-
ρ
v
-
s
A
s
ϕ
-
s
d
tSimilarly, 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:(20)
∫
t
t
+
Δ
t
∫
Δ
V
d
i
v
μ
g
r
a
d
ϕ
d
V
d
t
=
∫
t
t
+
Δ
t
μ
e
A
e
ϕ
-
E
-
ϕ
-
P
δ
x
P
E
-
μ
w
A
w
ϕ
-
P
-
ϕ
-
W
δ
x
W
P
+
μ
n
A
n
ϕ
-
N
-
ϕ
-
P
δ
y
P
N
-
μ
s
A
s
ϕ
-
P
-
ϕ
-
S
δ
y
S
P
+
-
μ
t
A
t
6
ϕ
-
P
h
-
μ
b
A
b
6
ϕ
-
P
h
d
tThe source term is approximated by means of a linear form:(21)
∫
t
t
+
Δ
t
∫
Δ
V
S
d
V
d
t
=
∫
t
t
+
Δ
t
S
c
+
S
P
ϕ
-
P
Δ
V
d
tThe fully implicit scheme is applied to conduct the time integral of Eq. (17). Then two variablesF andDare defined to represent the convective mass flux per unit area and diffusion conductance at cell faces, which are calculated with the following formulae:(22)
F
w
=
ρ
u
-
w
A
w
,
F
e
=
ρ
u
-
e
A
e
F
n
=
ρ
v
-
n
A
n
,
F
s
=
ρ
v
-
s
A
s
D
w
=
μ
w
A
w
δ
x
W
P
,
D
e
=
μ
e
A
e
δ
x
P
E
D
n
=
μ
n
A
n
δ
y
P
N
,
D
s
=
μ
s
A
s
δ
y
S
PBy substituting (18)-(22) into (17), the three-dimensional numerical calculation is transformed into two-dimensional calculation, and the corresponding discretized momentum equation can be rearranged as(23)
a
P
ϕ
-
P
=
a
E
ϕ
-
E
+
a
W
ϕ
-
W
+
a
N
ϕ
-
N
+
a
S
ϕ
-
S
+
bwhere(24)
a
P
=
a
W
+
a
E
+
a
S
+
a
N
+
μ
A
6
h
+
μ
A
6
h
+
F
e
-
F
w
+
F
n
-
F
s
+
a
P
0with(25)
a
P
0
=
ρ
Δ
V
Δ
t
,
b
=
S
Δ
V
+
a
P
0
ϕ
-
P
0The neighbour coefficients of (23) for the exponential scheme are shown as follows:(26)
a
W
=
D
w
max
⁡
0
,
1
-
0.1
F
w
D
w
5
+
max
⁡
0
,
F
w
a
E
=
D
e
max
⁡
0
,
1
-
0.1
F
e
D
e
5
+
max
⁡
0
,
-
F
e
a
S
=
D
s
max
⁡
0
,
1
-
0.1
F
s
D
s
5
+
max
⁡
0
,
F
s
a
B
=
D
b
max
⁡
0
,
1
-
0.1
F
b
D
b
5
+
max
⁡
0
,
F
bThe discretized form of continuity equation for the control volume gives(27)
ρ
u
A
e
-
ρ
u
A
w
+
ρ
v
A
n
-
ρ
v
A
s
=
0Based on SIMPLE algorithm, the discretized continuity equation will be turned into a pressure correction equation. A variable, pressure correctionp
′, is firstly defined as the difference between the correct pressure and the guessed pressure. The pressure correction values of nodal pointP and its neighbours are denoted as p
P
′, p
E
′, p
W
′, p
N
′, and p
S
′, the subscript of which indicates the position of the node in the three-dimensional geometry.The pressure correction equation is treated in a similar manner. Its discretized equation is written as(28)
a
P
′
p
P
′
=
a
E
′
p
E
′
+
a
W
′
p
W
′
+
a
N
′
p
N
′
+
a
S
′
p
S
′
+
b
P
′where(29)
b
P
′
=
ρ
u
-
A
w
-
ρ
u
-
A
e
+
ρ
v
-
A
s
-
ρ
v
-
A
n
,
a
P
′
=
a
E
′
+
a
W
′
+
a
N
′
+
a
S
′with the coefficients as follows:(30)
a
E
′
=
ρ
d
A
e
,
a
W
′
=
ρ
d
A
w
a
N
′
=
ρ
d
A
n
,
a
S
′
=
ρ
d
A
sFrom the above calculation process, it can be seen that the three-dimensional flow has been transformed into a two-dimensional problem by introducing the hypothesis of Hele-Shaw model and deducing the analytical expression for partial derivatives of velocity with respect toz. The present model taking the fracture aperture into account only needs to make the two-dimensional numerical calculation; therefore it is called the quasi-3D 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, Level-Set 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:(31)
λ
x
,
y
,
t
=
0
x
,
y
∈
Ω
1
1
x
,
y
∈
Ω
2where Ω
1 and Ω
2 represent the area occupied by air and grout, respectively. The scalar function satisfies the following equation:(32)
∂
λ
∂
t
+
u
∂
λ
∂
x
+
v
∂
λ
∂
y
=
0wheret is time.u andv are the velocities inx-andy-direction, respectively.The computational domain is divided into multiple noncoincident cells. The volume of a cell is denoted as△
V
i
,
j. The integral ofλ(x,y,t) over each cellI
i,j isF
i,j, which can be written as(33)
F
i
,
j
=
1
Δ
V
i
,
j
∫
I
i
,
j
λ
x
,
y
,
t
d
VThe functionF
i,j is called the fractional volume function, which obeys the following convection equation:(34)
∂
F
∂
t
+
u
∂
F
∂
x
+
v
∂
F
∂
y
=
0Equation (34) is the VOF equation and its solving results can be divided into three cases: ifF= 1, only the grout is present in the cell; ifF= 0, the cell is full of air; if 0
<
F
<
1, the phase interface can be in the cell. According to the value ofF, 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 Hirt-Nichols 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 ofF 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 functionF, velocity vectoru, and pressurep are initialized at all cells and the boundary conditions are defined.
(2)
Material properties at all cells are initialized according to the timet and fractional volume functionF.
(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 AppendixA.The fracture is 1 m long and 1 m wide with an aperture of 2 mm, as shown in Figure3. A no-slip 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/m3 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/m3. 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 quasi-3D numerical model and the commercial code ANSYS-Fluent. In the new model the element length in thez-direction 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 B150-PLUS
(iii)
Memory: 16 GB (DDR4 2401MHz), 400Figure 3
Geometric sizes of the fracture in Case 1.Figure 4
Computational domain of Case 1.In the quasi-3D 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 thez-direction 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 ANSYS-Fluent 17.0 under the same boundary condition. The simulated domain was divided into 500064 hexahedral cells with an edge length of 1 mm in thez-direction. 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.Figures5–7 display the advancement of the grout flow front as a function of time in the fracture. The diffusion front positions of grout att = 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.Figure 5
Simulated grout diffusion process by quasi-3D model in Case 1.
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sFigure 6
Simulated grout diffusion process by Fluent (2 layers model).
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sFigure 7
Simulated grout diffusion process by Fluent (6 layers model).
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sThe numerical results and analytical solution of the grout diffusion radius are shown in Table1. According to these data, the variation curves of grout diffusion radius with time are plotted in Figure 8. The average relative errors of the Quasi-3D 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.Table 1
Variations of the diffusion radius with time calculated by three models in Case 1.
Time
t(s)
Diffusion radiusr
Analytical(m)
Quasi-3D model(m)
Fluent_21(m)
Fluent_6l(m)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
5
0.1732
0.1718
0.1760
0.1740
0.81
1.62
0.46
10
0.2449
0.2431
0.2450
0.2460
0.73
0.04
0.45
15
0.3000
0.2981
0.3010
0.3020
0.63
0.33
0.67
20
0.3464
0.3444
0.3460
0.3480
0.58
0.12
0.46
25
0.3873
0.3862
0.3880
0.3890
0.28
0.18
0.44
30
0.4242
0.4236
0.4260
0.4250
0.16
0.40
0.16
35
0.4582
0.4579
0.4580
0.4590
0.09
0.07
0.15
40
0.4899
0.4899
0.4910
0.4910
0.00
0.22
0.22
Average relative errors (%)
0.41
0.37
0.38Figure 8
Comparison of diffusion radius variation with time among three models in Case 1.The grout pressure at different positions calculated by the three models is presented in Table2 att= 20 s and Table 3 att= 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% att= 20 s, respectively, while 4.25% and 1.90% att= 40 s, respectively. By contrast, 6 layers model calculated by ANSYS-Fluent has a lower accuracy. Att= 20 s, the relative errors of the pressure and velocity are up to 7.30% and 5.63% separately. Att= 40 s, they are 7.52% and 5.77%, respectively. Then the solution accuracy of 2 layers model calculated by ANSYS-Fluent is further reduced. The maximum relative error is more than 33%.Table 2
Grout pressure distribution over radius att= 20 s in Case 1.
Radius
r (m)
Pressure (p-p
0)
Analytical(kPa)
Quasi-3D model(kPa)
Fluent_2l(kPa)
Fluent_6l(kPa)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.025
0.9463
0.9372
0.6191
0.8772
0.96
34.58
7.30
0.050
0.6968
0.6896
0.4611
0.6527
1.03
33.83
6.33
0.075
0.5508
0.5415
0.3654
0.5167
1.69
33.66
6.19
0.100
0.4473
0.4383
0.2969
0.4193
2.01
33.62
6.26
0.125
0.3669
0.3574
0.2436
0.3436
2.62
33.62
6.38
0.150
0.3013
0.2921
0.2000
0.2816
3.05
33.62
6.54
0.175
0.2458
0.2375
0.1631
0.2292
3.38
33.65
6.75
0.200
0.1978
0.1906
0.1311
0.1838
3.64
33.72
7.08
0.225
0.1553
0.1494
0.1029
0.1447
3.80
33.74
6.83
0.250
0.1174
0.1127
0.0777
0.1092
4.00
33.82
6.98
0.275
0.0831
0.0796
0.0548
0.0776
4.21
34.06
6.62
0.300
0.0518
0.0495
0.0340
0.0482
4.44
34.36
6.95
Average relative errors (%)
2.90
33.86
6.68Table 3
Grout pressure distribution over radius att= 40 s in Case 1.
Radius
r (m)
Pressure (p-p
0)
Analytical(kPa)
Quasi-3D model(kPa)
Fluent_2l(kPa)
Fluent_6l(kPa)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.05
0.8216
0.8173
0.5450
0.7709
0.52
33.67
6.17
0.10
0.5720
0.5660
0.3807
0.5375
1.05
33.44
6.03
0.15
0.4261
0.4197
0.2838
0.3998
1.50
33.40
6.17
0.20
0.3225
0.3160
0.2150
0.3019
2.02
33.33
6.39
0.25
0.2422
0.2356
0.1615
0.2259
2.73
33.32
6.73
0.30
0.1765
0.1724
0.1178
0.1640
2.32
33.26
7.08
0.35
0.1211
0.1180
0.0808
0.1124
2.56
33.28
7.18
0.40
0.0730
0.0712
0.0488
0.0677
2.47
33.15
7.26
0.45
0.0306
0.0293
0.0206
0.0283
4.25
32.68
7.52
Average relative errors (%)
2.16
33.28
6.73Table 4
Average flow velocity at different radiuses att= 20 s in Case 1.
Radius
r (m)
Average flow velocityu
Analytical(cm/s)
Quasi-3D model(cm/s)
Fluent_2l(cm/s)
Fluent_6l(cm/s)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.025
12.0000
11.9976
8.0095
11.4062
0.02
33.25
4.95
0.050
6.0000
6.1148
3.9987
5.6874
1.91
33.36
5.21
0.075
4.0000
4.0639
2.6656
3.7897
1.6
33.36
5.26
0.100
3.0000
3.0323
1.9997
2.8418
1.08
33.36
5.27
0.125
2.4000
2.4126
1.5993
2.2731
0.52
33.36
5.29
0.150
2.0000
2.0105
1.3328
1.8940
0.52
33.36
5.30
0.175
1.7143
1.7205
1.1424
1.6230
0.36
33.36
5.33
0.200
1.5000
1.5061
0.9996
1.4196
0.41
33.36
5.36
0.225
1.3333
1.3377
0.8885
1.2613
0.33
33.36
5.40
0.250
1.2000
1.2035
0.7996
1.1344
0.29
33.37
5.47
0.275
1.0909
1.0893
0.7269
1.0303
0.15
33.37
5.56
0.300
1.0000
0.9813
0.6663
0.9437
1.87
33.37
5.63
Average relative errors (%)
0.76
33.35
5.34Table 5
Average flow velocity at different radiuses att=40 s in Case 1.
Radius
r (m)
Average flow velocityu
Analytical(cm/s)
Quasi-3D model(cm/s)
Fluent_2l(cm/s)
Fluent_6l(cm/s)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.05
6.0000
6.1140
3.9987
5.6874
1.90
33.36
5.21
0.10
3.0000
3.0308
1.9992
2.8418
1.03
33.36
5.27
0.15
2.0000
2.0083
1.3328
1.8943
0.42
33.36
5.28
0.20
1.5000
1.5035
0.9996
1.4204
0.23
33.36
5.31
0.25
1.2000
1.2017
0.7997
1.1360
0.14
33.36
5.33
0.30
1.0000
1.0007
0.6664
0.9462
0.07
33.36
5.38
0.35
0.8571
0.8570
0.5712
0.8103
0.01
33.36
5.46
0.40
0.7500
0.7494
0.4998
0.7081
0.08
33.36
5.59
0.45
0.6667
0.6729
0.4442
0.6282
0.93
33.37
5.77
Average relative errors (%)
0.53
33.36
5.40Figure 9
Grout pressure distribution over radius at different times in Case 1 (a)t= 20 s; (b)t= 40 s.
(a)
(b)Figure 10
Distribution of average flow velocity over the radius at different times in Case 1 (a)t= 20 s; (b)t= 40 s.
(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.Table6 shows an overview of the quasi-3D 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 quasi-3D 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.Table 6
Description of three models and average relative errors of the solution.
Model
Number of the grid cell
Size of grid cell (mm)
Total computing time (h)
Average relative errors (%)
Diffusion radius
Pressure distribution
Velocity distribution
Quasi-3D model
2500
10
0.5
0.41
2.58
0.66
2 layers model
500064
1<△x <3.5
48
0.37
33.61
33.36
6 layers model
1500264
0.3<△x <2.1
432
0.38
6.70
5.36
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 parallel-plate fracture under constant grouting pressure, and its analytical solution is listed in AppendixB.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/m3 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 quasi-3D model, the number of the grid unit is 2500.Figure11 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.Figure 11
Change of grout diffusion front positions with time.
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
s
(e)
t
=
50
sTable7 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 quasi-3D 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%.Table 7
Variations of the diffusion radius and the pumping rate with time in Case 2.
Time
t(s)
Pumping rateq
Diffusion radiusr
Analytical(mL/s)
Quasi-3D model(mL/s)
Relative errors(%)
Analytical(mL/s)
Quasi-3D model(mL/s)
Relative errors(%)
10
1.3761
1.3981
1.60
0.2246
0.2218
1.25
20
1.2747
1.2855
0.85
0.3040
0.3051
0.36
30
1.2216
1.2181
0.29
0.3634
0.3654
0.55
40
1.1864
1.1800
0.54
0.4127
0.4155
0.68
50
1.1604
1.1581
0.20
0.4557
0.4576
0.42
Average relative errors (%)
0.70
0.65Figure 12
Variations of the pumping rate with time in Case 2.Figure 13
Variations of the diffusion radius with time in Case 2.Tables8 and 9 show the contrast of the radial pressure and average velocity calculated by the quasi-3D model and the analytical model att= 20 s andt= 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.Table 8
Distribution of the pressure and average flow velocity over the radius att= 20 s in Case 2.
Radius
r (m)
Pressure (p-p
0)
Average flow velocityu
Analytical(kPa)
Quasi-3D model(kPa)
Relative errors(%)
Analytical(cm/s)
Quasi-3D model(cm/s)
Relative errors(%)
0.05
87.8813
89.1206
1.41
4.0575
4.1630
2.60
0.10
54.1318
54.3447
0.39
2.0288
2.0781
2.43
0.15
34.3895
34.2431
0.43
1.3525
1.3763
1.76
0.20
20.3822
20.0348
1.70
1.0144
1.0277
1.31
0.25
9.5172
9.2332
2.98
0.8115
0.8184
0.85
0.30
0.6399
0.6108
4.55
0.6763
0.6950
2.77
Average relative errors (%)
1.91
1.95Table 9
Distribution of the pressure and average flow velocity over the radius att= 50 s in Case 2.
Radius
r (m)
Pressure (p-p
0)
Average flow velocityu
Analytical(kPa)
Quasi-3D model(kPa)
Relative errors(%)
Analytical(cm/s)
Quasi-3D model(cm/s)
Relative errors(%)
0.05
97.9463
99.3736
1.46
3.6935
3.769
2.04
0.10
67.2244
67.8087
0.87
1.8468
1.8865
2.15
0.15
49.2533
49.5546
0.61
1.2312
1.2504
1.56
0.20
36.5025
36.6361
0.37
0.9234
0.9355
1.31
0.25
26.6123
26.6317
0.07
0.7387
0.747
1.12
0.30
18.5314
18.4711
0.33
0.6156
0.6213
0.93
0.35
11.699
11.5866
0.96
0.5276
0.531
0.64
0.40
5.7806
5.6527
2.21
0.4617
0.4596
0.45
0.45
0.5602
0.5306
5.28
0.4104
0.4202
2.39
Average relative errors (%)
1.35
1.40Figure 14
Pressure distribution over radius at different times in Case 2 (a)t= 20 s; (b)t= 50 s.
(a)
(b)Figure 15
Distribution of average flow velocity over radius at different times in Case 2 (a)t= 20 s; (b)t= 50 s.
(a)
(b)


### 3.2. Comparison with Experimental Results

To further verify the applicability of the quasi-3D 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 Figure16. 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.Figure 16
Schematic diagram of experimental set-up: (a) top view (b) lateral view.
(a)
(b)The quasi-3D 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 Figure17. The domain was divided into 50 cells in thex-direction and 30 cells in they-direction.Figure 17
Calculation region in Case 3.Figure18 shows the grout diffusion region and morphology att = 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 quasi-3D model are in good accordance with the experimental and analytical results, demonstrating that the presented model can accurately simulate the grouting diffusion processTable 10
Diffusion radius at different times for the three methods.
Time
t(s)
Diffusion radiusr(m)
Analytical
Quasi-3D model
Numerical relative errors (%)
Experimental
Experimental relative errors (%)
5
0.0892
0.0902
1.12
0.0908
1.79
10
0.1262
0.1273
0.87
0.1255
0.55
15
0.1545
0.1556
0.71
0.1558
0.84
20
0.1784
0.1792
0.45
0.1780
0.22
25
0.1995
0.2002
0.35
0.1976
0.95
30
0.2185
0.2192
0.32
0.2167
0.82
Average relative errors (%)
0.64
0.86Figure 18
Variation of the grout diffusion region with time in the experiment.
(a)
t
=
5
s
(b)
t
=
10
s
(c)
t
=
15
s
(d)
t
=
20
s
(e)
t
=
25
s
(f)
t
=
30
sFigure 19
Comparison of the grout diffusion region and morphology between numerical and experimental results.
(a)
t
=
5
s
(b)
t
=
10
s
(c)
t
=
15
s
(d)
t
=
20
s
(e)
t
=
25
s
(f)
t
=
30
s


## 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 AppendixA.The fracture is 1 m long and 1 m wide with an aperture of 2 mm, as shown in Figure3. A no-slip 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/m3 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/m3. 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 quasi-3D numerical model and the commercial code ANSYS-Fluent. In the new model the element length in thez-direction 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 B150-PLUS
(iii)
Memory: 16 GB (DDR4 2401MHz), 400Figure 3
Geometric sizes of the fracture in Case 1.Figure 4
Computational domain of Case 1.In the quasi-3D 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 thez-direction 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 ANSYS-Fluent 17.0 under the same boundary condition. The simulated domain was divided into 500064 hexahedral cells with an edge length of 1 mm in thez-direction. 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.Figures5–7 display the advancement of the grout flow front as a function of time in the fracture. The diffusion front positions of grout att = 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.Figure 5
Simulated grout diffusion process by quasi-3D model in Case 1.
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sFigure 6
Simulated grout diffusion process by Fluent (2 layers model).
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sFigure 7
Simulated grout diffusion process by Fluent (6 layers model).
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sThe numerical results and analytical solution of the grout diffusion radius are shown in Table1. According to these data, the variation curves of grout diffusion radius with time are plotted in Figure 8. The average relative errors of the Quasi-3D 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.Table 1
Variations of the diffusion radius with time calculated by three models in Case 1.
Time
t(s)
Diffusion radiusr
Analytical(m)
Quasi-3D model(m)
Fluent_21(m)
Fluent_6l(m)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
5
0.1732
0.1718
0.1760
0.1740
0.81
1.62
0.46
10
0.2449
0.2431
0.2450
0.2460
0.73
0.04
0.45
15
0.3000
0.2981
0.3010
0.3020
0.63
0.33
0.67
20
0.3464
0.3444
0.3460
0.3480
0.58
0.12
0.46
25
0.3873
0.3862
0.3880
0.3890
0.28
0.18
0.44
30
0.4242
0.4236
0.4260
0.4250
0.16
0.40
0.16
35
0.4582
0.4579
0.4580
0.4590
0.09
0.07
0.15
40
0.4899
0.4899
0.4910
0.4910
0.00
0.22
0.22
Average relative errors (%)
0.41
0.37
0.38Figure 8
Comparison of diffusion radius variation with time among three models in Case 1.The grout pressure at different positions calculated by the three models is presented in Table2 att= 20 s and Table 3 att= 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% att= 20 s, respectively, while 4.25% and 1.90% att= 40 s, respectively. By contrast, 6 layers model calculated by ANSYS-Fluent has a lower accuracy. Att= 20 s, the relative errors of the pressure and velocity are up to 7.30% and 5.63% separately. Att= 40 s, they are 7.52% and 5.77%, respectively. Then the solution accuracy of 2 layers model calculated by ANSYS-Fluent is further reduced. The maximum relative error is more than 33%.Table 2
Grout pressure distribution over radius att= 20 s in Case 1.
Radius
r (m)
Pressure (p-p
0)
Analytical(kPa)
Quasi-3D model(kPa)
Fluent_2l(kPa)
Fluent_6l(kPa)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.025
0.9463
0.9372
0.6191
0.8772
0.96
34.58
7.30
0.050
0.6968
0.6896
0.4611
0.6527
1.03
33.83
6.33
0.075
0.5508
0.5415
0.3654
0.5167
1.69
33.66
6.19
0.100
0.4473
0.4383
0.2969
0.4193
2.01
33.62
6.26
0.125
0.3669
0.3574
0.2436
0.3436
2.62
33.62
6.38
0.150
0.3013
0.2921
0.2000
0.2816
3.05
33.62
6.54
0.175
0.2458
0.2375
0.1631
0.2292
3.38
33.65
6.75
0.200
0.1978
0.1906
0.1311
0.1838
3.64
33.72
7.08
0.225
0.1553
0.1494
0.1029
0.1447
3.80
33.74
6.83
0.250
0.1174
0.1127
0.0777
0.1092
4.00
33.82
6.98
0.275
0.0831
0.0796
0.0548
0.0776
4.21
34.06
6.62
0.300
0.0518
0.0495
0.0340
0.0482
4.44
34.36
6.95
Average relative errors (%)
2.90
33.86
6.68Table 3
Grout pressure distribution over radius att= 40 s in Case 1.
Radius
r (m)
Pressure (p-p
0)
Analytical(kPa)
Quasi-3D model(kPa)
Fluent_2l(kPa)
Fluent_6l(kPa)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.05
0.8216
0.8173
0.5450
0.7709
0.52
33.67
6.17
0.10
0.5720
0.5660
0.3807
0.5375
1.05
33.44
6.03
0.15
0.4261
0.4197
0.2838
0.3998
1.50
33.40
6.17
0.20
0.3225
0.3160
0.2150
0.3019
2.02
33.33
6.39
0.25
0.2422
0.2356
0.1615
0.2259
2.73
33.32
6.73
0.30
0.1765
0.1724
0.1178
0.1640
2.32
33.26
7.08
0.35
0.1211
0.1180
0.0808
0.1124
2.56
33.28
7.18
0.40
0.0730
0.0712
0.0488
0.0677
2.47
33.15
7.26
0.45
0.0306
0.0293
0.0206
0.0283
4.25
32.68
7.52
Average relative errors (%)
2.16
33.28
6.73Table 4
Average flow velocity at different radiuses att= 20 s in Case 1.
Radius
r (m)
Average flow velocityu
Analytical(cm/s)
Quasi-3D model(cm/s)
Fluent_2l(cm/s)
Fluent_6l(cm/s)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.025
12.0000
11.9976
8.0095
11.4062
0.02
33.25
4.95
0.050
6.0000
6.1148
3.9987
5.6874
1.91
33.36
5.21
0.075
4.0000
4.0639
2.6656
3.7897
1.6
33.36
5.26
0.100
3.0000
3.0323
1.9997
2.8418
1.08
33.36
5.27
0.125
2.4000
2.4126
1.5993
2.2731
0.52
33.36
5.29
0.150
2.0000
2.0105
1.3328
1.8940
0.52
33.36
5.30
0.175
1.7143
1.7205
1.1424
1.6230
0.36
33.36
5.33
0.200
1.5000
1.5061
0.9996
1.4196
0.41
33.36
5.36
0.225
1.3333
1.3377
0.8885
1.2613
0.33
33.36
5.40
0.250
1.2000
1.2035
0.7996
1.1344
0.29
33.37
5.47
0.275
1.0909
1.0893
0.7269
1.0303
0.15
33.37
5.56
0.300
1.0000
0.9813
0.6663
0.9437
1.87
33.37
5.63
Average relative errors (%)
0.76
33.35
5.34Table 5
Average flow velocity at different radiuses att=40 s in Case 1.
Radius
r (m)
Average flow velocityu
Analytical(cm/s)
Quasi-3D model(cm/s)
Fluent_2l(cm/s)
Fluent_6l(cm/s)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.05
6.0000
6.1140
3.9987
5.6874
1.90
33.36
5.21
0.10
3.0000
3.0308
1.9992
2.8418
1.03
33.36
5.27
0.15
2.0000
2.0083
1.3328
1.8943
0.42
33.36
5.28
0.20
1.5000
1.5035
0.9996
1.4204
0.23
33.36
5.31
0.25
1.2000
1.2017
0.7997
1.1360
0.14
33.36
5.33
0.30
1.0000
1.0007
0.6664
0.9462
0.07
33.36
5.38
0.35
0.8571
0.8570
0.5712
0.8103
0.01
33.36
5.46
0.40
0.7500
0.7494
0.4998
0.7081
0.08
33.36
5.59
0.45
0.6667
0.6729
0.4442
0.6282
0.93
33.37
5.77
Average relative errors (%)
0.53
33.36
5.40Figure 9
Grout pressure distribution over radius at different times in Case 1 (a)t= 20 s; (b)t= 40 s.
(a)
(b)Figure 10
Distribution of average flow velocity over the radius at different times in Case 1 (a)t= 20 s; (b)t= 40 s.
(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.Table6 shows an overview of the quasi-3D 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 quasi-3D 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.Table 6
Description of three models and average relative errors of the solution.
Model
Number of the grid cell
Size of grid cell (mm)
Total computing time (h)
Average relative errors (%)
Diffusion radius
Pressure distribution
Velocity distribution
Quasi-3D model
2500
10
0.5
0.41
2.58
0.66
2 layers model
500064
1<△x <3.5
48
0.37
33.61
33.36
6 layers model
1500264
0.3<△x <2.1
432
0.38
6.70
5.36
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 parallel-plate fracture under constant grouting pressure, and its analytical solution is listed in AppendixB.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/m3 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 quasi-3D model, the number of the grid unit is 2500.Figure11 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.Figure 11
Change of grout diffusion front positions with time.
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
s
(e)
t
=
50
sTable7 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 quasi-3D 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%.Table 7
Variations of the diffusion radius and the pumping rate with time in Case 2.
Time
t(s)
Pumping rateq
Diffusion radiusr
Analytical(mL/s)
Quasi-3D model(mL/s)
Relative errors(%)
Analytical(mL/s)
Quasi-3D model(mL/s)
Relative errors(%)
10
1.3761
1.3981
1.60
0.2246
0.2218
1.25
20
1.2747
1.2855
0.85
0.3040
0.3051
0.36
30
1.2216
1.2181
0.29
0.3634
0.3654
0.55
40
1.1864
1.1800
0.54
0.4127
0.4155
0.68
50
1.1604
1.1581
0.20
0.4557
0.4576
0.42
Average relative errors (%)
0.70
0.65Figure 12
Variations of the pumping rate with time in Case 2.Figure 13
Variations of the diffusion radius with time in Case 2.Tables8 and 9 show the contrast of the radial pressure and average velocity calculated by the quasi-3D model and the analytical model att= 20 s andt= 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.Table 8
Distribution of the pressure and average flow velocity over the radius att= 20 s in Case 2.
Radius
r (m)
Pressure (p-p
0)
Average flow velocityu
Analytical(kPa)
Quasi-3D model(kPa)
Relative errors(%)
Analytical(cm/s)
Quasi-3D model(cm/s)
Relative errors(%)
0.05
87.8813
89.1206
1.41
4.0575
4.1630
2.60
0.10
54.1318
54.3447
0.39
2.0288
2.0781
2.43
0.15
34.3895
34.2431
0.43
1.3525
1.3763
1.76
0.20
20.3822
20.0348
1.70
1.0144
1.0277
1.31
0.25
9.5172
9.2332
2.98
0.8115
0.8184
0.85
0.30
0.6399
0.6108
4.55
0.6763
0.6950
2.77
Average relative errors (%)
1.91
1.95Table 9
Distribution of the pressure and average flow velocity over the radius att= 50 s in Case 2.
Radius
r (m)
Pressure (p-p
0)
Average flow velocityu
Analytical(kPa)
Quasi-3D model(kPa)
Relative errors(%)
Analytical(cm/s)
Quasi-3D model(cm/s)
Relative errors(%)
0.05
97.9463
99.3736
1.46
3.6935
3.769
2.04
0.10
67.2244
67.8087
0.87
1.8468
1.8865
2.15
0.15
49.2533
49.5546
0.61
1.2312
1.2504
1.56
0.20
36.5025
36.6361
0.37
0.9234
0.9355
1.31
0.25
26.6123
26.6317
0.07
0.7387
0.747
1.12
0.30
18.5314
18.4711
0.33
0.6156
0.6213
0.93
0.35
11.699
11.5866
0.96
0.5276
0.531
0.64
0.40
5.7806
5.6527
2.21
0.4617
0.4596
0.45
0.45
0.5602
0.5306
5.28
0.4104
0.4202
2.39
Average relative errors (%)
1.35
1.40Figure 14
Pressure distribution over radius at different times in Case 2 (a)t= 20 s; (b)t= 50 s.
(a)
(b)Figure 15
Distribution of average flow velocity over radius at different times in Case 2 (a)t= 20 s; (b)t= 50 s.
(a)
(b)


## 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 AppendixA.The fracture is 1 m long and 1 m wide with an aperture of 2 mm, as shown in Figure3. A no-slip 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/m3 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/m3. 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 quasi-3D numerical model and the commercial code ANSYS-Fluent. In the new model the element length in thez-direction 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 B150-PLUS
(iii)
Memory: 16 GB (DDR4 2401MHz), 400Figure 3
Geometric sizes of the fracture in Case 1.Figure 4
Computational domain of Case 1.In the quasi-3D 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 thez-direction 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 ANSYS-Fluent 17.0 under the same boundary condition. The simulated domain was divided into 500064 hexahedral cells with an edge length of 1 mm in thez-direction. 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.Figures5–7 display the advancement of the grout flow front as a function of time in the fracture. The diffusion front positions of grout att = 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.Figure 5
Simulated grout diffusion process by quasi-3D model in Case 1.
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sFigure 6
Simulated grout diffusion process by Fluent (2 layers model).
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sFigure 7
Simulated grout diffusion process by Fluent (6 layers model).
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
sThe numerical results and analytical solution of the grout diffusion radius are shown in Table1. According to these data, the variation curves of grout diffusion radius with time are plotted in Figure 8. The average relative errors of the Quasi-3D 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.Table 1
Variations of the diffusion radius with time calculated by three models in Case 1.
Time
t(s)
Diffusion radiusr
Analytical(m)
Quasi-3D model(m)
Fluent_21(m)
Fluent_6l(m)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
5
0.1732
0.1718
0.1760
0.1740
0.81
1.62
0.46
10
0.2449
0.2431
0.2450
0.2460
0.73
0.04
0.45
15
0.3000
0.2981
0.3010
0.3020
0.63
0.33
0.67
20
0.3464
0.3444
0.3460
0.3480
0.58
0.12
0.46
25
0.3873
0.3862
0.3880
0.3890
0.28
0.18
0.44
30
0.4242
0.4236
0.4260
0.4250
0.16
0.40
0.16
35
0.4582
0.4579
0.4580
0.4590
0.09
0.07
0.15
40
0.4899
0.4899
0.4910
0.4910
0.00
0.22
0.22
Average relative errors (%)
0.41
0.37
0.38Figure 8
Comparison of diffusion radius variation with time among three models in Case 1.The grout pressure at different positions calculated by the three models is presented in Table2 att= 20 s and Table 3 att= 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% att= 20 s, respectively, while 4.25% and 1.90% att= 40 s, respectively. By contrast, 6 layers model calculated by ANSYS-Fluent has a lower accuracy. Att= 20 s, the relative errors of the pressure and velocity are up to 7.30% and 5.63% separately. Att= 40 s, they are 7.52% and 5.77%, respectively. Then the solution accuracy of 2 layers model calculated by ANSYS-Fluent is further reduced. The maximum relative error is more than 33%.Table 2
Grout pressure distribution over radius att= 20 s in Case 1.
Radius
r (m)
Pressure (p-p
0)
Analytical(kPa)
Quasi-3D model(kPa)
Fluent_2l(kPa)
Fluent_6l(kPa)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.025
0.9463
0.9372
0.6191
0.8772
0.96
34.58
7.30
0.050
0.6968
0.6896
0.4611
0.6527
1.03
33.83
6.33
0.075
0.5508
0.5415
0.3654
0.5167
1.69
33.66
6.19
0.100
0.4473
0.4383
0.2969
0.4193
2.01
33.62
6.26
0.125
0.3669
0.3574
0.2436
0.3436
2.62
33.62
6.38
0.150
0.3013
0.2921
0.2000
0.2816
3.05
33.62
6.54
0.175
0.2458
0.2375
0.1631
0.2292
3.38
33.65
6.75
0.200
0.1978
0.1906
0.1311
0.1838
3.64
33.72
7.08
0.225
0.1553
0.1494
0.1029
0.1447
3.80
33.74
6.83
0.250
0.1174
0.1127
0.0777
0.1092
4.00
33.82
6.98
0.275
0.0831
0.0796
0.0548
0.0776
4.21
34.06
6.62
0.300
0.0518
0.0495
0.0340
0.0482
4.44
34.36
6.95
Average relative errors (%)
2.90
33.86
6.68Table 3
Grout pressure distribution over radius att= 40 s in Case 1.
Radius
r (m)
Pressure (p-p
0)
Analytical(kPa)
Quasi-3D model(kPa)
Fluent_2l(kPa)
Fluent_6l(kPa)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.05
0.8216
0.8173
0.5450
0.7709
0.52
33.67
6.17
0.10
0.5720
0.5660
0.3807
0.5375
1.05
33.44
6.03
0.15
0.4261
0.4197
0.2838
0.3998
1.50
33.40
6.17
0.20
0.3225
0.3160
0.2150
0.3019
2.02
33.33
6.39
0.25
0.2422
0.2356
0.1615
0.2259
2.73
33.32
6.73
0.30
0.1765
0.1724
0.1178
0.1640
2.32
33.26
7.08
0.35
0.1211
0.1180
0.0808
0.1124
2.56
33.28
7.18
0.40
0.0730
0.0712
0.0488
0.0677
2.47
33.15
7.26
0.45
0.0306
0.0293
0.0206
0.0283
4.25
32.68
7.52
Average relative errors (%)
2.16
33.28
6.73Table 4
Average flow velocity at different radiuses att= 20 s in Case 1.
Radius
r (m)
Average flow velocityu
Analytical(cm/s)
Quasi-3D model(cm/s)
Fluent_2l(cm/s)
Fluent_6l(cm/s)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.025
12.0000
11.9976
8.0095
11.4062
0.02
33.25
4.95
0.050
6.0000
6.1148
3.9987
5.6874
1.91
33.36
5.21
0.075
4.0000
4.0639
2.6656
3.7897
1.6
33.36
5.26
0.100
3.0000
3.0323
1.9997
2.8418
1.08
33.36
5.27
0.125
2.4000
2.4126
1.5993
2.2731
0.52
33.36
5.29
0.150
2.0000
2.0105
1.3328
1.8940
0.52
33.36
5.30
0.175
1.7143
1.7205
1.1424
1.6230
0.36
33.36
5.33
0.200
1.5000
1.5061
0.9996
1.4196
0.41
33.36
5.36
0.225
1.3333
1.3377
0.8885
1.2613
0.33
33.36
5.40
0.250
1.2000
1.2035
0.7996
1.1344
0.29
33.37
5.47
0.275
1.0909
1.0893
0.7269
1.0303
0.15
33.37
5.56
0.300
1.0000
0.9813
0.6663
0.9437
1.87
33.37
5.63
Average relative errors (%)
0.76
33.35
5.34Table 5
Average flow velocity at different radiuses att=40 s in Case 1.
Radius
r (m)
Average flow velocityu
Analytical(cm/s)
Quasi-3D model(cm/s)
Fluent_2l(cm/s)
Fluent_6l(cm/s)
Relative errors (%)
Quasi-3D model
Fluent_2l
Fluent_6l
0.05
6.0000
6.1140
3.9987
5.6874
1.90
33.36
5.21
0.10
3.0000
3.0308
1.9992
2.8418
1.03
33.36
5.27
0.15
2.0000
2.0083
1.3328
1.8943
0.42
33.36
5.28
0.20
1.5000
1.5035
0.9996
1.4204
0.23
33.36
5.31
0.25
1.2000
1.2017
0.7997
1.1360
0.14
33.36
5.33
0.30
1.0000
1.0007
0.6664
0.9462
0.07
33.36
5.38
0.35
0.8571
0.8570
0.5712
0.8103
0.01
33.36
5.46
0.40
0.7500
0.7494
0.4998
0.7081
0.08
33.36
5.59
0.45
0.6667
0.6729
0.4442
0.6282
0.93
33.37
5.77
Average relative errors (%)
0.53
33.36
5.40Figure 9
Grout pressure distribution over radius at different times in Case 1 (a)t= 20 s; (b)t= 40 s.
(a)
(b)Figure 10
Distribution of average flow velocity over the radius at different times in Case 1 (a)t= 20 s; (b)t= 40 s.
(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.Table6 shows an overview of the quasi-3D 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 quasi-3D 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.Table 6
Description of three models and average relative errors of the solution.
Model
Number of the grid cell
Size of grid cell (mm)
Total computing time (h)
Average relative errors (%)
Diffusion radius
Pressure distribution
Velocity distribution
Quasi-3D model
2500
10
0.5
0.41
2.58
0.66
2 layers model
500064
1<△x <3.5
48
0.37
33.61
33.36
6 layers model
1500264
0.3<△x <2.1
432
0.38
6.70
5.36
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 parallel-plate fracture under constant grouting pressure, and its analytical solution is listed in AppendixB.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/m3 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 quasi-3D model, the number of the grid unit is 2500.Figure11 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.Figure 11
Change of grout diffusion front positions with time.
(a)
t
=
10
s
(b)
t
=
20
s
(c)
t
=
30
s
(d)
t
=
40
s
(e)
t
=
50
sTable7 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 quasi-3D 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%.Table 7
Variations of the diffusion radius and the pumping rate with time in Case 2.
Time
t(s)
Pumping rateq
Diffusion radiusr
Analytical(mL/s)
Quasi-3D model(mL/s)
Relative errors(%)
Analytical(mL/s)
Quasi-3D model(mL/s)
Relative errors(%)
10
1.3761
1.3981
1.60
0.2246
0.2218
1.25
20
1.2747
1.2855
0.85
0.3040
0.3051
0.36
30
1.2216
1.2181
0.29
0.3634
0.3654
0.55
40
1.1864
1.1800
0.54
0.4127
0.4155
0.68
50
1.1604
1.1581
0.20
0.4557
0.4576
0.42
Average relative errors (%)
0.70
0.65Figure 12
Variations of the pumping rate with time in Case 2.Figure 13
Variations of the diffusion radius with time in Case 2.Tables8 and 9 show the contrast of the radial pressure and average velocity calculated by the quasi-3D model and the analytical model att= 20 s andt= 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.Table 8
Distribution of the pressure and average flow velocity over the radius att= 20 s in Case 2.
Radius
r (m)
Pressure (p-p
0)
Average flow velocityu
Analytical(kPa)
Quasi-3D model(kPa)
Relative errors(%)
Analytical(cm/s)
Quasi-3D model(cm/s)
Relative errors(%)
0.05
87.8813
89.1206
1.41
4.0575
4.1630
2.60
0.10
54.1318
54.3447
0.39
2.0288
2.0781
2.43
0.15
34.3895
34.2431
0.43
1.3525
1.3763
1.76
0.20
20.3822
20.0348
1.70
1.0144
1.0277
1.31
0.25
9.5172
9.2332
2.98
0.8115
0.8184
0.85
0.30
0.6399
0.6108
4.55
0.6763
0.6950
2.77
Average relative errors (%)
1.91
1.95Table 9
Distribution of the pressure and average flow velocity over the radius att= 50 s in Case 2.
Radius
r (m)
Pressure (p-p
0)
Average flow velocityu
Analytical(kPa)
Quasi-3D model(kPa)
Relative errors(%)
Analytical(cm/s)
Quasi-3D model(cm/s)
Relative errors(%)
0.05
97.9463
99.3736
1.46
3.6935
3.769
2.04
0.10
67.2244
67.8087
0.87
1.8468
1.8865
2.15
0.15
49.2533
49.5546
0.61
1.2312
1.2504
1.56
0.20
36.5025
36.6361
0.37
0.9234
0.9355
1.31
0.25
26.6123
26.6317
0.07
0.7387
0.747
1.12
0.30
18.5314
18.4711
0.33
0.6156
0.6213
0.93
0.35
11.699
11.5866
0.96
0.5276
0.531
0.64
0.40
5.7806
5.6527
2.21
0.4617
0.4596
0.45
0.45
0.5602
0.5306
5.28
0.4104
0.4202
2.39
Average relative errors (%)
1.35
1.40Figure 14
Pressure distribution over radius at different times in Case 2 (a)t= 20 s; (b)t= 50 s.
(a)
(b)Figure 15
Distribution of average flow velocity over radius at different times in Case 2 (a)t= 20 s; (b)t= 50 s.
(a)
(b)


## 3.2. Comparison with Experimental Results

To further verify the applicability of the quasi-3D 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 Figure16. 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.Figure 16
Schematic diagram of experimental set-up: (a) top view (b) lateral view.
(a)
(b)The quasi-3D 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 Figure17. The domain was divided into 50 cells in thex-direction and 30 cells in they-direction.Figure 17
Calculation region in Case 3.Figure18 shows the grout diffusion region and morphology att = 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 quasi-3D model are in good accordance with the experimental and analytical results, demonstrating that the presented model can accurately simulate the grouting diffusion processTable 10
Diffusion radius at different times for the three methods.
Time
t(s)
Diffusion radiusr(m)
Analytical
Quasi-3D model
Numerical relative errors (%)
Experimental
Experimental relative errors (%)
5
0.0892
0.0902
1.12
0.0908
1.79
10
0.1262
0.1273
0.87
0.1255
0.55
15
0.1545
0.1556
0.71
0.1558
0.84
20
0.1784
0.1792
0.45
0.1780
0.22
25
0.1995
0.2002
0.35
0.1976
0.95
30
0.2185
0.2192
0.32
0.2167
0.82
Average relative errors (%)
0.64
0.86Figure 18
Variation of the grout diffusion region with time in the experiment.
(a)
t
=
5
s
(b)
t
=
10
s
(c)
t
=
15
s
(d)
t
=
20
s
(e)
t
=
25
s
(f)
t
=
30
sFigure 19
Comparison of the grout diffusion region and morphology between numerical and experimental results.
(a)
t
=
5
s
(b)
t
=
10
s
(c)
t
=
15
s
(d)
t
=
20
s
(e)
t
=
25
s
(f)
t
=
30
s


## 4. Conclusion

Numerical simulation is always an important way to study the diffusion mechanism of grout in the rock fracture. Based on the traditional three-dimensional FVM, the aim of this paper is to develop a quasi-3D 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 toz at fracture walls were deduced by introducing the Hele-Shaw 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 three-dimensional diffusion of the grout in the parallel fracture can be transformed into a two-dimensional problem. Combined with the VOF method of tracking the grout front, a quasi-3D simulation method was established to predict the flow and diffusion behavior of grout in the fracture.
(ii)
The quasi-3D 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 ANSYS-Fluent 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 quick-setting 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 quasi-3D model may be further extended to a three-dimensional intersected fracture network in the next phase; thus it is possible to simulate the grout injection in a more complex environment [36, 37].



---

*Source: 4139616-2019-07-25.xml*