#### Abstract

The flow of cement-based slurry within concentric annular geometry is a major problem of study in the field of engineering, especially regarding the prevention of water ingress in an annulus formed during shaft construction utilizing the artificial freezing technique in China. In this study, an analytical governing equation of motion for the axial flow of an incompressible Newtonian fluid in a long concentric annulus is derived under the condition of high groundwater pressure. A stepwise calculation method is proposed to describe the grouting process based on two injection modes, namely, flow rate control and pressure control. The injection time is divided into a series of time segments; correspondingly, the grouted zone is subdivided into infinitesimal elements. Some key parameters, such as the location, dimension, slurry viscosity, and pressure gradient of each element, can be obtained using the developed MATLAB program. On this basis, the distribution of pressure and slurry viscosity in the grouted zone and the variations in injection pressure at the grouting point and grouting flow rate are determined. Two injection mode cases are investigated to reveal the grout propagation in the concentric annulus. Finally, numerical simulations are conducted and employed to validate and calibrate the calculated results. The results obtained by the present stepwise calculation method show good agreement with the numerical results.

#### 1. Introduction

Many water-burst hazards have occurred in coal mines in China during the past few decades, causing casualties and significant property losses [1–17]. To prevent water ingress, the grouting technique has been extensively and effectively applied to water-bearing strata in underground constructions for more than two hundred years [18–21].

In the northwestern area of China, many typical North China coalfields, such as the Ordos Coalfield, include large-scale coalmines. In the Dongsheng area, the Ordos basin coal deposit developed in the late Jurassic and is covered by Cretaceous sandstones and fine sandstones [22, 23]. Due to their short period of diagenesis, these strata are characterized by low strengths, weak cementation, and water-rich rock, and full-depth artificial freezing technology is widely used to construct deep mine shafts in this area [24–26]. During the full-depth artificial freezing process, the depth of the frozen barrier is greater than that of the shaft to prevent water ingress, as shown in Figure 1. After the shaft construction is finished, the frozen rock thaws, and the annular space between the freezing borehole and freezing pipe opens again. The upper aquifers directly come into contact with the lower aquifers, inducing frozen wall breakage and shaft flooding [27, 28]. In addition, during the freezing-thawing process, the surrounding rock is damaged due to the frost force from the groundwater, and water movement occurs in the damage zone, causing the formation of fissures and flow of water around the annular space. The long-term interaction of water and rock aggravates the shaft water hazard, possibly leading to shaft lining collapse. Pregrouting or grouting has been one of the main effective options used to control groundwater inrush from annular space [24, 25].

In China, dozens of engineering cases of grouting annular channels around shafts have been implemented successfully [25, 26]. Based on this rich engineering experience, researchers have proposed several effective grouting techniques, mainly including back wall grouting and perforation grouting. Due to the maturity of back wall grouting technology, its convenient operation, and its high adaptability, the method of back wall grouting is employed widely in the treatment of frozen boreholes. The following design for back wall grouting was introduced as follows: (1) gathering in situ data; (2) determining the appropriate stratum for a grouting field trial based on the analysis of accessible geological exploration data; (3) calculating the spatial position of the freezing holes according to the inclinometer data of the boreholes; (4) defining the grouting parameters such as the slurry recipe, injection pressure and injection borehole position; (5) implementing the grouting process under the guidance of the grouting design; and (6) stopping grouting and utilizing the geophysical exploration and/or hydrogeological analysis to evaluate the grouting effect.

Although the efficacy of the back wall grouting method has been verified by many field trials, some of the main operation parameters, such as injection pressure, penetration length, and grouting flow rate, are difficult to select in the grouting design stage. This limitation stems from a lack of understanding of the propagation mechanism of the slurry used in annular conduit grouting.

Over the past few decades, many researchers have studied the mechanical characteristics of fluid flow in coaxial cylinders used across various industries, including oil and gas exploration, geothermal heating, and transportation of suspensions. Various scholars have carried out an abundance of research on the axial laminar flow of fluids, injected under pressure, within the annular conduit during the drilling of oil wells. The fluid velocity and pressure drop of rheological fluids in an ideal concentric conduit have been widely discussed. Numerical and analytical solutions of axial annular flow have been presented for Newtonian and non-Newtonian fluids [29–33]. These research findings have important implications for the analysis of slurry propagation in annular channels. However, the annulus grouting design is determined by three aspects: the rheological properties of the slurry, the geometry of the annulus, and the injection modes, i.e., flow rate control grouting or pressure control grouting. In the stages of grouting design and execution, the fundamental parameters, including the injection pressure, grout advance distance, and injection hole layout, are governed by these factors. Therefore, describing the control effect of these influencing factors on slurry propagation in a concentric annulus becomes a critical problem to prevent water inrush from freezing holes.

In this paper, first, theoretical foundations of the basic equations of a rheological slurry and the governing equation of axial laminar flow in a concentric annulus are expressed. Second, a stepwise calculation method is proposed to describe the grouting process according to pressure control and flow control conditions, and some crucial parameters, such as slurry viscosity, injection pressure, and pressure distribution in the grouted zone, are analyzed through case studies. Then, the results of the developed calculation method are validated and compared with numerical simulations; finally, the conclusions are drawn.

#### 2. Mathematical Formulation of Grout Propagating in a Concentric Annulus

##### 2.1. Assumptions

For the mathematical formulation, the change in the grout viscosity with injection time is considered with the following assumptions:(1)The grout penetrates axially into a concentric annulus as a laminar flow(2)The grout is assumed to be an incompressible Newtonian fluid(3)The concentric annulus is vertical, with constant radii of the inner cylinder and the outer cylinder, and the inertial force is negligible(4)A nonslip condition is assumed for both surfaces of the annulus(5)It is assumed that the groundwater pressure in the annular grouting area is constant and high due to a deep burial condition, and the gravity force of the grout is negligible(6)The chemical reaction time of the grout is assumed to be equal to the injection time, that is, the time required for the grout to move from the grouting pump to the annular passage is negligible(7)There is no physical mixing and/or chemical reaction between the slurry and water

##### 2.2. Axial Laminar Flow Equation

To simulate the flow state of a grout in a rock mass, it is important to calculate the rheology of the injected fluid. In most cases, the material behavior of the injected fluid can be described with various rheological models, which are divided into two general categories: Newtonian fluids and Bingham fluids [34]. In general, a cement-based grout can be considered a Newtonian fluid when the cement slurry is mixed with a sodium silicate slurry or the relevant water to cement (w/c) ratio of the cementitious grout ranges from 2 to 10 [21, 35, 36]. Therefore, the relationship between the shear stress and velocity is [37]where *τ* is the shear stress, is a time-dependent function of the dynamic viscosity, *t* is the injection time, is the velocity, *h* is the distance perpendicular to the flow direction, and *dv/dh* is the velocity gradient.

It is assumed that grout flows axially in a concentric annulus. The geometry of cement-based grout in a concentric annulus is depicted in Figure 2. As shown in Figure 2(a), the radii of the inner cylinder and outer cylinder are *R*_{i} and *R*_{o}, respectively, and the radius of the injection borehole is *r*_{0}. The amount of grout in the annulus is limited by its pressure front. The pressure at the grout front *p*_{f} is assumed to be equal to the initial water pressure *p*_{w}, and the pressure at the intersection of the annulus and a horizontal plane through the injection borehole axis is the grouting pressure *p*_{g}.

**(a)**

**(b)**

In general, the grout is propelled by pumps, mixed in the mixer, and injected into the concentric annulus through the injection borehole. It takes a relatively short time for the slurry to move through the borehole; hence, the chemical reaction time of the grout mix could be considered equal to the injection time. In addition, this study focused on the axial flow in a concentric annulus; thus, the slurry is assumed to instantaneously fill the initial annular space with a height of 2*r*_{0} once the grout injection begins. Consequently, the slurry travel time in the annular channel is assumed to be equivalent to the injection time, namely, the chemical reaction time of the grout mix.

As shown in Figure 2(b), a section that crosses through the injection hole and the upward flow in the annulus is selected and a rectangular coordinate system is established. It is assumed that the axial flow in the annulus is symmetrical across the center of the open annulus between the pipe and hole surfaces, i.e., , if the concentric annular channel is rather narrow [38, 39]. During grouting of a vertical concentric annulus under the condition of high groundwater pressure, it is assumed that the grout axially penetrates the annulus as a laminar flow. A grout element in the center of the open annulus between the pipe and hole surfaces is selected to analyze the force state, as shown in Figure 2(b). At a position *z*, i.e., the distance between the injection point and the grout element, the equilibrium equation of the grout element, can be expressed aswhere is the length of the grout element, is the pressure, is the pressure increment, and is the location of the symmetric plane of slurry flow.

In addition, the velocity profile of Figure 2(b) shows that the axial flow is divided into an inner section and an outer section .

###### 2.2.1. Average Grout Velocity in Area

The condition of exists in area , and the relationship between the grout velocity and shear rate satisfies the following formula:

Combining equations (2) and (3), the velocity gradient of the grout in area can be expressed as

It is considered that grout spread in the concentric annulus is incompressible and laminar. The velocity of the grout front is obtained by integrating equation (4). The integration constants are eliminated using the boundary condition due to a nonslip condition imposed on the inner cylinder surface. Then, the velocity profile of the grout front in area can be calculated as

The mean slurry velocity in area can be expressed by integrating equation (5) in the interval as follows:where is the average slurry velocity in the grouted zone .

###### 2.2.2. Average Grout Velocity in Area

Based on a similar calculation process, the average slurry velocity in area can be expressed as follows:where is the average slurry velocity in the grouted zone .

Combining equations (6) and (7), the average slurry velocity in the concentric annulus can be given as

In general, the annular space around the freezing holes is filled with high-pressure water due to the deep burial condition. In some coalmines in the Ordos area, the groundwater pressure in the grouting area is approximately 5 MPa [24, 40], which varies with the buried depth. The variation of ambient water pressure complicates the analysis of the pressure at the grout front and makes it difficult to study the axial grout flow. In addition, when the slurry moves vertically, gravity is an important factor affecting grout distribution in the initial grout stage, resulting in asymmetric grout distribution in the vertical two directions, which hinders the time-independent study of the entire grout process. Consequently, the gravity force of the groundwater and slurry can be neglected during analysis of the grouting process, i.e., the water pressure remains constant in the grouting area. Thus, the slurry injected into the annular channel would be distributed into two parts with equal mass, one moving upward and one downward. For concentric annular grouting, the principle of mass conservation is obeyed, and the relationship between the mean grout velocity and flow rate can be given aswhere *q* is the grouting flow rate of the upward flow and downward flow along the concentric annulus.

Inserting equation (9) into equation (8) , the governing equation of the pressure gradient along the grout propagation direction in the grouted zone is calculated as follows:

#### 3. Stepwise Calculation Method for Concentric Annulus Grouting

##### 3.1. Injection Modes

According to equation (10), it is concluded that the main aspects influencing the grouting process include the geometry of the annulus and the rheology, flow rate, and injection pressure of the slurry material. In terms of a specific project case, the former two factors could be determined based on testing the rheological properties of the grout materials and reviewing the shaft construction data; however, the latter two factors would vary with the injection time, which is determined by adjusting the grouting sequence. This adjustment is achieved by simultaneously controlling the injection flow rate and pressure [41], as illustrated in Figure 3.

The grout evolution in Figure 3, considering first flow rate control and then pressure control with a break of a few minutes between the end of the flow rate control and the beginning of flow rate control, shows that the grout is injected into a fracture. Some key parameters, such as the flow rate, injection pressure, and penetration length, are calculated based on realistic values of the theological grout parameters. If the grout is injected under the flow rate control condition, the grouting flow rate is constant. With increasing grouting time, the grout advances and the penetration rate declines. Simultaneously, similar phenomena are revealed in the injection pressure. However, the grout propagates dramatically when the grout is injected under the pressure control condition. With the further increase in grouting time, the grout continues to advance. Conversely, the grouting flow decreases with a variable decay rate. Thus, the grouting process is difficult to describe by an analytical solution for rheological slurry materials. Kim et al. proposed a stepwise numerical calculation method to analyze the grouting process in a porous medium, considering the variable viscosity and the decrease in porosity resulting from filtration. Zhang et al. proposed a stepwise method to describe the grouting process of quick-setting slurry in a rock fracture. Based on the theory of iterative study [42–44], a stepwise calculation method is established to describe the annular channel grouting process, considering the two different grouting modes, i.e., the flow rate control condition and injection pressure control condition.

##### 3.2. Stepwise Calculation Method Corresponding to Injection Modes

According to the injection time, at each time step , the grouted area could be subdivided into a series of infinitesimal elements that fill the size of the concentric cylindrical annulus. Each element could be represented by , where *j* represents the number of time steps corresponding to the total injection time *T* (, *j* = 1, 2, 3, 4, …) and *i* indicates the time step when the slurry of element is injected (*i* = 1, 2, 3, …, *j*) [44]. For any grouting time , the element neighboring the injection borehole is expressed as and that at the grout front is marked as , as shown in Figure 4.

Moreover, the location of an element could be restricted by a starting boundary and end boundary, which are far from and close to the injection point and represented by and , respectively. Thus, the height of the infinitesimal element could be expressed as [44]

To describe the propagation of the slurry in the grouted zone, the variations in the dimension, pressure, and slurry viscosity in every element should be calculated first, and then the grouting overpressure (grouting pressure minus groundwater pressure) at the injection point, grouting flow rate, and grout advance could be determined. When the grouting process is performed under flow rate control conditions, the grouting flow rate is constant. The position of a grout particle at a given injection time can be obtained according to the mass conservation equation. Furthermore, the location of a fixed point in the grouted zone can be expressed as a function of the injection time and grouting flow rate, and the viscosity of the slurry can be calculated according to the rheology of the slurry. Based on the governing equation of grout propagation, the pressure can be calculated at any point in the grouted zone [44]. However, as the grout is injected under the pressure control condition, the injection pressure remains constant and the flow rate decreases gradually, meaning that the location and geometry of the annular element cannot be directly estimated from the mass conservation law. The force analysis would be used to determine the flow rate first, and the location and geometry of an annular element can be calculated [44]. Then, the grouting overpressure can be determined by the sum of pressure segments of all elements in the grouted zone at any injection time, and the overpressure distribution along the grout advance would be obtained according to the axial laminar flow equation.

As presented previously, the determination of the grout take of each element is the critical step to describe the process of grout propagation [44]. With known flow rate control conditions, the grouting flow rate can be expressed as a function of injection time, and the grout take during the time interval can be obtained as follows:where is the grouting flow rate function of injection time.

In the time interval , the grout mix is injected and fills element . According to the mass conservation law, the grout take in the grouting durability time is equal to the element volume. The relationship between the starting position and end position of element can be determined by using

As the time step Δ*t* is set to be rather small, the slurry viscosity within an infinitesimal element is assumed to linearly increase. In other words, the viscosity of element can be used in place of the average value of the viscosity of the grout particles at the starting and end boundaries and . The viscosity of element can then be expressed as [44]where and are the slurry viscosity at the starting and end boundaries of element , respectively. Because the travel time of the grout mix in the injection borehole is negligible, the times required for the grout particles to move from the injection point to the starting and end boundaries can be calculated as and , respectively. The slurry viscosity at the advance is calculated from

Consequently, the slurry viscosity at the injection point is considered to be always equal to the initial value. In other words, the slurry viscosity at the injection point can be written as .

Similarly, the location of element can also be indicated by the average value of the location of the two boundaries, expressed as [44]where and are the propagation distances to the starting boundary and end boundary of element , respectively. For the element neighboring the injection point, the penetration length at the end boundary is always equal to the radius of the injection borehole and can be written as .

To drive grout movement, the pressure at the end boundary must be greater than that at the starting boundary . The pressure difference is equal to the pressure segment of element . In this case, at each computational step, the relationship between the boundary pressure and pressure segment can be expressed aswhere is the pressure segment and can be written as the following equation by combining equations (12)–(17):

The grout is confined with the area from the injection point to the grout front along the penetration length. Moreover, the pressure at the grout front and the pressure at the injection point are assumed to be equal to the groundwater pressure and grouting pressure, respectively:where is the pressure at the grout front and is the groundwater pressure.where is the pressure at the injection point and is the grouting pressure.

If the grout was injected into the annular channel under constant pressure, the grouting flow rate decreases with injection time; therefore, the location and dimension of the element cannot be calculated directly by using mass balance equations. Thus, the grouting flow rate at any injection time should be obtained by a force analysis of the axial flow, and then the location and size of the element can be calculated.

In the initial stage of pressure control grouting, the flow rate is generally quite high to maintain a constant injection pressure. Thus, the slurry would be injected into the annulus at a very high speed to fill a region around the injection point. The width of the region, which is called the initial element, is greater than the radius of the injection hole and can be expressed as [44]. It is considered that the injection time of the initial element is neglected. Consequently, the viscosity of the slurry in this region can be assumed to be equal to the initial value, which can be written as . In the initial element, the pressure linearly decreases with grouting advance; thus, the initial grouting flow rate can be expressed as the following formula according to (10):where is the initial grouting flow rate under the pressure control injection mode.

Moreover, the flow rate at injection time can be described as the average grouting flow rate in the time interval , and thus the relationship of the grout take and grouting flow rate can be given as

At injection time , the pressure gradient in any annular element satisfies equation (18), and the sum of the pressure segments in each element is equal to the difference between the grouting pressure at the injection point and the pressure of the surrounding water and can be written as [44]where is the grouting pressure at injection time .

Combining equations (18) and (23), the grouting flow rate at injection time can be expressed as follows:where is the grouting flow rate when the slurry is injected into the annular channel at time .

As presented previously, the crucial parameters, such as the grout take, the slurry viscosity, and the pressure segment in any element, of element in the grouted area can be calculated according to equations (12)–(25) under the two injection modes of flow rate control and pressure control, based on the stepwise calculation method developed by using MATLAB. Hence, the grouting process under different injection modes in a concentric annulus can be described accurately.

Figure 5 provides a flow chart of the stepwise calculation method, in which the flow rate control mode and the pressure control mode are both used. First, some key initial parameters, such as the initial slurry viscosity and rheological parameters, constant flow rate under flow control mode or constant grouting pressure under pressure control mode, total injection time, time step, radius of the injection borehole, and geometry of the annular channel and groundwater pressure, are input into the stepwise calculation process. Second, the grout take during the time interval is calculated by using (12) when the grout is injected with a constant flow rate, or the initial flow rate and grout take are determined according to equations (21) and (22), respectively, when the condition of injection pressure is adapted. Third, the initial boundary conditions, e.g., the initial grouted zone position , the pressure at the grout front , and the slurry viscosity of the initial grouted zone , are set in the stepwise algorithm. Fourth, the slurry viscosities at the boundaries of the element and their average value are calculated. Moreover, the grout flow rate is subsequently obtained using (24) if the pressure control injection mode is applied. Then, the average grout penetration length, location, and pressure segment of element are calculated. Finally, the calculated step is updated, and the grouting overpressure at the injection point, the penetration length, and the distributions of overpressure and slurry viscosity in the grout propagation direction at each injection time are acquired.

#### 4. Application of the Proposed Stepwise Method

##### 4.1. Calculation Parameters

To verify the applicability of the proposed stepwise calculation method describing the grouting process in a concentric annulus, in this study, a computer program was developed using MATLAB, considering the two injection modes of flow control and pressure control. The slurry types selected for this calculation include two cement-based grouts: cement and sodium silicate grout (CSG) and cement slurry (CG). In general, CSG is a quick-setting slurry and will gel within tens of seconds to several minutes, depending on two main parameters, namely, the w/c ratio and the volume ratio of the cement slurry to the sodium silicate solution (C : S ratio) [21].

In this study, the w/c ratio and C : S ratio are chosen as 1 and 1 : 1, respectively. The rheological properties of a cement grout and its stability are discussed in many papers based on experimental tests [45, 46]. As a reference grout material, the cement slurry is assumed to behave as a Newtonian fluid, maintaining a constant viscosity throughout the grouting process. The calculation parameters of the two calculation conditions are summarized in Table 1.

##### 4.2. Calculation Results of the Flow Rate Control Injection Mode

The calculation results of the flow rate control injection mode are shown in Figure 6.

**(a)**

**(b)**

**(c)**

Figure 6(a) presents the distribution of the slurry viscosity with grout advance at three different injection times: 15 s, 35 s, and 55 s. For the CSG, the viscosity of the slurry in the grouted zone varies significantly with the penetration length, generally increasing from the region close to the injection point to that far from the injection point. That is, the slurry viscosity quickly increases with the increase in grout advance when the quick-setting grout is injected into the annular channel. In addition, the slurry viscosity and penetration distance curves for different injection times are partially coincident near the injection hole, indicating that the slurry viscosity is stable at a given position in the injected region [21]. This phenomenon can be attributed to the constant grouting flow rate and the fixed section size of the annulus. The injection time of the grout suspensions moving from the injection point to a fixed position is constant, and the viscosity of the grout suspensions corresponding to the fixed position maintains a certain value. In terms of the CG selected in this study, the viscosity of the slurry does not vary with injection time and grout advance, so that the corresponding curve of the slurry viscosity with penetration length is a straight line with a constant slope.

Figure 6(b) presents the overpressure distribution curves of the CSG in the grouted zone at three different times: 15 s, 35 s, and 55 s. The spatial distributions of overpressure for the CG are shown for comparison. For the grout with time-varying viscosity (i.e., CSG), the pressure attenuates far from the injection point, but after a certain penetration length, the decay rate increases, which is controlled by the change law of the slurry viscosity [21]. As the chemical reaction time between the cement slurry and sodium silicate solution is less than the twist time, approximately 16 s, the rate of increase in the slurry viscosity is low, which leads to the slow increase in the resistance and low corresponding pressure gradient in the grout penetration direction. Once the chemical reaction time exceeds the twist time, the rate of increase in the slurry viscosity increases and the driving force required for the slurry to move increases; this change is characterized by the increase in the pressure gradient. The larger the area of grout with a rapidly growing viscosity within the concentric annular channel, the greater the pressure gradient that exists near the grout front. For constant-viscosity grout (i.e., CG), with increasing penetration length, the pressure decreases linearly in the grouted zone, which is very different from the results of the CSG and is ascribed to the superimposed effect of axial diffusion and its constant viscosity.

Figure 6(c) indicates the variations in the grouting overpressure with injection time for the time-dependent viscosity slurry and constant-viscosity slurry. Between the two types of grouts, the grouting overpressure changes throughout the grouting process in very different ways. With increasing grouting durability time, the injection pressure of the CSG increases and the pressure rate increases greatly after a certain grouting time due to the constant flow rate and time-varying viscosity of CSG. Conversely, the injection pressure increases linearly with the grouting time when CG is used, but the injection pressure would reduce tens to hundreds of times compared to that of the quick-setting slurry.

##### 4.3. Calculation Results of Pressure Control Injection Mode

In the grouting process with the pressure control mode, the constant pressure was selected as 10 kPa. To avoid an excessive flow rate in the initial grouting period when the grout was injected, an appropriate treatment method was used: the injection pressure was manually and gradually increased from 0 to 10 kPa in the injection time interval 0∼1 s.

The calculation results of the pressure control injection mode are illustrated in Figure 7. Figure 7(a) presents the penetration length and grouting flow rate with injection time. As the slurry with a time-varying viscosity is injected into the concentric annulus, the grouting flow rate decreases sharply from a high initial value, approximately 0.25 m^{3}/h, and asymptotically decreases to less than 0.001 m^{3}/h at approximately 16 s, i.e., the switch time. Moreover, the variation in flow rate is significantly influenced by the grout advance. When the grouting time is less than the switch time, the slurry spreads effectively, and the penetration length increases slightly during the grout process. However, as the grouting time exceeds the switch time, the grouting process would enter an “ineffective stage,” in which the grout diffusion almost stops. In terms of the constant-viscosity slurry, the grouting rate also decreases significantly with increasing injection time but asymptotically decreases to a constant value of approximately 0.0016 m^{3}/h. This causes a continuous increase in grout advance during the injection process.

**(a)**

**(b)**

**(c)**

**(d)**

The distributions of grout viscosity at 15 s, 35 s, and 55 s are illustrated in Figure 7(b). The grout viscosity increases in the diffusion direction in the grouted zone; moreover, the viscosity at a fixed position increases with injection time. In other words, the three viscosity curves in Figure 7(b) are not coincident, which is strikingly different from the results of the grouting process adopting the constant flow rate mode. Additionally, with increasing injection time, the grouting flow rate decreases; consequently, the migration time of the grout mix from the injection borehole to a given position increases, causing the grout viscosity to increase at the same point.

The overpressure distributions of the injected CSG and CG in the grouted zone at 15 s, 35 s, and 55 s are shown in Figures 7(c) and 7(d), respectively. The overpressure decreases approximately linearly with increasing distance from the injection point, and the grouting pressure at the injection point remains constant throughout the grouting process, as presented by Zhang et al. [44]. Therefore, the overpressure decreases in approximately the propagation direction when the cement-based slurry is injected under constant pressure conditions, regardless of the type of diffusion path.

#### 5. Verification of the Numerical Simulation of Concentric Annulus Grouting

##### 5.1. Numerical Methodology and Choice of Numerical Parameters

The commercial finite element software Comsol Multiphysics was used to perform the numerical simulation of the upward flow in the vertical concentric annulus grouting process. A two-phase Darcy’s law model was adopted. In this model, the density of the water, the viscosity of the water, the density of the grout, and the space-dependent function of the slurry dynamic viscosity must be provided by the user. The calculation parameters used in the numerical analysis are presented in Table 2.

Analysis of the grouting within the concentric annulus was treated as a two-dimensional axisymmetric problem. A triangular mesh element was employed in this analysis. The mesh was extremely refined, and the maximum size of the mesh was 0.01 m. In addition, the right vertical boundary and the top boundary were subjected to pressure control and no-flux boundary condition, respectively. The bottom boundary was subjected to pressure control when constant pressure grouting was adopted or flow control when a constant flow velocity was applied, as illustrated in Figure 8. The transient solving method was adopted in the calculation process, and the variations in injection pressure at the inlet boundary, slurry penetration length, and pressure distribution in the grouted zone were simulated under constant flow rate conditions. Moreover, variations in the grouting flow rate, grout penetration length, and pressure distribution in the grouted zone were calculated during the pressure control grouting process.

**(a)**

**(b)**

In the initial stage of injection, the concentric annulus was filled with groundwater at a high pressure of 5 MPa. The volume fraction method was used to describe the distribution of slurry and water in the annulus, and the relationship between the volume fraction of the groundwater and that of the slurry at an injection time can be expressed as follows [46]:where is the volume fraction of the water and is the volume fraction of the slurry.

##### 5.2. Comparisons of the Calculated and Numerical Results

###### 5.2.1. Injection Mode of Flow Rate Control

Figure 9 illustrates the comparisons of the calculated and numerical simulation results under flow rate control conditions for the CSG.

**(a)**

**(b)**

The comparisons of the grouting overpressure at the injection point are shown in Figure 9(a), indicating that the result of the stepwise approach agrees with the numerical result. In the initial grouting process, some deviation between the calculated curve and the numerical curve is observed, and the maximum simulation error for CSG is 10%. However, with increasing injection time, the calculation error of the grouting pressure becomes less than 1%, which indicates that the stepwise algorithm can accurately describe the slurry propagation in the concentric annular channel. Figure 9(b) presents the overpressure (pressure in the grouted zone minus groundwater pressure) distribution of the CSG in the grouted zones at 15, 35, and 55 s. The predictions of the stepwise algorithm are in good agreement with the numerical results in the regions close to the injection point. However, a slight deviation between the calculated pressure and the simulated results is observed in the regions close to the grout front, where the numerical data are much larger, causing a maximum error of 5%. This deviation is because the boundary between the slurry and water is modeled as a transition zone rather than a sharp delimitation in the numerical simulation process.

###### 5.2.2. Injection Mode of Pressure Control

Figure 10 illustrates comparisons of the prediction and numerical simulation results under constant grouting pressure conditions for the cement slurry (CG).

**(a)**

**(b)**

The distributions of overpressure in the grouted zone filled with CG slurry at 15, 35, and 55 s are presented in Figure 10(a). Although there is good agreement between the predicted curves obtained by the stepwise algorithm and numerical data, with increasing distance from the injection point, insignificant deviations occur. Figure 10(b) shows that the predicted curves of the grouting penetration curve and flow rate obtained by the stepwise algorithm match the obtained curves obtained by the simulated results well; however, a slight deviation exists between the grouting penetration curves. The underestimated overpressure and grout advance in the grouted zone result from the differences between the assumptions in the analytical prediction and the processes that occurred in the numerical simulation. For example, a nonreaction or mixing assumption was made regarding the slurry and water in the step-time prediction; however, physical mixing may have occurred during the numerical simulation, producing a transition zone at the grout front and an underestimation of penetration length.

#### 6. Conclusions

The wall-back grouting technique has been extensively employed to prevent shaft flooding accidents induced by annular channel in freezing boreholes. In this paper, governing equation of concentric annulus grouting was deduced, and two kinds of slurry, cement, and sodium silicate grout (CSG) with a time-independent viscosity and cement grout (CG) with a constant viscosity were selected to study the grouting process in a concentric annulus.

The stepwise computational method, in which the grouted zone was divided into infinitesimal elements according to predetermined time step, was proposed to describe the concentric annulus grouting process. Some key parameters, such as the location, dimension, and slurry viscosity and pressure gradient of each element, were obtained by using the developed MATLAB program. On this basis, the distribution of pressure and slurry viscosity in the grouted zone and the variations in injection pressure at the grouting point and the grouting flow rate were determined.

Two grouting modes (flow rate control injection and pressure control grouting) were investigated by applying the proposed stepwise algorithm. And, the grout propagation in the concentric annulus was profoundly revealed. The results of the prediction of the proposed algorithm and the results of the numerical simulation performed in this study are compared and validated. In terms of the two grouting modes, the calculation data are in good agreement with those of the numerical simulation. However, slight deviations occur between the results of the two methods in the region close to the grout front because of the difference between the assumption in the stepwise calculation process and the processes that occur in the numerical simulation. For example, there exists a transition zone at the grout front in the numerical simulation due to the volume fraction method adopted, while this transition zone is not included in the calculation.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was supported by the Scientific Research Foundation of Shandong University of Science and Technology for Recruited Talents (no. 2017RCJJ030), the National Natural Science Foundation of China (no. 51509148), and the Shandong Provincial Natural Science Foundation of China (no.BS2015NJ010).