Multiphase flow coupling analysis of shield synchronous grouting filling law with time-dependent viscosity

The expansion of urban transportation systems frequently clashes with above-ground structures. Underground tunnels offer an effective solution to this challenge. Shield tunneling, as a mechanized construction technique, can mitigate disruptions to surface transportation and social activities in densely populated urban areas. During shield tunneling construction, an annular gap exists between the detached segmental lining and the soil,¹ During shield tunneling construction, an annular gap exists between the detached segmental lining and the soil. Synchronous grouting in shield tunneling significantly impacts the success of the construction,^2,3 playing crucial roles in controlling strata deformation,^4–6 enhancing tunnel impermeability,^7,8 and optimizing the stress state of segments.^9,10

The diffusion process of slurry during shield construction can be divided into filling diffusion,¹¹ infiltration diffusion,^12,13 compaction diffusion,¹⁴ and other stages. The main part of the synchronous grouting slurry is used to fill the diffusion. Infiltration and compaction diffusion primarily develop after the basic completion of filling diffusion. Additionally, the force exerted by the filling slurry on the segment is intimately related to the segment’s buoyancy. Thus, research on filling grouting can offer a more scientifically grounded rationale for analyzing the interaction between grouting and soil, as well as the stress state of the segments during subsequent stages of slurry infiltration and compaction diffusion.

In recent years, scholars have employed computational fluid dynamics, the principle of equilibrium, and slurry constitutive models to theoretically derive the slurry filling diffusion model, analyzing the pattern of slurry pressure distribution and its determinants.¹⁰ Bai et al.¹⁵ deduced the slurry pressure distribution model in the cross section of the shield tail gap through the Newtonian fluid model. Based on Bingham fluid, Liu¹⁶ derived the circumferential distribution model of slurry filling pressure. Based on simplified assumptions, Hu et al.¹⁷ analyzed the spatial mechanics of the circumferential fluid micro-element in the shield tail and established the circumferential filling diffusion model of the slurry with two-dimensional flow. Li et al.¹⁸ further supplemented the longitudinal filling model of shield synchronous grouting, yet overlooked the intricate coupling between longitudinal and circumferential diffusion of slurry. Chen¹⁹ employed the circumferential independent model and the overall filling diffusion model to calculate and analyze the slurry pressure distribution model, underscoring the superiority of the latter through the L5 Nanjing rail transit project. Li and Liu^20,21 devised a theoretical model for the longitudinal and circumferential overall diffusion in shield tail synchronous grouting, introducing the computational fluid dynamics (CFD) method for simulating the overall filling and diffusion process. However, the crucial impact of the slurry’s time-dependent viscosity was neglected.

As derived from the aforementioned research acknowledging the integral coupling of the slurry in both longitudinal and circumferential directions and its evolving viscosity²² is more conducive to the accurate description of the filling process. However, due to the difficulty of monitoring the flow process of synchronous grouting filling in shield tunnel construction,²³ previous studies have not been able to explore the diffusion mechanism of synchronous grouting filling through accurate description. Present models harbor extensive assumptions and oversimplifications, which regard the filling process as a two-dimensional flow process (circumferential filling and longitudinal filling). However, these studies ignore the actual motion law of the slurry in the shield tail gap. In three-dimensional space–time, the synchronous slurry moves in the shield tail gap, and the velocity vector at any time has both circumferential and longitudinal components. A holistic interpretation of the synchronous slurry’s diffusion demands a three-dimensional model that contemplates both the integral coupling of the slurry's motion and its time-dependent viscosity. As an important research method, the numerical analysis method has also been widely used in the study of grouting behind the shield wall. Yet, predominant studies hone in on strata deformation²⁴ and the structural stresses of the segment.²⁵ The solid element with elastic modulus changing with time is used to simulate the shield tail grouting layer, focusing on the influence of synchronous grouting on the strata environment and lining structure, and the filling and diffusion law of slurry in the shield tail gap is not analyzed. Given the feasibility of finite element simulation in calculating the slurry diffusion law in rock fractures,²⁶ this research uses the finite element software COMSOL Multiphysics to simulate the shield synchronous grouting. The model considers the longitudinal and circumferential overall coupling of the slurry, distinguishes the significant differences in the viscosity development stages of the aged and fresh slurry, and considers it as a multiphase flow transient problem for analysis.

Combined with a shield project in Hangzhou, China, this research systematically analyzes the influence of grouting construction parameters (grouting pressure, speed) and slurry physical properties (slurry viscosity, density) on slurry diffusion motion characteristics and buoyancy of segments, which has certain scientific research significance and engineering application value.

As illustrated in Fig. 1, the slurry used in a Hangzhou shield construction project was tested for density and viscosity. The grouting materials, sourced directly from the construction site, had their compositional design and mass proportion outlined in Table I. The slurry was thoroughly mixed using a JJ-5 slurry mixer. The slurry’s density was measured at 1823 kg/m³ using a mortar density gauge (JGJ/T70-2009). Viscosity tests on the prepared slurry were conducted using an American Brookfield DV3T rheometer (outfitted with LV1-5 rotors) at a temperature of 20 ± 1 °C. From the experimental data, a fitting provided the time-dependent expression for slurry viscosity as μ(t) = 10.97e^0.0147t.

Owing to the slurry’s time-dependent viscosity, when progressing longitudinally to a new section, the fresh slurry introduced via the synchronous grouting pipe will replace the previously existing slurry within the shield tail’s annular gap. The fresh and aged slurries are at distinct stages of viscosity development, exhibiting notable viscosity differences, so they are considered as two phases. According to preliminary calculations, the Reynolds number (Re) of the two-phase slurry is relatively low, and the inertial force cannot be neglected, while the viscous force is limited to the boundary layer, shear layer, and wake region. The fluid can be regarded as a laminar flow with regular and smooth flow patterns. Consequently, this study integrates the laminar flow interface with the two-phase flow interface, utilizing the level set method for tracking the interface between these phases. During the numerical simulation of slurry fill and diffusion, the distribution of the two phases across the entire computational realm is depicted using the volume fraction approach.

Although two-dimensional numerical simulations are widely used for studying engineering issues related to tunnel excavation,^29,30 recent research has indicated that three-dimensional numerical modeling is irreplaceable for the computational analysis of understanding grouting.^31–34 In this research, a three-dimensional grout-filling model for shield tunneling is constructed. Due to the computational intensity of three-dimensional finite element models, the computational cost is reduced by utilizing the symmetry of the computational domain and modeling only half of the model’s geometry.

The numerical model is based on an 8-hole grouting shield tunnel project in Hangzhou, modeled at a 1:1 full scale. The annular shield tail gap is enclosed by the shield tail boundary, the outer side of the tunnel segment, the strata interface, and the artificial boundary at the tail. The inner diameter of the shield tail gap is 14.5 m (equivalent to the outer diameter of the tunnel segment), the outer diameter is 15.01 m, and the thickness of the shield tail gap is 25.5 cm. The longitudinal width of the shield tail gap is 2.0 m (equivalent to the longitudinal length of one tunnel segment). The centers of the grouting holes align with the gap's central axis. As depicted in Fig. 2, these holes are evenly distributed at 40° intervals and are sequentially numbered in a clockwise direction (①-④). Each grouting hole has an inner diameter of 8 cm.

The grouting inlet is configured as a constant pressure inlet boundary, while the artificial boundary at the tail adopts a constant velocity outlet boundary. The chosen constant velocity corresponds to the average speed of shield excavation, which is 0.6 mm/s, thus simulating the shield excavation process and the physical process of gradually moving the shield tail boundary away from the solidified grout layer. All other boundaries are defined with no-slip conditions. Initially, the shield tail gap is occupied by the previously injected grout from the advancement of the prior ring. Concurrently, the grouting tubes are prepared for the injection of fresh grout for the upcoming ring. In the model’s initial computational state, the fill phase within the shield tail gap is designated as the aged grout, injected synchronously during the progression of the previous ring. Meanwhile, the initial fill phase in the grouting pipes corresponds to the fresh grout synchronously injected for the upcoming ring advancement. Consequently, the initial interface of the two phases is situated at the interface between the grout pipes and the cavity awaiting filling. Under a constant grouting inlet pressure, the fresh grout steadily injects and occupies the shield tail gap. As the tunneling shield advances to form a new annular gap, the aged grout is successively expelled from the computational domain of the single annular gap. Once the tunneling shield concludes its longitudinal progression for the new ring distance, the fresh grout concurrently achieves the filling process within the shield tail gap.

To enhance the analysis of computational results, six observational lines, as depicted in Fig. 3, were established. Observational lines 1–5 are all located on the axis plane of the shield tail gap, while observational line 6 is positioned on the outer side of the tunnel segment. Lines 1–4 align with the centers of grouting holes 1–4 and are parallel to the lateral surface. Observational lines 5 and 6 are semi-circular and parallel to the front wall of the shield tail gap, with a distance of 100 mm from it. Lines 1–4 will delineate both the pressure and velocity field distributions in the direction of grout injection at the inlet. Observations from line 5 provide insights into the circumferential velocity field distribution within the shield tail gap, while the observations from observational line 6 reveal the circumferential pressure field distribution on the outer side of the tunnel segment.

As shown in Table II, to investigate the effects of grouting pressure, time-dependent slurry viscosity, and slurry density on the grouting filling of the shield tail void, seven computational cases were established for comparative analysis. Case 1 was based on the actual parameters from a specific tunnel project in Hangzhou and served as the control case.

To analyze the impact of grouting pressure on the filling process, Cases 2 and 3 utilized different grouting pressure combinations relative to Case 1. Case 2 had a grouting pressure of 50 kPa lower than Case 1, while Case 3 had a grouting pressure of 50 kPa higher.

To assess the influence of slurry viscosity’s time-dependency on grouting filling, Cases 4 and 5 used different time-dependent viscosity expressions compared to Case 1. Case 4 employed a rapid-setting slurry viscosity expression μ₂(t) = 32.27e^0.017t mPa s. Case 5 applied a retarding slurry viscosity expression μ₃(t) = 8 mPa s.³⁵ Given that each ring advancement in shield tunneling averages a 2-h interval, the viscosity value for the aged slurry is determined at time T = t + 120 min.

To analyze the effect of slurry density on the grouting filling, Cases 6 and 7 used different slurry densities relative to Case 1. The slurry density in Case 6 was set 300 kg/m³ lower than in Case 1. The slurry density in Case 7 was set 300 kg/m³ higher than in Case 1.³⁵ 

Combining the parameter settings of all aforementioned cases, the influence of various grouting factors on the slurry filling process is systematically analyzed by the research method of controlling variables.

Under varying grouting conditions, the slurry diffusion showcases some similar patterns. As depicted in Fig. 4, typical Case 1 is selected to analyze the diffusion pattern of slurry. The side views of the slurry volume fraction are captured at 9 representative moments (T = 0, 10, 20, 30, 40, 100, 150, 300, 400 s) to depict the entire diffusion progression of the slurry. The contour plot illustrates the spatial distribution of volume fractions for both fresh and aged slurries. The boundary, defined where the volume fraction of each phase is 0.5, represents the two-phase interface: the blue zone indicates the aged slurry within the shield tail gap, while the red zone signifies the freshly injected slurry.

The entire slurry filling and diffusion process can be divided into two stages: the independent filling stage of each hole and the combined filling stage of multiple holes. During the initial grouting phase, each of the four grout holes, driven by the grouting pressure, progressively fills the gap in the shield tail with a diffusing fan-like pattern of the slurry. Owing to the slurry’s self-weight, it exhibits a downward movement tendency while diffusing and filling. Given that slurries injected through holes 2 and 3 diffuse around surfaces with slopes of 1.732 and 5.671, respectively, their self-weight effect is pronounced. Using an isosurface with a two-phase volume fraction of 0.5 as the diffusion boundary, at T = 30 s, the slurry from grouting holes 2 and 3 merges first with that from hole 4. Conversely, the slurry from hole 1, spreading around a surface with a mere slope of 0.364, exhibits a more gradual downward movement. By T = 40 s, all slurries from all four grouting holes transition to the combined filling phase. Notably, at T = 400 s, an incompletely filled region remains on the upper part of the shield tail gap, suggesting this area is prone to underfilling. Thus, during the actual project construction, it is essential to intensify the monitoring of the filling efficacy in this region. Should underfilling be detected, prompt secondary grouting is crucial to ensure the shield tail gap is adequately filled. This approach mitigates the engineering risks associated with disturbances in the surrounding strata and excessive surface subsidence due to uneven filling.

Figure 5 depicts the velocity distribution of the slurry along the semi-circular observational Line 5. Before T = 40 s, noticeable velocity peaks can be observed at four angles (20°, 60°, 100°, 140°) corresponding to the independent filling stage of each hole. Within the four fan-shaped grouting diffusion zones, the slurry exhibits a relatively high flow rate. Outside these zones, the displaced aged slurry recedes slowly by the pushing of the newly injected slurry, resulting in a relatively slower flow rate. The peak velocity increases with a wider angle between the tangential line of the hole’s position and the horizontal axis, attributed mainly to the slurry’s weight. As illustrated in Fig. 5, with the progression of time and the evolution of grouting into the combined filling stage of multiple holes, the pronounced regions of peak slurry velocity begin to diminish, initiating from the bottom and gradually extending upward. Concurrently, velocity across all areas diminishes, moving toward a uniform distribution. This characteristic evolution of the circumferential velocity field results from the diffusion pattern of the slurry filling, which progress from the bottom to the top. In the vicinity of grouting hole 1, the velocity exhibits a prominent peak due to a substantial pressure difference between the grouting and surrounding pressures. A more comprehensive analysis of this pressure field is elaborated in the sections that follow.

As illustrated in Fig. 6, the pressure distribution along line 6 showcases an eccentric quasi-circular pattern, with a gravity-dominated characteristic, being less pronounced at the top and more accentuated at the bottom. As the grouting proceeds, the pressure surrounding the segment consistently rises, albeit at a diminishing rate. Additionally, pressure augmentations are particularly evident in the lower regions of the shield tail gap. Notably, the pressure at the bottom of the shield gap (θ = 180°) surges from 309.94 to 376.29 kPa, marking an increase of 66.35 kPa. Conversely, the pressure at the segment’s top (θ = 0°) remains stable, ranging between 133 and 147 kPa, with no evident fluctuations.

Given the symmetry of the computational domain, the horizontal components resulting from the circumferential pressure cancel each other out. Thus, by solely integrating the vertical forces of the radial pressure, one can ascertain the cumulative buoyancy exerted on the tunnel lining at different times.

The relationship between the buoyancy evolution curve in Fig. 7 and the contour plot suggests a direct link between the buoyancy acting on the tunnel lining and the development of filling. With the rapid escalation of the filling rate, the buoyancy imparted on the tunnel lining intensifies swiftly. Upon complete filling of the shield tail gap, the buoyancy per unit length of the lining reaches its peak value of 2960 kN/m.

As illustrated in Fig. 8, under different grouting pressure conditions, there are significant differences in the numerical value of the circumferential distribution curve family of the outer pressure field of the segment, but there are common characteristics in the variation law with time. Specifically, with the increase in grouting pressure, the circumferential pressure of the segment increases, and only the pressure at the top of the segment (θ = 0°) changes little with time. At the bottom of the segment (θ = 180°), the pressure increases most obviously with time, and Cases 1, 2, and 3 increase by 68.44, 68.05, and 68.424 kPa, respectively. It can be seen that the slurry is seriously deposited at the bottom of the shield tail gap.

In this study, simulations were conducted on conventional, rapid-setting, and retarded slurries under the same operational conditions (1, 4, 5). It was observed that the rate of development of slurry viscosity plays a significant role in the effectiveness of slurry filling. Analyzing from a diffusion pattern perspective, slower development of slurry viscosity results in faster diffusion of slurry filling. As shown in Fig. 9, the retarded slurry can fill the shield tail void within 200 s. In contrast, the conventional slurry struggles to completely fill the top part of the shield tail gap before 200 s, and the rapid-setting slurry leaves a significant unfilled area at the bottom of the shield tail gap within the same timeframe. Moreover, the rapid-setting slurry displays a distinct characteristic of accelerating first and decelerating later during the diffusion process. As depicted in Fig. 10, the synchronous grouting filling effect in the 420th ring of the project site was scanned using the LTD-2100 geological radar. The radar detected areas at the top of the shield tail gap that were challenging to fully fill. This observation aligns with the calculated results of Case 1, which employed onsite parameter settings. Consequently, the non-destructive testing results from the geological radar confirm the model's capability to effectively predict the pattern of filling and diffusion. This provides insightful references for adjusting subsequent construction parameters and positioning for secondary grouting.

From the perspective of buoyancy development (Fig. 11), slurries with faster viscosity enhancement provide a slower rate of buoyancy increase. Since the rapid-setting slurry does not achieve complete filling even at 800 s, its buoyancy never reaches the maximum value. Corresponding to the previous perspectives, the slurry velocity is faster when the viscosity development is slower. As depicted in Fig. 12, the slurry velocity of the retarded type at T = 20 s on line 6 is significantly higher than the other two types. A higher slurry viscosity indicates greater internal shear resistance during flow, thus consuming more slurry kinetic and potential energy. All other conditions being equal, a rapid increase in slurry viscosity is detrimental to slurry diffusion.

In this study, high-density, low-density, and conventional slurries were simulated. As depicted in Fig. 13, the findings suggest that as the slurry density increases, both the gradient of buoyancy growth of the lining during the filling process and the maximum limit of buoyancy increase. Moreover, under identical grouting conditions, a higher density of the slurry is beneficial for the displacement of the aged slurry by the fresh one, which can be attributed to the greater inertia of the denser slurry. Taking the velocity distribution on observational line 6 at T = 20 s as an example (Fig. 14), under the same spatiotemporal conditions, a higher slurry density corresponds to a faster diffusion rate.

This study establishes a three-dimensional filling diffusion model for synchronous grouting based on multiphase flow coupling with consideration of the time-dependent viscosity of the grout, thereby unveiling the kinetic laws governing the grout diffusion process in large-diameter shield tunneling. A systematic analysis of the time-dependent properties of grout viscosity, grouting pressure, and grout density on diffusion patterns and segment buoyancy has been conducted, offering theoretical insights for the optimization of synchronous grouting parameters, the design of secondary grouting strategies, and the mitigation of surface settlement. From this research, the main conclusions can be drawn as follows:

1. Incorporating multiphase flow coupling and considering the time-varying viscosity of the slurry, this research actualized a three-dimensional filling diffusion simulation for large-diameter shield tunnels during synchronous grouting. Based on the spatiotemporal evolution of the slurry filling pattern, the diffusion process is divided into two phases: the independent filling stage of each hole and the combined filling stage of multiple holes.

2. It is observed that as the density of the slurry increases, there is a corresponding enhancement in the buoyancy growth gradient of the segment during the filling process, leading to a higher peak buoyancy. Concurrently, with an increase in slurry viscosity, there is an elevation in its internal shear resistance, necessitating greater consumption of both kinetic and potential energy of the slurry. Utilizing a low-density, rapid-setting slurry for tail void filling under high-pressure grouting is advantageous, ensuring a commendable filling rate while mitigating the adverse effects of segment buoyancy.

3. Due to the influence of gravity, there is an inherent downward movement tendency in the slurry. In steeper gradients, the impact of gravity becomes increasingly pronounced, leading to accelerated diffusion rates. For conventional slurries, the upper gap in the shield tail is identified as an area prone to incomplete filling, making it a focal point for secondary grouting in shield grouting projects.

4. Higher slurry viscosities entail increased internal shear resistance during flow. All other conditions being constant, a rapid increase in slurry viscosity is identified as a detrimental factor to slurry diffusion.

5. The simulation operates on the assumption that synchronous slurry diffusion is predominantly driven by filling diffusion, excluding considerations of slurry infiltration diffusion and compaction diffusion into the surrounding soil. Consequently, the simulated outcomes are most pertinent to shield tunneling projects where the strata maintain excellent self-stability, ideal shield tail gap formation, and low permeability.

This paper was supported by the China Postdoctoral Science Foundation (Grant No. 2022M722956).

The authors have no conflicts to disclose.

G.L. and H.C. contributed equally to this work.

Guanglun Li: Resources (equal); Supervision (equal); Writing – review & editing (equal). Hongtao Cao: Conceptualization (equal); Data curation (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Jian Wu: Conceptualization (equal); Supervision (equal). Bo Wang: Investigation (supporting); Writing – review & editing (supporting). Xuetao Zhou: Investigation (supporting); Writing – review & editing (supporting). Wei Zhao: Investigation (supporting); Writing – review & editing (supporting). Xiaolan Cai: Investigation (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.