Use of finite element 3D analyses of tunnelling induced settlements

tunnelling Results of a series of 3D and 2D analyses of open-face tunnelling in marly-clayey deposits are presented in the paper. A parametric study was performed to investigate the influence of initial stress state, soil deformability, shear strength parameters, and soil anisotropy, on settlement prediction. Comparison of 2D and 3D results shows that the settlement trough predicted by 2D analysis agrees well with 3D results when an appropriate amount of unloading prior to lining installation is adopted. It is demonstrated that the load reduction factor significantly depends on shear strength parameters of soil.


The main issue of open-face tunnelling (NATM -New Austrian
Tunnelling method or open-face shield) in urban environment is the control of surface settlements. Predicting the settlements induced by tunnelling in soft ground is an important task of geotechnical design, and as such it has been analysed in the scope of many studies. There are three categories for approaching this problem: empirical methods, analytical methods, and numerical methods. Empirical methods, which are based on assumption that the settlement trough can be described by a Gaussian function [1], are often used in engineering practice. They provide reasonable prediction of surface settlements if site conditions are well known, and if design parameters are appropriately calibrated. Analytical methods [2][3][4][5] provide simple solutions in closed form, but their practical application is limited, since they are often based on idealized assumptions. Advanced numerical methods are required for modelling complex tunnelsoil interaction problems. The finite element method (FEM) is a flexible tool that has been adopted by many authors. Various studies have shown that a crucial role in FE prediction of tunnelling-induced settlement is assumed by the 3D effects of tunnelling, the initial stress state (particularly the coefficient of lateral earth pressure at rest K 0 ) and the stress-strain behaviour of soil [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. It has been noted by several authors that the finite element analysis of tunnelling predicts excessively wide transverse settlement troughs, when compared with field data, especially in stiff clays under high K 0 conditions (K 0 > 1). It has been suggested that improved prediction can be achieved by modelling soil anisotropy and small-strain nonlinearity. Simpson et al. [8] concluded from FE analyses of tunnelling in London Clay that the width of the calculated surface settlement trough is substantially influenced by the shear modulus anisotropy, while it is much less influenced by non-linearity (both the orthotropic linear elastic perfectly plastic Mohr-Coulomb model and the orthotropic non-linear model gave settlement troughs that were much closer to the observed shape than the non-linear but isotropic model). Lee and Rowe [6] concluded from linear elastic perfectly plastic FE analyses of tunnelling that the effect of elastic anisotropy, in particular the ratio of the independent shear modulus G vh to vertical modulus E v , should be considered if reasonable predictions of settlement are to be obtained. Tunnel construction is clearly a three-dimensional process. A full 3D FE analysis is required to simulate progress of tunnelling work and stress-strain changes that take place at the temporary working face. However, because of computational effort involved in the case of 3D modelling, 2D models are still commonly used in routine geotechnical design. When the process of tunnel construction is modelled in plane strain, at least one assumption must be made in order to account for stress relief and ground movements occurring at the tunnel face prior to lining installation. Various methods that take into account 3D effects of tunnelling within simplified 2D plane strain analysis have so far been proposed in literature (review of methods can be found in [21][22]). The most commonly used method for 2D modelling of open-face tunnel construction is the stress (load) reduction method, which is actually a FE utilisation of the convergence-confinement method [23]. A series of 2D and 3D FE elastoplastic analyses of openface tunnelling works in Belgrade marly-clayey deposits was performed in this study. The transverse surface settlement troughs predicted by 3D and 2D FE modelling were compared and values of the stress reduction factor were calibrated. A parametric study was performed to investigate the influence of the coefficient of earth pressure at rest K 0 , soil stiffness, shear strength parameters (cohesion and friction angle), and soil anisotropy on settlement prediction. Within the conducted parametric studies, some conclusions are given about the stress reduction factor. The FE analysis of settlement above the Dedinje Tunnel, which is part of the Belgrade railway junction, is presented as a case study.

Finite element analysis
The finite element mesh adopted for full 3D modelling of open-face tunnelling works in Belgrade marly-clayey deposits is shown in Figure 1. These deposits are calcareous silty clays of stiff consistency and high plasticity. The adopted soil profile was assumed to be representative of soil in Belgrade area. It is composed of a 5 m thick layer of loess and a 15m thick layer of weathered brownish-yellow marly clay overlaying grey unaltered marly clay and marl. The tunnel was assumed to be 6m in diameter with the tunnel lining thickness of 0.35 m (adopted from the Feasibility Study for Rapid Public Urban Transport in Belgrade [24]). The tunnel axis was assumed to be at a depth of 15 m.

Figure 1. FE mesh adopted for 3D analyses
The finite element analyses presented in this paper were carried out using DIANA Finite Element Analysis code (TNO DIANA BV, Delft). DIANA is a multi-purpose finite element program (threedimensional and nonlinear) that enables phased construction modelling [25].

Modelling of tunnel construction
The 3D process of tunnel construction was simulated using a stepby-step approach [26][27][28]. The analysis started from the in-situ GRAĐEVINAR 72 (2020) 8, 673-680 Use of finite element method for 2D and 3D analyses of tunnelling induced settlements stress state. The tunnel construction was modelled by successive removal of elements in front of the tunnel face to simulate an unsupported excavation with round length of 2m, while successively installing lining elements to support the previous excavation. The process of tunnel construction was simulated in 40 steps, over a length of 80 m. The soil was modelled by 20-node isoparametric solid brick elements, whereas 8-node quadrilateral isoparametric curved shell elements were used to represent the tunnel lining. No interface elements were used between the lining and the soil (the assumption was that there is no relative movement between the lining and adjacent soil). Normal horizontal movements were restricted at all vertical boundaries of the mesh, whereas movements in all directions were restricted for the base of the mesh. An additional condition involving prevention of rotation around longitudinal axis in the plane of symmetry was applied for shell elements of the tunnel lining. All 2D FE analyses were performed in this study using the stress (load) reduction method, which is actually a FE utilisation of the convergence-confinement method [23]. In this method, the lining is installed after a certain percentage of the initial stresses is released. The stress reduction factor l (the proportion of unloading before lining installation) is the controlling parameter that has to be prescribed. Figure 2 shows the FE mesh adopted for 2D analyses, which corresponds to a cross section of the 3D mesh shown in Figure 1. The soil was modelled by 8-node quadrilateral isoparametric plane strain elements, whereas tunnel lining was modelled by 3-node curved infinite shell elements.

Figure 2. FE mesh adopted for 2D analyses
Starting from the initial geostatic stress state, the soil elements within the tunnel boundary were removed and the lining was installed at a prescribed value l on the deformed tunnel boundary. The parameter l depends on a number of factors such as soil behaviour, tunnel geometry, construction method and round length. It can be calibrated based on the comparison of 2D and 3D results. In practice, parameter l is often estimated based on engineering experience with similar tunnelling conditions or monitoring data. In this study, values of the stress reduction factor l were calibrated to obtain the same magnitude of surface settlement above the tunnel axis by 2D and 3D analysis.

Initial conditions, constitutive model and model parameters
Initial soil conditions were established at the start of each FE analysis. Initial stresses were specified using the values of the mass density and the coefficient of lateral earth pressure at rest K 0 listed in Table 1. Given the hydro-geological characteristics of the terrain in Belgrade, the presence of a groundwater table was not considered. In order to investigate the influence of K 0 on surface settlement, FE analyses were performed with K 0 = 0.65, 0.85 and 1.0 for the second ground layer (given that the value of K 0 at tunnel depth has the greatest influence on the results). The drained soil behaviour was modelled in this study using a linear elastic perfectly plastic soil model with Mohr-Coulomb (MC) failure criterion [25]. The intention was to investigate the influence of soil modulus of elasticity E and shear strength parameters (cohesion c' and friction angle f') on predicted surface settlements. The parameters of the MC model, assumed to be representative of Belgrade soil, are listed in Table 1 (review of geotechnical properties of these soils can be found in [29]). It is adopted that the angle of dilation y' is equal to zero for all soil layers. In the parametric study, the values of E, f' and c' were varied for the second soil layer (weathered brownish-yellow marly clay), in which the tunnel was located. The tunnel lining was modelled by elastic shell elements with the thickness of 0.35 m, r = 2.5 g/cm 3 , E = 15 GPa, and n = 0.15. Additional FE analyses were performed in this study to investigate the effects of elastic anisotropy on predicted surface settlements. The cross-anisotropic (orthotropic) material behaviour is fully described by five independent elastic parameters: vertical Young's modulus E v , horizontal Young's modulus E h , two Poisson's ratios n vh and n hh , and shear modulus are often specified when using these parameters. Lee and Rowe [6] have shown that the ratio of independent shear modulus G vh to vertical modulus E v has a major influence on the shape of the tunnel-induced settlement trough and that the

Results of analysis and discussion
A 3D analysis was initially performed with input parameters given (without parentheses) in Table 1. Figure 3 shows the longitudinal ground surface settlement profile above the tunnel axis, computed after 40 excavation phases (tunnel face at y = -80 m). The steady-state settlement conditions (horizontal part of the longitudinal profile) were developed at approximately 30 m (5D) behind the tunnel face. Figure 4 presents the transverse ground surface settlement profile obtained from 3D analysis at y = -50 m (steady-state conditions) and the profile obtained from equivalent 2D analysis using the stress reduction method with l = 0.63. The stress reduction factor l was calibrated based on the steady-state settlement from 3D analysis (as will be discussed in Section 3.5). It can be observed that 3D and 2D analyses give similar surface settlement troughs when an appropriate value of the stress reduction factor is adopted. The transverse settlement troughs obtained from FE analyses were compared with empirical Gaussian curves. Extensive field measurements have shown that a transverse settlement trough can be well described by a Gaussian function [1], with vertical settlement in transverse direction given by Eq. (1): where S vmax is the maximum settlement measured above the tunnel axis and i is the horizontal distance from the tunnel axis to the point of inflection of the settlement trough.   Use of finite element method for 2D and 3D analyses of tunnelling induced settlements derived). It can be seen that the width of the settlement trough was somewhat overestimated by the finite element method.

Influence of the coefficient of lateral earth pressure
Additional 3D analyses were performed in this study to investigate the influence of the coefficient of lateral earth pressure at rest (K 0 ) on the tunnel-induced settlement trough. Figure 5 shows the transverse surface settlement profiles obtained from 3D analyses (at y = -50 m) with K 0 of 0.65, 0.85 and 1.0 (for the second ground layer in which the tunnel was located). It can be observed that lower values of K 0 give narrower and deeper transverse settlement troughs, which is in line with the findings of other authors.

Influence of soil stiffness
Surface settlement troughs obtained from 3D analyses (at y = -50m) using Young modulus values of E = 10, 15 and 30 MPa (for the second ground layer) are shown in Figure 6. It can be seen that the elastic modulus has a significant influence on the calculated surface settlement. As expected, the lower values of E give higher settlements.   The given graphs indicate that strength parameters, cohesion and angle of shearing resistance, have a significant influence on calculated subsidence. Given that the value of cohesion obtained from laboratory tests commonly vary within wide limits, a particular attention should be given to the determination of the appropriate value of this parameter.

Influence of soil anisotropy
The transverse settlement profile for the isotropic model is compared in Figure 9 with profiles for the anisotropic model with parameters n = E h /E v = 1.2 and m = G vh /E v = 0.2, 0.3 and 0.4 (for the second ground layer). The presented results show that the ratio of the independent shear modulus G vh to vertical modulus E v has a significant effect on the shape of the settlement trough. The surface settlement trough becomes narrower and deeper with an increase in the level of anisotropy.

Findings about the stress (load) reduction factor
2D models are still commonly used in routine geotechnical design. Additional 2D analyses were performed in order to investigate the influence of the initial stress state K 0 , soil stiffness, shear strength parameters, and soil anisotropy, on the magnitude of the stress reduction factor. The stress reduction factor l was calibrated via 3D analysis based on the steady-state settlement.
In the calibration process, a series of 2D analyses was performed with different (progressively increasing) values of l and the obtained surface settlement above tunnel axis was monitored.
The value of parameter lwas adopted to match the magnitude of surface settlement above the tunnel axis obtained by 2D analysis with settlement magnitudes obtained by 3D analysis. Derived values of l are presented in Table 2.

Table 2. Results of parametric study
As is evident from the presented table, the stress reduction factor l significantly depends on shear strength parameters of soil and, at that, the stress reduction factor reduces with an increase in the cohesion and angle of shear strength. This once again points to the importance of determining appropriate values of strength parameters, with particular attention being paid to the determination of cohesion. It was established that the factor l is not substantially influenced by the initial stress state K 0 , soil stiffness, and soil anisotropy.

Case study -Dedinje Tunnel
The Dedinje Tunnel is a twin-tube single-track 2.9 km long tunnel located between railway stations Prokop (Belgrade-Centre) and Rakovica, within the Belgrade railway junction. During the tunnel construction, the vertical ground surface movements were measured at a 220 m long section between chainages 0+850 km and 1+070 km [33]. The tunnel was constructed using the Belgian method (the underpinning arch method) with steel arch ribs and timbers as temporary support. This paper presents the 2D analysis of the greenfield surface settlement above the left tunnel tube, which was constructed first. The subsequent construction of the right tunnel tube is not included in the analysis. Figure 10 shows the 2D finite element mesh used for modelling the left tunnel tube at chainage 1+019 km. The tunnel tube is 5.4 m wide and 7.1 m high (without lining). Its crown is located 14 m below the ground level. The soil profile is composed of marly-clayey sediments covered by loess: completely degraded marls and degraded marls in weathering zone and fresh, unaltered marl under the weathering zone. Since the problem is unsymmetrical, the full geometry was modelled. The soil was modelled by the six-node triangular isoparametric plane strain elements while the lining was modelled by the three-node infinite shell elements. On vertical sides of the mesh, horizontal movements were restricted whereas horizontal and vertical movements were restricted for the base of the mesh. The phased excavation of the tunnel cross section was modelled. Starting from the initial stress state, the top heading was excavated first, and the lining was installed on the deformed tunnel boundary at a prescribed value l. After that, the bench was excavated and the lining was installed. The drained analysis was carried out using input parameters of the MC model given in Table 3 (parameters assumed to be representative of marly clayey sediments of weathering zone were adopted for layers 2 and 3 [29]). The tunnel lining was modelled by elastic shell elements with a thickness of 0.60 m (0.40 m invert), r = 2.4 g/ cm 3 , E = 20 GPa and n = 0.20. Figure 11 shows the surface settlement profiles obtained from FE analyses with l = 0.6 and 0.65, the Gaussian curve (i = 0.5z 0 ) and field data recorded one month after tunnel excavation (from [33]). The stress reduction factor was calibrated based on monitoring data. It can be observed that the shape of the surface settlement trough was predicted with reasonable accuracy.

Conclusion
Tunnel construction is a three-dimensional process. Consideration of stress relief and ground movements that occur at tunnel face prior to the installation of lining is essential if reasonable predictions of settlement induced by tunnelling are to be obtained. The application of full 3D FE modelling enables prediction of surface settlement without the need for an initial assumption regarding volume loss or proportion of unloading before installation of lining. Full 3D FE elastoplastic analysis of an open-face tunnelling work in Belgrade marly clay was performed in this study. The 3D process of tunnel construction was simulated in steps over a length of 80m and the steady-state settlement conditions were developed at approximately 30 m (5D) behind the tunnel face. However, due to the computational effort of 3D modelling, 2D methods are still commonly used In engineering practice. A controlling parameter that accounts for 3D effects of tunnelling must be assumed when the tunnel construction process is modelled in plane strain. The transverse settlement profiles obtained from 3D and 2D FE analyses were compared in this study, and the stress reduction factor values were calibrated. It can be concluded that the settlement trough predicted by 2D analysis agrees well with the 3D results when an appropriate value of the stress reduction factor is adopted. By comparing settlement troughs obtained via FE analyses with the Gaussian curve it can be seen that the width of the settlement trough is somewhat overestimated by the finite element prediction. A series of 3D FE elastoplastic analyses was performed to investigate the influence of the coefficient of earth pressure at rest K 0 , soil modulus of elasticity, strength parameters, and soil anisotropy on settlements induced by tunnelling. The results show that lower values of K 0 give narrower and deeper transverse settlement troughs, which is in line with the findings of other authors. Also, it is shown that the modulus of elasticity and strength parameters have a substantial effect on the predicted settlement magnitude. Soil anisotropy, i.e. the ratio of the independent shear modulus G vh to vertical modulus E v , was found to have significant influence on the shape of settlement trough. It can be concluded that the surface settlement trough becomes narrower and deeper with an increase in the level of anisotropy. Some conclusions about the stress reduction factor can be given in the scope of parametric studies, based on systematic comparison of the 3D and 2D results. The presented results suggest that the stress reduction factor significantly depends on shear strength parameters of soil, while its value increases with a decrease in the cohesion or angle of shearing resistance. This is because the weaker soil leads to higher deformations at tunnel face prior to the installation of lining. Given that the value of cohesion obtained from laboratory tests commonly varies within wide limits, a particular attention should be paid to the determination of an appropriate value of this parameter. The initial stress state K 0 , soil stiffness, and soil anisotropy, were not found to influence l substantially. Finally, the 2D FE analysis of settlements above the Dedinje Tunnel, which is part of the Belgrade railway junction, was performed and the load reduction factor was calibrated based on monitoring data. It is shown that the shape of the settlement trough obtained by 2D analysis is in a reasonable agreement with the Gaussian curve and monitoring data. This work indicates that FE analysis can produce reasonable prediction of tunnel-induced surface settlements, provided that appropriate values of soil parameters are determined.