Lingcho Zou, Ulf H?knsson, Vldimir Cvetkovic

a Division of Resources,Energy and Infrastructure,Department of Sustainable Development,Environmental Science and Engineering,Royal Institute of Technology,Stockholm,10044, Sweden

b Division of Soil and Rock Mechanics, Department of Civil and Architectural Engineering, Royal Institute of Technology, Stockholm,10044, Sweden c Skanska Sweden AB, Stockholm,11274, Sweden

Keywords:Rock grouting Radial flow of Bingham fluids Plug flow region Force balance Energy dissipation Analytical solution

ABSTRACT Solutions for radial flow of a Bingham fluid are analyzed in this paper.It aims to eliminate confusions in the literature concerning the plug flow region in different solutions for analysis and design of grouting in rock fractures. The analyses based on the force balance equation reveal that the plug flow region in Bingham radial flow is independent of the fracture radius,and is not a growth function adapted from the solution of one-dimensional(1D)slit flow according to‘similarity’.Based on the shear stress distribution,we analytically proposed that a non-uniform plug flow region cannot exist. The Bingham fluid (grout)penetration and flowrate evolution as functions of grouting time are given using the correct expression for the plug flow region.The radius-independent plug flow region and the presented flowrate evolution equation are also verified numerically. For radial flow, the relative penetration length is equal to the relative width of plug flow region, which is the same as that for 1D channel flow. Discrepancies in analytical solutions for grout penetration and flowrate evolution were also illustrated.The clarification of the plug flow region and evaluation of discrepancies in analytical solutions presented in this work could simplify modeling and design of grouting in rock engineering applications.

1. Introduction

Cement grouting is widely used in rock engineering to seal rock fractures and limit groundwater flow.Modeling and analysis of the grout flow are important in design, execution and monitoring of grouting in fractured rocks (Wallner,1976; Lombardi,1985; Zhang et al., 2014; Stille, 2015; Sui et al., 2015; Li et al., 2016; Zhang et al., 2018). The mostly used cement grouts in practice are typically non-Newtonian fluids. In particular, the Bingham model has been frequently used to model grouting, due to its simplicity and physically based parameters(H?ssler,1991;H?kansson et al.,1992;H?kansson,1993; Rahman et al., 2015; Zou et al., 2018, 2019).

At present, simulation of rock grouting process remains a challenge, since the rock grouting process involves a complex fluid propagation in rock fractures with complex geometrical structures.In practice, design of rock grouting generally relies on analytical solutions based on simplified rheological properties for grouts and idealized geometrical conditions for single rock fractures. In order to derive analytical solutions for the Bingham type of grout propagating in single fractures,from an intersecting borehole, the process is generally idealized as two-dimensional (2D) radial flow between parallel planar disks (e.g. Wallner,1976; Lombardi,1985;Stille,2015;Sui et al.,2015;Zhang et al.,2018).Dai and Bird(1981)derived a solution of flowrate by adapting solution of onedimensional (1D) plane channel flow according to similarity. It was concluded that the plug flow region in the Bingham radial flow is an increasing function of the radial distance. By using the flowrate equation presented by Dai and Bird (1981), Gustafson et al.(2013) derived a solution for the relative penetration length of grout as a function of time,which lays a theoretical foundation for the real time grouting control (RTGC) method for grouting design(Kobayashi et al., 2008; Stille et al., 2009). This solution has been extended and widely used in grouting (e.g. Fransson et al., 2007;Kobayashi et al.,2008;Stille et al.,2012;Funehag and Th?rn,2014;Rafi and Stille, 2014; Stille, 2015). However, the aforementioned conclusion of the Bingham plug flow region presented by Dai and Bird (1981) is not directly derived from the mass and momentum equations and has led to discussions in the literature.The validity of the solution by Gustafson et al. (2013) was questioned by El Tani(2013).El Tani(2012)derived a solution for radial flow of a Bingham fluid and relative grout penetration with time,and determined the plug flow region using the energy dissipation equation.The results of El Tani(2012)showed that the plug flow region is independent of fracture radius at each time step,and the relative width of the plug flow region (i.e. ratio of the width of the plug flow region to the fracture aperture) equals the relative length of penetration (i.e.ratio of the penetration length to the maximum penetration length)when the grout is propagating, which is the same as that for 1D channel flow.

The diverging perceptions of the plug flow region directly lead to the two alternative analytical solutions for grout penetration(El Tani, 2012; Gustafson et al.,2013).Rafi and Stille (2015) compared two solutions with the results of the relative penetration length only for relative time at the range of less than 0.6(a dimensionless time defnied asby Gustafson et al.(2013),where t is the injection time,τ0is the yield stress,μ is the plastic viscosity,Pgis the injection pressure, and Pwis the in situ groundwater pressure).It was found that when relative time is less than 0.2,the two results are almost the same. Recently, El Tani and Stille (2017)discussed the two solutions without providing any verification or direct comparison of them. Ignoring the fundamentally different perceptions of the plug flow region and different solutions for the pressure,they conjectured that the two solutions are equal,which contradicts with El Tani’s perspectives (e.g. El Tani, 2012, 2013). In fact,a comprehensive verification of the two solutions has not been done yet,and the issue regarding the shape of plug flow region has not been fully resolved.

The objective of this work is to clarify the shape of the plug flow region in 2D radial flow of a Bingham fluid theoretically and numerically, and to evaluate the discrepancy between the two analytical solutions for theoretical analyses of grouting in rock fractures. More importantly, the grout volume and flowrate evolution as functions of grouting time and relative penetration length are presented, using a correct expression of the plug flow region.The originality of this study lies in the following three aspects:(1)we numerically verified the steady state solution of flowrate and the shape of plug flow region,which has not been done before;(2)we showed the simplest steps to obtain the correct expression of the plug flow region using the force balance equation which implies a radius-independent plug flow region; and (3) we made a direct comparison between the two solutions in the literature to assert the discrepancy.

2. Problem formulation

Fig.1. Schematic of grouting with 2D radial flow.

Fig.1 presents the physical system of grouting in 2D radial flow for an idealized planar fracture represented by a pair of parallel disks. The fracture aperture is 2B, which is basically characterized by hydraulic test in applications.The cement grout is considered as a Bingham fluid injected through a perpendicular borehole with a constant pumping pressure Pg.The radius of the injection borehole is r0.The pressure gradient between the pumping pressure and the in situ groundwater pressure Pwdrives the grout spreading from the borehole. A radial cross-section A-A shows the grout propagating in the fracture. The distance between the edge of the borehole and the grout front represents the grout penetration length I(t), which is a function of the grouting time.

It is assumed that the viscosity of the cement grout (approximately 0.025 Pa s) is higher than that of groundwater (normally around 0.001 Pa s),and that the pressure drop of the groundwater flow is negligible (i.e. the groundwater pressure is constant).Therefore, at each snapshot in the grouting process, the flow of grout can be simplified as steady-state radial flow. It is also hypothesized that the grout is incompressible,the inertial effects are negligible,and the fracture aperture is much smaller than its lateral dimensions (i.e. the pressure gradient across the aperture is negligible by adopting the lubrication approximation).

Given the above assumptions, mass and momentum conservations lead to the following governing equations for 2D radial flow(Bird et al.,1960):

where r is the radius,vris the radial velocity,P is the pressure and τzris the shear stress. For a Bingham fluid,τzris expressed as

where τ0is the yield stress and z is the vertical coordinate along the fracture aperture.The boundary conditions for this system are

where zpis half of the plug flow region;r0and rgare radius of the injection borehole and the grout front,respectively(Fig.1). Eq. (4)represents the known pressure boundary conditions; Eq. (5) denotes the non-slip boundary condition at the fracture walls;and Eq.(6)represents the zero-shear rate condition in the plug flow region.

3. Analytical solutions

Solving Eqs. (1) and (2) with Eqs. (3)-(6) yields pressure and fulid velocities at steady-state (see Appendix A). Since= 0 is assumed,the pressure is a function of r only.The momentum Eq.(2)is integrated over z to obtain

where C is an integration constant determined by introducing the boundary condition at z = zp(the interface between the solid and fulid regions), i.e. Eq. (6), as C =For the fluid part, the shear stress is

For the solid(plug)region(i.e.0 ≤z ＜zp),where the shear rate is 0,the shear stress is explicitly given by the definition of Bingham model, i.e.τzr＜τ0. Determining the plug flow region is an important step to obtain a complete analytical solution for radial flow of Bingham fluids. The simplest and most straightforward way to determine the plug flow region is by using force balance equation.Specifically, the difference between the force at the injection borehole and that at the grout penetration front is equal to the total friction force on the wall surfaces. This can be expressed as

where the shear stress on the wall τwis

The integration of shear stress on the wall gives τ0(rg- r0)+(Pg- Pw)(B - zp). Half of the plug flow region zpcan then be obtained as

Eq.(11) implies that the plug flow region is independent of the radius r,i.e.it is constant along the radial direction.Since the shear stress and pressure are independent of the circumferential angle,the radial cross-section A-A is provided(see Fig.1)to illustrate the force balance in Fig.2.Note that the same radius-independent plug flow region can also be determined by using the energy dissipation equation (El Tani, 2012), see Appendix B.

Quantifying the shear stress is a key in identifying different plug flow regions and therefore providing different analytical solutions in the literature,when following different ways of determining the integration constant C in Eq.(7).One common mistake is assumed when C = 0, i.e.τzr= 0 for z = 0; such a case, however, is applicable only for fluids without yield stress, e.g. Newtonian and power-law fluids. Dai and Bird (1981) essentially adapted this assumption from the solution of 1D slit flow by a similarity or analogy argument. For Bingham fluid flow in a 1D channel, half of the plug flow region is

Fig. 2. Demonstration of shear stress, plug flow region, and force balance in a crosssection.

where ΔP and L are the pressure difference and distance between the injection inlet and the fracture length, respectively(Bird et al.,1960).In the solution for 2D radial flow by Dai and Bird(1981),half of the plug flow region is obtained by replacing(12).The plug flow region thenbecomes a function of the radius r,sinceis a nonlinear function of the radius,according to Eq.(A4).This yields a different result for the flowrate and grouting penetration from the case when C≠0 in Eq. (7). Note that the flowrate obtained from Eq. (12) is radius-dependent, which violates mass conservation. A schematic illustration of radius-independent and radius-dependent plug flow regions for 2D radial flow is shown in Fig. 2. Based on Eqs. (8) and (11) or (12), the flowrate and velocity profiles can be computed (see Appendix A).

4. Numerical verification

In addition to above theoretical analyses,numerical simulations are conducted to verify the radius-independent plug flow region for steady-state radial flow of Bingham fluids and the analytical solution for the flowrate evolution. It is worth mentioning that the numerical simulations are not restricted by the assumptions required for deducing the analytical solution in Section 2.

The general governing equations of mass and momentum conservations for steady-state flow process of incompressible fluids can be expressed as

where u is the fluid velocity vector, ρ is the fluid density, τ is the shear stress, p is the pressure, I is the identity matrix, and g is acceleration due to gravity.

The commercial finite element software COMSOL Multiphysics 5.4 (COMSOL, 2019) is used to solve the governing equations and conduct numerical simulations for comparison with the analytical solution.For the Bingham model,the shear stress is discontinuous due to the yield stress, which causes numerical difficulties in simulations. Papanastasiou (1987) proposed a modified Bingham model(referred to as the Bingham-Papanastasiou model)by adding an exponential term to regularize the yield stress,which is written as

where m is an exponential index;τijis the shear stress tensor;and ˙γijis the shear rate(rate-of-strain) tensor, which is given by

As m becomes close to infinity(m →+∞),the modified model Eq. (15) approaches the Bingham model. In practice, taking m as a large number,e.g.m=100 in this study,is adequate to approximate the Bingham model.

In COMSOL, the modified rheological model, i.e. Eq. (15), can easily be introduced into the non-Newtonian flow solver by using the apparent viscosity, expressed as

where μais the apparent viscosity.A parallel disk’s geometry model with r0= 0.025 m,r = 0.125 m m and aperture 2B=0.001 m was built for numerical simulations. The inlet boundary is set with constant pressure, Pg= 5000 Pa, and the outlet boundary is set with zero pressure,Pw= 0 Pa.The upper and lower disks are set as no-slip walls,i.e.u=0.The viscosity(μ)and yield stress(τ0)of the Bingham fluid is 0.025 Pa s and 5 Pa,respectively.According to the analytical solution,i.e.Eq.(11),the size of plug flow region is 2zp=0.0002 m.

Fig. 3 presents the velocity profiles on different vertical crosssections at different locations along radial direction, i.e. r =0.06 m, 0.07 m, 0.08 m, and 0.09 m. The velocity profiles at locations near the middle of the fracture model were chosen to avoid potential boundary effects(i.e.entries and exits).It shows that the analytical results match perfectly with the numerical results. The plug flow region in the middle of the fracture is independent of the radius r. The numerical results confirmed that the analytical solution and the formula of radius-independent plug flow region are correct.

Fig. 4 presents the comparison of the flowrate against the pressure difference between analytical and numerical results. The flowrates calculated by the analytical solution, i.e. Eq. (A9), based on the radius-independent plug flow region match well with that obtained by numerical simulations (Fig. 4a). Following the mass conservation,the flowrate obtained both by Eq.(A9)and numerical simulation is shown as radius-independent plug flow region.However, the flowrates calculated by the equation of Dai and Bird(1981) are dependent on the fracture radius (Fig. 4b), which is in contradiction to the mass conservation. This comparison result again confirms that the flowrate equation, i.e. Eq. (A9), is correct and that the plug flow region is radius-independent.

5. Real time grouting control

Using the flowrate equation by Dai and Bird(1981),Gustafson et al. (2013) was the first to develop the grout penetration function with relative grouting time. It identified the steering parameters for rock grouting, i.e. pressure, grout properties and fracture aperture, which lay theoretical foundation for the RTGC method (Gustafson and Stille, 2005; Stille, 2015). For clarification, the detailed solution for the grout penetration based on the correct expression of plug flow region is presented in Appendix A.

Fig. 3. Comparison of velocity profiles on different cross-sections along the r-axis,showing a constant plug flow region in the middle of fracture.

Fig.4. Comparison of flowrate between analytical and numerical results:(a)Flowrate against pressure difference, and (b) Flowrate against radius.

In addition to the grout penetration length, the injected grout volume and flowrate, as functions of grouting time, are of great importance in grouting design and practice, since they are often taken as stop criteria for grouting (Gustafson and Stille, 2005;Kobayashi et al., 2008; Rafi and Stille, 2014; Stille, 2015). The injected grout volume Vgcan be obtained through the grout penetration length I (Gustafson et al., 2013), written as

The maximum injection volume Vg，maxarrives when I = Imax,which is

where γ = Imax／r0. The relative volume of injected grout is

where ID= I／Imaxis the relative penetration length. Using the correct formula of relative penetration rate, i.e. Eq. (A15), the evolution of grouting flowrate Qgcan be determined according to its definition, written as

where tDis a dimensionless time,

A simplified empirical formula of relative penetration as function of dimensionless grouting time was proposed by Gustafson and Stille (2005) to approximate the evolution of flowrate, due to the complexity of the original solution of grout penetration developed by Gustafson et al. (2013). However, by using the exact solution of grout penetration, i.e.Eq. (A15),a closed-form solution of flowrate as function of grouting time can explicitly be obtained for the first time without any approximation;the result is shown in Eq. (22). This closed-form solution of the flowrate can be directly applied to controlling and monitoring grouting in practice.

In order to verify the flowrate evolution equation presented in this study, a series of simulations for steady state flow with increasing size of radius (i.e. different penetration lengths at different time steps) is conducted. Due to the fact that a large number of meshes are needed for cases when γ is relatively large,e.g.γ ＞100,and because the small mesh sizes are controlled by the aperture(around 0.1 mm),the presented simulations are limited to cases with small γ values. Nevertheless, the small γ values are adequate to verify the solution of flowrate evolution.

The same fluid properties,i.e.μ=0.025 Pa s and τ0= 5 Pa,the same basic geometry parameters, i.e.r0= 0.025 m,and the same aperture 2B = 0.001 m, are adopted for all the simulations. Two cases with the driving pressure Pg-Pw= 2500 Pa and 25,000 Pa are simulated,so that we have γ = 10 and 100,and Imax= 0.25 m and 2.5 m, respectively. The simulated flowrates for different relative penetration lengths are compared with the analytical solutions, as presented in Fig. 5. It shows that the simulation results match well with the curves of the analytical solution, indicating that the expression of flowrate evolution presented in this study is correct.

Note that the real grout penetration is a transient multiphase(grout-water) flow process. A direct analytical solution is unavailable and direct transient numerical simulation remains difficult due to the discontinuous yield stress and complex transient multiphase flow process. The analytical solutions for grout penetration and flowrate evolution are based on the assumption of quasi-steady state approximation,where the flow at each time step is assumed in a quasi-steady state. Therefore, the transient flow process can be approximated by using the steady-state solution at each time step. Instead of a direct transient simulation, the numerical results presented in Fig.5 are obtained from simulation of the steady-state flow with corresponding penetration length at each time step. The mean velocity at the penetration front is subsequently used to calculate the penetration length for the next time step, following the assumption of quasi-steady state approximation. Examining the validity of such a quasi-steady state approximation remains an open issue; however, it is beyond the scope of the present study.

Fig.5. Comparison of analytical and numerical results for the flowrate evolution with relative penetration length.

In view of the theoretical discrepancy discussed above,we shall refer to the solution by Gustafson et al. (2013) as an additional approximation solution. Since this approximate solution is widely used in the literature,it is important to evaluate its differences from the exact solution. The results of penetration length and flowrate are compared to the exact solutions within the full range of relative penetration length, as shown in Figs. 6 and 7. According to the theoretical analysis, the steering parameter in the functions of relative penetration length, injected grout volume, and flowrate is γ. The relative penetration and flowrate evolution curves for different γ values are calculated and compared in Figs.6 and 7.The adopted parameters (typically used in grouting practice) for the calculation of flowrate are summarized in Table 1.

As shown in Fig. 6a, the grout penetration is increasing but decelerating with elapsed time, regardless of the parameter γ. In contrast, the grout flowrate is decreasing with grouting time(Fig.7a)as well as the relative penetration length(Fig.7b).With the increase of γ, the penetration becomes relatively slower and the corresponding flowrate becomes relatively larger. By using the penetration and flowrate functions, the grouting process can be fully defined for design, monitoring and quality assurance in grouting applications.

Fig.6. Comparison between exact and approximate analytical solutions of the relative penetration length with (a) dimensionless grouting time and (b) log-log plot of the complementary relative penetration length (1 - ID), highlighting the discrepancy in relatively longer grouting time (i.e. tD ＜1).

Fig. 7. Comparison between exact and approximate analytical solutions of the grout flowrate evolution with (a) dimensionless grouting time and (b) relative penetration length.

Table 1Grouting parameters for the illustration example.

The results of the approximate solution are in good agreement with those of the exact solution within a shorter penetration distance(i.e.ID＜0.5)and a lower dimensionless time(i.e.tD＜1),but gradually deviate with the increases of penetration distance and grouting time,especially for higher values of γ.Particularly,a loglog plot of the complementary relative penetration length (1-ID)is presented in Fig.6b,to highlight the discrepancy in the relatively longer grouting time (i.e. tD＞1). Such deviations demonstrated that the two solutions are not equal, since they are based on different expressions of the plug flow region. The reason for that the differences are relatively small is because the values of the termsare relatively close. Therefore,the differences of the plug flow region in fractures with such small apertures(0.1 mm)are relatively small.In addition,since

6. Conclusions

Analytical solutions for steady-state,2D radial flow of a Bingham fluid is important in grouting applications, which constitutes a basis for design tools using quasi-steady state assumption. The design tools and the quasi-steady state assumption are based on the work of Gustafson et al.(2013)which in turn builds on the work of Dai and Bird(1981)that the plug flow region for 2D radial flow of a Bingham fluid is analogous or similar to the 1D channel case. In the present study, we revisited the analytical solutions for steadystate, 2D radial flow of a Bingham fluid and the following conclusions can be drawn:

(1) For 2D radial flow of Bingham fluid, the plug flow region is independent of the radius at any time.The constant C in Eq.(7)is critical for solving Eqs.(1)-(9)as it determines the plug flow region zp.It shows that the simplest steps to obtain the plug flow region zpis from the bulk force balance, as illustrated in Fig. 2; the bulk energy balance leads to the same value of zp(El Tani,2012).It also suggests that neglecting the constant C in Eq.(7)or adapting the plug flow region for 2D radial flow from the 1D channel case (Dai and Bird, 1981)yields zpthat is dependent on the radius which is unphysical and violates fluid mass balance.

(2) Steady-state numerical simulations for radial flow of Bingham fluids in a parallel fracture verify the theoretical analyses, showing a good agreement. In the present numerical simulations,the Bingham-Papanastasiou model is adopted to regularize the discontinuous yield stress, which shows adequate precision. Due to computational capacity, current numerical simulations are limited to relatively small values of γ, which constrains its applicability in rock grouting practice.

(3) The quasi-steady state solution of Gustafson et al. (2013),which is based on the conjecture that C equals 0 in Eq.(7)and zpis adapted from the 1D channel case(Dai and Bird,1981),is in good agreement with the exact solution,i.e.Eqs.(A15)and(22), within a shorter relative distance (i.e. ID＜0.5), but gradually deviates with the increases of penetration length and grouting time,especially for higher values of γ(i.e. γ =1000).From a practical point of view,this solution seems to be sufficiently accurate.

(4) Using the correct expression of the plug flow region, the grout penetration length, injected volume and flowrate evolution with elapsed grouting time under constant pumping pressure are presented.The closed-form solution of flowrate evolution as a function of grouting time is presented in this study for the first time. These evolution functions provide fundamental inputs for rock grouting design and monitoring.

The analytical solutions discussed in this work are at best applicable for scoping calculations since they are based on several restrictive assumptions: idealized geometry of a fracture (i.e.smooth parallel disks), the simplest non-Newtonian rheological model (i.e. the Bingham model) for grouts, and the condition of constant injection pressure.Rock fractures consist of rough surfaces that aggregate into complex networks with variable apertures.Moreover, grouting is fundamentally a two-phase flow process in which non-Newtonian fluids(grouts)displace groundwater in rock fractures (Zou et al., 2018). Implementation of the analytical solutions discussed in this work for operation and design of real systems should therefore always be expected to result in some level of uncertainty. Developing numerical tools for more realistic rheological properties of grouts in realistic structures of rock fractures and associated networks,as well as for the two-phase flow process,will be an important step towards the development of quantitative tools for design and operation of rock grouting applications (e.g.Zou et al., 2018, 2019).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The funding for this study is provided by the BeFo Rock Engineering Research Foundation (Grant No. 392). We would like to thank Professor H?kan Stille, Dr. Lars H?ssler, Professor Johan Claesson, and Dr. Johan Funehag for their helpful discussion.

Appendices A and B. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2019.12.021.