3D Modelling for Mechanized Tunnelling in Soft Ground-Influence of the Constitutive Model

The main purpose of this study was to provide a three-dimensional numerical model which would allow to show the influence of the constitutive model of the ground on the tunnel lining behaviour and displacement field surrounding the tunnel. Most of the processes that occur during mechanized excavation have been simulated in this model. Besides the well-known elastic linear perfectly plastic model with Mohr Coulomb failure criteria (MC), the Cap-Yield model (CYsoil), which is a strain hardening constitutive model that takes into consideration the hyperbolic behaviour of the soil has also been adopted. The results have shown a considerable effect of the constitutive model on the tunnel behavior and the ground displacement field. Generally, the CYsoil model produces higher structural force and ground settlement values than those predicted by the MC model.


INTRODUCTION
Rapid progress in the development of user friendly computer codes and thelimitations of analytical methods have led to an increase in the use of numerical methods for the design of tunnel lining. A summary of the finite element models used for tunnelling analyses prior to 2000 was given by Negro and Queiroz (2000) and by Farias et al. (2004). They showed, after an examination of more than 65 recently published papers, that the Finite Element Method (FEM) is the most popular approach, accounting for 96% of the published cases, while the remaining ones used the Finite Differences Method (FDM) or others. They also noted that 92% of the published analyses were still performed in two-dimensions (2D), under the hypothesis of plane strain conditions. Another relevant fact that emerged was that most analyses still use simple models, mostly elastic-perfectly plastic.
The purpose of a numerical mechanized Tunnelling (TBM) model is to take into consideration the large number of processes that take place during tunnel excavation, such as the applied face pressure, shield overcutting, shield conicity, the annular void behind the shield, grout injection in this void and its consolidation process, segmental concrete lining. In order to conduct a rigorous analysis, a three-Dimensional (3D) model should be used to take into account all these processes.
Tunnel structure behaviour is a complex phenomenon in which the behaviour of the surrounding ground is one of the main aspects of a tunnel excavation that should be taken into account (Oreste, 2007). Consequently, a realistic constitutive model is crucial to estimate the structural force induced in a tunnel structure and ground displacement.

AJAS
The influence of the constitutive model of the soil on the ground displacement has been studied by many authors in the past. However, less works have been devoted to its effect on the tunnel lining behaviour. Generally, simple constitutive models lead to a shallower and wider settlement trough than the one observed experimentally (Bolton et al., 1994;Addenbrooke et al., 1997;Masin and Herle, 2005;Hejazi et al., 2008;Rodriguez, 2008;Lambrughi et al., 2012). Their studies indicate that in order to take into account some of the fundamental aspects of soil behavior, such as the variation inthe stiffness modulus, according to the stress state and the difference in the modulus during the unloading and loading phases, it is at least necessary to use an elasto-plastic model with strain hardening.
The main purpose of this study was to provide a 3D model which would allowto showthe influence of the soil constitutive model on the tunnel lining behaviour and displacement of the ground surrounding the tunnel. Most of the main elements of a mechanized excavation are simulated in this model: the conical geometry of the shield, the face pressure, circumferential pressure acting on the soil surface in the working chamber behind the tunnel face, circumferential pressure due to the migration of the grout acting on the soil surface and at the shield tail, grouting pressure acting simultaneously on the soil surface and tunnel structure after the shield tail, progressive hardening of the grout, the jacking force and the segmentation of the tunnel lining in particular. The Bologna-Florence railway high speed line project has been adopted in this study as a reference case.

The Bologna-Florence Railway Line Project
The Bologna-Florence railway line project is a part of the Italian high speed railway network. The project involves the excavation of two tunnels in Bologna with a space distance of 15 m between the two tunnel centres. The tunnel has an external excavation diameter of 9.4m and an internal diameter of 8.3 m for a useful section of 46 m 2 .The tunnels have been excavated to a depth of between 15m and 25 m below the ground surface. The two tunnels were excavated through two main formations: alluvial deposits of the late Pleistocene -Plioceneera, which ismostly alluvial deposits from the SavenaRiver with deposits of clay and sandy soil (clayey sands and Pliocene clays). Some typical parameters of the ring 582 in the first tunnel, which have been adopted in this study as a reference case, are summarized in Table 1 and 2 (Croce, 2011).
The lining of the tunnels is composed of precast segments made of reinforced concrete. Each 1.5 m long circular ring consists of 6 conical regularly shaped blocks and a small sized key block. Each precast concrete ring has an extrados diameter of 9.1 m and is 0.4 m thick. The mechanical characteristics of the tunnel lining are shown in Table 1.

MC Model
Even though the stiffness remains constant during the loading and unloading phases, the MC model is able to describe the basic features of a frictional material with just a few parameters. For this reason it is commonly used in geotechnical engineering. Five parameters are necessary to model the linear-elastic perfectly plastic behaviour: the elastic parameters, E (Young's modulus) and ν (Poisson's ratio) and the plastic parameters, φ (friction angle), c (cohesion) and ψ (dilatancy angle).

CYsoil Model
The CYsoil model is a strain-hardening constitutive model that is characterized by a frictional Mohr-Coulomb shear envelope (zero cohesion) and an elliptic volumetric cap in the (p', q) plane. Apart from the caphardening law and the compaction/dilation law, which allow the volumetric power law behaviour observed in isotropic compaction tests and the irrecoverable volumetric strain that occurs as a result of soil shearing to be captured, the friction hardening law in the CYsoil model offers the possibility of alternatively expressing the hyperbolic behaviour. In the CYsoil model, the stiffness is adopted as a function of the effective confinement and it leads to a higher value for unloadingreloading stiffness. For the CYsoil model, if friction hardening behavior is adopted, the input parameters are: • Failure ratioR f which is a constant and smaller than 1 (0.9 in most cases) • Ultimate friction angle φ f • Calibration factor β

Parameter Calibration
Numerical models of drained triaxial tests corresponding to the different depths (half of the overburden, top of the tunnel and bottom of the tunnel) of the reference case have been simulated in order to calibrate the parameters of the CYsoil model. The calibration results are presented in Table 2.   The model has been developed using the FLAC 3D software (ICG, 2009) which is based on the Generalised Finite Difference Method. The analyses have been performed using small strain calculations.

AJAS
The tunnel construction process is modelled using a step-by-step approach, as this has been done by other authors (Jenck and Dias, 2003;Melis et al., 2002;Swoboda et al., 2004;Kasper and Meschke, 2004;Mroueh and Shahrour, 2008;Boubou, 2010;Mollon et al., 2013;Oreste, 2003;Oreste, 2013). Accordingly, each excavation step corresponds to an advancement of the tunnel face of 1.5 m, which is equal to the width of a lining ring.

Simulation of Different Phases of the Mechanized Tunnelling Process
A 3D numerical model has been developed which allows the advancement of the TBM in the soft ground to be simulated, taking into account the components and procedures that can occur in actual tunnel excavation as much as possible.
In general, a tunnelling process consists of three main phases (Bernat and Cambou, 1998;Jenck and Dias, 2004;Boubou, 2010;Dias and Kastner, 2013): • Excavating the ground at the tunnel face and simultaneously applying a confinement to ensure tunnel face stability • Installing the tunnel lining, applying the jacking force and injecting the grout behind the segments in order to fill the voids created at the shield tail • The TBM continues to advance and the ground begins to become stabilized, which is expressed by a consolidation phase A schematic view of the present model is provided in Fig. 2.
In this study, the face pressure has been modelled by applying a pressure distribution to the excavation face (Melis et al., 2002;Kasper and Meschke, 2004;Mroueh and Shahrour, 2008;Lambrughi et al., 2012;Mollon et al., 2013), using a trapezoidal profile in order to account of the slurry density. The unit weight of the slurry is taken equal to 11 kN/m 3 . The method to choose the value of the face pressure is based on the vertical stress σ v in the ground mass. The value of the average face pressure (applied at the tunnel axis) is thus generally set equal to the horizontal ground pressure (Mollon et al., 2013) Equation (1): where, K 0 is the lateral earth pressure coefficient at rest and σ v is the vertical ground overburden pressure at the tunnel axis. The applied pressure is taken equal to σ t = 150 kPa for all of the presented numerical analyses. Due to a slight overcutting, a possible slurry migration could occur over a short distance behind the cutting wheel (Dias et al., 2011;Mollon et al., 2013). Therefore, in addition to the pressure acting on the tunnel face, a pressure, caused by the slurry solution, has also been applied to the cylindrical surface just behind the tunnel face. In this study, the cylindrical surface corresponds to the excavation chamber situated between the excavated face and the shield. This slurry is simulated by means of a uniform pressure diagram over a length of 1.5 m and is followed by a reduced pressure acting over a 1.5 m length with a triangular distribution (Fig. 2).
For the sake of simplicity, the "fictive" shield introduced by (Dias et al., 2011;Jenck and Dias, 2004;Simic, 2006;Mollon et al., 2013), that is used in this study has not been simulated considering volume elements. Instead, it is replaced by a condition on the soil displacements which states that the soil cannot ''penetrate'' into the shield (even though this shield does not exist in the model). Thus, at each excavation stage, the position of the envelope of the ''fictitious'' shield is calculated and each node of the soil located on the contour of this virtual shield is artificially blocked when this point began to be in contact with the "fictive" shield. The geometrical parameters of the shield are presented in Fig. 2.
During mechanized tunnelling, the jacking force is one of the primary segment loads in the construction stage (JSCE, 1996;ITA, 2000). The distribution of the jacking force is uneven over the tunnel height. The jacking force is higher at the bottom of the tunnel than at the top. This is caused by the moments the jacks have to exert (Rijke, 2006). The distribution of the jacking force used in this study has been assumed to be linear over the height of the tunnel. In the present model, the jacking forces were simulated by concentrated forces that act Science Publications AJAS directly on the nodes at the edge of the segments. These jacking forces were set on each segment considering three plates located at 1/6, 1/2 and 5/6 of the segment length. Atotal jacking forces of about 40 MN was adopted in the present model on the basis of the theoretical method proposed by Rijke (2006).
In general, after being injected into the void behind the shield tail, the grouting action is modelled in two phases: (1) the liquid state (state 1) represented by a certain pressure acting on the ground surface and on the tunnel lining; (2) the solid state (state 2) (Melis et al., 2002;Mollon, 2010;Mollon et al., 2013). The distributional radial pressure used to simulate this kind of pressure in this study has been adopted by other authors (Melis et al., 2002;Kasper and Meschke, 2004;Mollon, 2010;Mollon et al., 2013). The grouting pressure applied to the tail void is generally set to (Mollon, 2010;Mollon et al., 2013) Equation (2): inj ν σ = 1.2σ ( 2) where, σ v is the ground overburden pressure at the tunnel crown. Aσ inj value of 180 kPa at the tunnel crown has been adopted in the present model. The grout was simulated by adopting a uniform pressure which was applied to both the cylindrical surface of the excavated ground and external surface of the tunnel lining. As for the face pressure, the annular void between the outside surface of the shield and the excavated ground made the migration of some grout towards the shield possible. This migration was simulated by means of a triangular pressure over a certain length (Mollon, 2010;Mollon et al., 2013). In this study, the length of the grouting pressure acting on the lining and the migration length of the grouting pressure on the annulus void were arbitrarily chosen equal to the length of one ring (1.5 m). Beyond this length, the grout was assumed to harden and was simulated by means of volume elements with perfect elastic behaviour and with the elastic characteristics E grout = 10 MPa and σ grout = 0.22 (Mollon, 2010;Mollon et al., 2013) (Table 1).
In the present model, the tunnel segments have been modelled using a linear-elastic embedded liner element. This element was used to model a thin liner (based on the classical Kirchhoff plate theory) for which both normal-directed compressive/tensile interaction and shear-directed frictional interaction with the host medium occur (ICG, 2009;Do et al., 2013).
The segment joints have been simulated using double node connections. The stiffness characteristics of the joint connection have been represented by a set composed of a rotational spring (K θ ), an axial spring (K A ) and a radial spring (K R ) Table 3 (Do et al., 2013).
On the basis of empirical data (Cavalaro and Aguado, 2012), the behaviour of axial springs has been represented by means of a linear relation using a constant coefficient. The radial stiffness and rotational stiffness of a segment joint were instead modelled by means of a bilinear relation that is characterized by a stiffness factor and maximum bearing capacity (Do et al., 2013).
In the same way, as for the segment joint, the ring joint has also been simulated using double connections. The rigidity characteristics of the ring joint connection have been represented by a set composed of a rotational spring (K θR ), an axial spring (K AR ) and a radial spring (K RR ). The interaction mechanism of each spring is the same as the one applied for a segment joint.
Owing to the arbitrary distribution of the joints in the lining perimeter, a full 3D model is necessary.
A full model with a height of 60 m and a width of 120 m has been adopted in the present model. The mesh length of the model is equal to 120 m. This number has not been chosen arbitrarily and will be explained later on. The excavation step length is equal to 1.5 m and corresponds to the width of a lining ring. A schematic view of the present model is provided in Fig. 2. The model introduced in Fig. 3 is composed of around 505,000 grid points and 480,000 zones.
The positions of the segment joints in each ring are presented in Fig. 4 and Table 4. Finally, it should be mentioned that the average time requested for one calculation is approximately 190 hours when using a 2.67 GHz core i7 CPU computer. Figure 5 show a great impact of the constitutive model on the predicted settlements developing on the ground surface and tunnel crown. A difference of about 74% between the maximum surface settlements estimated by means of the CYsoil model and MC model can be seen (Fig. 5a), thus indicating the importance of taking into account the variation in the stiffness modulus. The corresponding difference for the tunnel crown settlement is about 40% (Fig. 5b).

Influence of the Constitutive Model
In order to highlight the influence of constitutive model on the tunnel excavation, the plastic zones developed around the tunnel are illustrated in Fig. 6. As can be seen, the plastic zone surrounding the tunnel in case of using the Cysoil model is approximately twice larger than the one observed in case of using the MC model.    Both the constant stiffness of the MC model and yield characteristic at the small strain of the CYsoil model are the main reasons for the difference in the ground settlement. Figure 5c indicates the impact of the constitutive model on the surface settlement in a transverse vertical plane. Large settlement differences can be observed between the two models at the shield tail and far behind the tunnel face. However, the surface movements in the two models are quiet similar at the tunnel face. Unlike the surface settlement, the movement at the tunnel crown indicates a considerable difference between the two models at each stage (Fig. 5d).
In the same way as for the settlements, Fig. 7 indicates n important difference between the structural forces, especially the bending moment, longitudinal forces and normal displacement induced in a tunnel lining, using the MC model and those obtained in the case of the CYsoil model. All these values are determined at the final state of the tunnel lining.
As expected, due to the higher development of the plastic zones surrounding the tunnel in case of using the CYsoil model, the normal forces induced in the tunnel lining are larger than the corresponding ones obtained in the case of using the MC model (Fig. 7b).
This could be attributed to the higher external loads which act on the tunnel lining in the case of the CYsoil model, caused by the self-weight of the failure ground in the surrounding plastic zone. Also, due to the higher magnitude of the external loads originating from the surrounding ground, a higher valueof the longitudinal force induced the lining is observed. This is caused by the partial restraint of the transversal deformation (Poisson effect).Lower values of the longitudinal forces in the case of the CYsoil model that appear locally at the tunnel crown could be explained by the larger vertical displacement of the tunnel lining which have occurred at the tunnel crown (Fig. 7d).
Similarly, the CYsoil model produces higher bending moments and normal displacements than those predicted by the MC model (Fig. 7a, d). The largerdifferences are observed at the crown and bottom of the tunnel and at the spring lines.
The above results reveal that the MC model does not adequately reflect realistic ground behavior during tunneling, due to the fact that it results induce lower structural forces. Then the design of the tunnel structure would therefore be unsafe. More results are illustrated in Table 5. The maximum differences in the results can be seen for the normal displacement and for the bending moment. The minimum difference is obtained for the normal forces. Considering the disadvantages of the MC model introduced in previous studies (Bolton et al., 1994;Addenbrooke et al., 1997;Masin and Herle, 2005;Hejazi et al., 2008), the above results show the importance of the complexity of aconstitutive model in tunnelling engineering.

CONCLUSION
A three-dimensional numerical model of the mechanized tunnelling process has beendeveloped. This model can be used to predict the ground movements and structural forces induced in a tunnel lining.
The results have shown agreatinfluence of the constitutive model of the ground on the tunnel behaviour and ground displacement.Generally, the CYsoil model provides higher structural forces and ground settlement than those predicted by the MC model. These results indicate the importance of using a sufficiently complex constitutive model to take into account the real soil behaviour in tunnelling engineering.