FEM modeling requirements for accurate analysis of highly nonlinear shallow tunnels

Abstract

Modern tunnel design in urban areas heavily relies on numerical modeling to assess excavation stability and predict ground movement. Recent advancements in soil modeling, hardware, and software have facilitated the development of sophisticated 3D models within tight schedules. Urban tunnels are often shallow and excavated in challenging ground conditions, with proximity to existing structures and infrastructure. Consequently, numerical modeling of such tunnels involves highly nonlinear analyses with complex boundary conditions. Despite the widespread use of numerical modeling in tunnel research and design, there is a lack of publications addressing modeling procedures to ensure accurate and reliable results for highly nonlinear shallow tunnel analyses. This paper investigates the requirements for accurate results for highly nonlinear shallow tunnel analyses. The Finite Element Method (FEM) is employed with different mesh refinements and element types. The study focuses on the hypothetical excavation stability scenario explored by Carranza-Torres et al. (2013). Tunnel stability is assessed using Caquot's analytical solution based on the lower bound theorem of plasticity, as well as FEM modeling with the strength reduction method. The FEM numerical solution, which approaches the exact solution for the problem, provided a factor of safety slightly larger (2.3%) than Caquot's lower-bound solution. The results of the FEM modeling indicate that a significantly less refined mesh is required to achieve accurate results for highly nonlinear shallow tunnel analyses when adopting 2^nd-order elements (i.e., quadratic interpolation) instead of 1^st-order elements (i.e., linear interpolation). This study improves our understanding of FEM modeling requirements and provides practical insights regarding the numerical modeling of highly nonlinear shallow tunnels in urban areas.

Keywords:
Shallow tunnel; Finite Element Method; Numerical modeling; Face stability; Tunnel stability

1. Introduction

The increasing demand for underground space in urban areas has posed significant challenges in designing tunnels and deep excavations due to unfavorable geological conditions and their proximity to existing infrastructure. To tackle these challenges, modern tunnel design in urban areas relies heavily on numerical modeling, leveraging recent advancements in soil modeling, software, and hardware. However, modeling procedures and verification guidelines for complex problems are scarce in this context. Consequently, numerical predictions can vary significantly among consulting firms, as Schweiger (2002) discussed, highlighting the need for comprehensive guidelines that are often lacking or not followed.

Ensuring numerical accuracy is crucial, particularly for highly nonlinear problems with complex boundary conditions. Shallow tunnel modeling, characterized by low cover-to-diameter ratios and weak soils, is typically highly nonlinear, leading to large plastic strains around the tunnel. Stability analysis for shallow tunnels is commonly conducted using the strength reduction method (SRM), where soil strength is gradually reduced until failure. In SRM, failure is indicated by the divergence of the numerical solution, signifying a violation of equilibrium. The strength reduction factor (SRF) at failure corresponds to the traditional factor of safety in geotechnical design. However, near collapse, the numerical analysis becomes highly nonlinear, posing challenges for numerical accuracy.

This paper addresses the modeling requirements to ensure accuracy in highly nonlinear shallow tunnel analyses using the implicit finite element method (FEM). The investigation focuses on a hypothetical scenario previously studied numerically and analytically by Carranza-Torres et al. (2013). The tunnel stability is assessed using Caquot's analytical solution, which is based on the lower bound theorem of plasticity and encompasses the limit state of equilibrium. Carranza-Torres et al. (2013) extended Caquot's solution and conducted Finite Difference Method (FDM) analyses for such hypothetical scenario. In this study, the FEM is utilized to explore the same scenario and to provide insights into the requirements for accurate numerical solutions.

Mesh refinement and the utilization of high-order elements have been recognized as methods to enhance accuracy in numerical simulations, albeit at the cost of increased computational resources. In previous studies by Vitali et al. (2018a, 2021a), mesh optimization investigations were conducted for deep tunnel modeling, leading to the recommendation of employing higher-order elements due to their efficiency in handling highly nonlinear problems. Taking into account the dimensions and mesh refinement guidelines proposed by Vitali et al. (2018a), a similar approach was adopted for modeling both shallow and deep tunnels under complex anisotropic conditions (e.g. Vitali et al., 2018b, 2019a, b, c, d, 2020a, b, c, 2021b, c, 2022). By adhering to these recommendations, the authors achieved a remarkable agreement between numerical and analytical solutions, as well as with the data obtained from instrumentation in real tunnel scenarios. However, their study did not address the mesh requirements for simulating the failure of shallow tunnels, where the numerical problem becomes significantly nonlinear, with extremely large plastic strains nearing tunnel collapse.

Therefore, this paper fills the existing research gap by examining the requirements for accurate numerical solutions in highly nonlinear shallow tunnel analyses. The work provides valuable insights derived from the implicit FEM analysis of the hypothetical scenario of Carranza-Torres et al. (2013). The accuracy of the numerical solutions is evaluated against Caquot's analytical solution and the FDM model conducted by Carranza-Torres et al. (2013), allowing for a comparison and discussion of the discrepancies between analytical and numerical approaches.

2. FEM modeling and results

In this investigation, the hypothetical scenario of an infinite tunnel previously examined by Carranza-Torres et al. (2013) was selected. The scenario, including its dimensions, loads, and soil properties, is illustrated in Figure 1. The problem is treated as a 2D plane-strain analysis. The tunnel has a circular shape with a diameter of 2 m, and its crown is situated 2 m below the ground surface. A surcharge of 72 kPa is applied at the ground surface, while a uniform pressure of 36 kPa is exerted along the tunnel perimeter. The soil properties chosen for this scenario represent soft ground conditions, with a friction angle of 30 degrees, cohesion of 9 kPa, unit weight of 18 kN/m³, Young's modulus of 3 MPa, and Poisson's ratio of 0.25.

Figure 1
Hypothetical scenario from Carranza-Torres et al. (2013).

Based on the analytical solution developed by Caquot (1934) and extended by Carranza-Torres et al. (2013), the factor of safety for this specific scenario is determined to be 1.71. In their study, Carranza-Torres et al. (2013) utilized the strength reduction method within the finite difference method (FDM) software FLAC3D and obtained a factor of safety of 2.01. They attributed the discrepancy between the analytical and numerical results to the conservative nature of Caquot's analytical solution, which is based on the lower bound theorem of plasticity, while the numerical solution should approach the exact solution.

In the current study, the hypothetical scenario was modeled using the implicit finite element method (FEM) with the software Midas GTS NX. The analysis was conducted in a 2D plane strain configuration, taking advantage of the problem's symmetry by modeling only half of the problem. Four different structured mesh refinements were investigated, namely: 8 divisions at half tunnel perimeter (Figure 2a), 16 divisions (Figure 2b), 32 divisions (Figure 2c), and 64 divisions (Figure 2d). The mesh refinement near the tunnel, as shown in Figure 2b, corresponds to the recommendation made by Vitali et al. (2018a). The models were analyzed using 1^st and 2^nd-order elements, representing linear and quadratic interpolation. To facilitate comparisons, the size of the models was kept the same as that adopted by Carranza-Torres et al. (2013). Mesh discretization far from the tunnel was not investigated as it is likely to have minimal influence on the numerical accuracy of the model.

Figure 2
Model size and mesh refinement. (a) 8 divisions; (b) 16 divisions; (c) 32 divisions; (d) 64 divisions.

The numerical calculations were carried out in four steps. In the first step, the geostatic stress was generated, assuming a coefficient of earth pressure at rest (K₀) of 1, consistent with Carranza-Torres et al. (2013). In the second step, a surface surcharge of 72 kPa was applied, and the displacements in the medium were reset to zero. The tunnel excavation was simulated in the third step by deactivating the elements inside the tunnel while maintaining the nodal forces balancing the stresses at the tunnel perimeter. In the same step, the forces at the tunnel perimeter were gradually reduced to achieve a uniform internal pressure of 36 kPa. In the fourth step, the strength properties of the ground were progressively reduced until failure. This approach, known as the Strength Reduction Method (SRM), was implemented using the Midas GTS NX software. In the SRM method, failure is associated with the inability of the software to balance stresses, leading to divergence in the numerical solution due to the violation of equilibrium. The ground was modeled using an elastic-perfectly plastic model with the Mohr-Coulomb failure criterion with an associated flow rule, as both lower and upper-bound solutions were based on this flow rule.

Figure 3 illustrates the relationship between the strength reduction factor (SRF) and the normalized vertical displacement at the ground surface. Figure 3a displays the results obtained with the 2^nd-order elements using the meshes shown in Figure 2, while Figure 3b presents the results with 1^st-order elements. The factors of safety obtained from the extended Caquot's analytical solution (FS = 1.71) and the FDM model by Carranza-Torres et al. (2013), FS = 2.01, are also included in Figure 3.

Figure 3
Strength reduction factor (SRF) plotted against the normalized vertical displacement at the ground surface above the crown (Sy) relative to the displacement for SRF=1 (Sy, SRF=1). Results are shown for meshes using 2^nd-order elements (a) and 1^st-order elements (b).

In Figure 3a, it is evident that the displacements at the surface increase non-linearly as the shear strength properties of the ground are progressively reduced. The SRF at failure obtained with the 16 and 32 division meshes using 2^nd-order elements are slightly higher than Caquot's lower bound solution (i.e., 1.77 and 1.75, respectively). Further refining the mesh beyond 16 divisions (Figure 2b) does not yield significant improvements, as indicated by the negligible difference observed between the 16-division and 32-division meshes using 2^nd-order elements (1% difference). Consequently, the 64-division mesh (Figure 2d) was not analyzed with 2^nd-order elements. The 16-division mesh with 2^nd-order elements aligns with the mesh refinement recommendation from Vitali et al. (2018a). The SRF at failure for the 8-division mesh with 2^nd-order elements (Figure 2a) is 1.83, which is slightly larger (4.6% larger) than the SRF at failure for the 32-division mesh.

In contrast, when using 1^st-order elements, the SRF at failure varies considerably with the mesh refinement, as shown in Figure 3b. As the 1^st-order element mesh becomes finer, the SRF at failure approaches Caquot's lower bound solution. For the coarser mesh (8 divisions, Figure 2a), the SRF increases non-linearly up to 2.5. Beyond this point, larger displacement increments occur as the strength is reduced, leading to the inability of the software to achieve equilibrium for an SRF of 2.93 (67% larger than the SRF at failure obtained with the 32-division mesh using 2^nd-order elements). Using 2^nd-order elements, the same mesh provides more reliable results, with an SRF at failure of 1.83. The SRF at failure for the 16-division mesh with 1^st-order elements (Figure 2b) is 2.18, which is 24.6% larger than the SRF at failure obtained with the 32-division mesh using 2^nd-order elements.

Interestingly, the most refined mesh using 1^st-order elements (i.e., 64 divisions, Figure 2d) yields an SRF at failure of 1.80, which is close to the value obtained with the 16-division mesh using 2^nd-order elements (i.e., 1.77, Figure 2b). These results indicate that significantly less refined meshes are sufficient to achieve accurate numerical results for highly nonlinear shallow tunnel analyses when 2^nd-order elements are employed. Additionally, it is worth noting that even very refined meshes using 1^st-order elements may not provide accurate results for highly nonlinear shallow tunnel analyses (e.g., the 32-division mesh with 1^st-order elements, Figure 2c, yielded a factor of safety of 1.90, 8.6% larger than the value of 1.75 obtained with the 32-division mesh using 2^nd-order elements).

Carranza-Torres et al. (2013) employed a highly refined grid with 30 divisions around the half-tunnel perimeter. The FDM solution from Carranza-Torres et al. (2013) yielded a larger factor of safety compared to Caquot's lower bound solution and the FEM solution with 32 divisions using 2^nd-order elements (i.e., 2.01 instead of 1.71 and 1.75, respectively). Factors of safety around 2 were obtained with 16 and 32 divisions using 1^st order elements (Figures 2b and 2c), indicating that the FDM analysis may be comparable to the FEM analysis using 1^st-order elements.

Figure 4 depicts the distribution of equivalent plastic strains for various strength reduction factors. The plots were generated from the mesh with 32 divisions (Figure 2c) using 2^nd-order elements. The plastic strains are limited to 10% to enhance visualization of the plastic zone. Equation 1 provides the expression for the equivalent plastic strain.

Figure 4
Evolution of the plastic strain field with the strength reduction factor (SRF) up to 1.7525, representing the maximum SRF value before divergence.

${\epsilon }_p={\displaystyle \int \sqrt{\frac{2}{3}{\dot{\epsilon }}_{ij}{\dot{\epsilon }}_{ij}}}dt$ (1)

For an SRF of 1, most of the ground surrounding the tunnel remains within its elastic regime. However, a localized plastic zone appears at the tunnel springline. As the strength properties are reduced by a factor of 1.2, the plastic strains intensify, and the size of the plastic zone expands. Figures 4c and 4d illustrate that as the strength reduction factor increases to 1.4 and 1.6, the plastic zone propagates from the tunnel springline toward the ground surface. Once the plastic zone reaches the ground surface, it continues to grow upwards, extending above the tunnel crown, as the ground strength properties further diminish.

Beyond an SRF of 1.75, the non-linear solution necessitates a significantly larger number of iterations to converge. Additionally, even minor reductions in strength lead to substantial displacement increments, as observed in Figure 3a. The non-linear solution fails to converge when the SRF reaches 1.7525. Consequently, the factor of safety of 1.75 likely approximates the exact solution of the problem, which surpasses Caquot's lower bound solution by 2.3%. These findings suggest that the lower bound solution for plane strain shallow tunnel stability analysis closely approximates the exact solution.

3. Conclusion

In conclusion, this study employed the finite element method (FEM) to investigate the hypothetical scenario of tunnel instability initially explored by Carranza-Torres et al. (2013). A comparison was made between FEM models utilizing different mesh refinements and element orders and the analytical solution for shallow circular tunnels proposed by Caquot (1934) and extended by Carranza-Torres et al. (2013). Additionally, the results were compared with those obtained from the finite difference method (FDM) model used by Carranza-Torres et al. (2013).

The analytical solution, based on the lower bound theorem of plasticity, provides the factor of safety for shallow tunnels considering the ground's strength properties, a uniform surcharge at the surface, and a constant internal pressure at the tunnel perimeter. To determine the factor of safety, the strength reduction method (SRM) implemented in the FEM software was employed. Four different mesh refinements were studied, consisting of 8, 16, 32, and 64 divisions along half of the tunnel perimeter.

The findings indicate that the numerical solution using FEM and SRM, with appropriately refined meshes and 2^nd-order elements, yields a factor of safety slightly larger than Caquot's lower bound solution (2.3% larger, to be specific). With 2^nd-order elements, the select meshes demonstrate similar strength reduction factors (SRF) at failure, with the 32-divisions mesh resulting in an SRF at failure of 1.75 and the 8-divisions mesh yielding an SRF at failure of 1.83, a difference of 4.6%.

However, when 1^st-order elements are used, the SRF at failure shows substantial variation depending on the mesh refinement. This contrasts with the consistent behavior observed with 2^nd-order elements. For instance, the 64-divisions mesh with 1st-order elements produces an SRF at failure of 1.80, while the 8-divisions mesh results in an SRF of 2.93. Notably, the most refined mesh with 1^st-order elements (i.e., 64 divisions along the half tunnel perimeter) generates results similar to those obtained with the 16-division mesh utilizing 2^nd-order elements.

These findings emphasize the criticality of utilizing 2^nd-order elements to achieve precise and reliable analyses for highly nonlinear shallow tunnels. Furthermore, the FEM results presented in this study suggest that the exact solution for shallow tunnel instability analysis closely aligns with Caquot’s lower bound solution.

List of symbols

cCohesion

EYoung’s modulus

FDMFinite Difference Method

FEMFinite Element Method

FSSafety Factor

K₀ Coefficient of earth pressure at rest

p_s Uniform internal pressure at tunnel perimeter

q_s Surface surcharge

SRFStrength Reduction Factor

SRMStrength Reduction Method

ε_ijStrain tensor

ε_pEquivalent plastic strain

φFriction angle

γUnit weight

νPoisson’s ratio

Data availability

The datasets generated and analyzed in the course of the current study are available from the corresponding author upon request.

Acknowledgements

The authors acknowledge the support provided by the research funding agency of the Brazilian Government, CNPq (Conselho Nacional de Desenvolvimento Cientifico), and CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior). The authors also thank the Midas company for providing the license for the Midas GTS NX software utilized in this research.

• Discussion open until May 31, 2024.

References

• Caquot, A. (1934). Équilibre des massifs a frottement interne. Gauthier-Villars.
• Carranza-Torres, C., Reich, T., & Saftner, D. (2013). Stability of shallow circular tunnels in soils using analytical and numerical models. In 61st Minnesota Annual Geotechnical Engineering Conference Falcon Heights: University of Minnesota.
• Schweiger, H.F. (2002). Results from numerical benchmark exercises in geotechnics. In 5th European Conference on Numerical Methods in Geotechnical Engineering (pp. 305-314). Paris: Presses Ponts et chaussees.
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2018a). 3D finite element modeling optimization for deep tunnels with material nonlinearity. Underground Space, 3(2), 125-139. http://dx.doi.org/10.1016/j.undsp.2017.11.002
  » http://dx.doi.org/10.1016/j.undsp.2017.11.002
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2018b). Analytical solution for tunnels not aligned with geostatic principal stress directions. Tunnelling and Underground Space Technology, 82, 394-405. http://dx.doi.org/10.1016/j.tust.2018.08.046
  » http://dx.doi.org/10.1016/j.tust.2018.08.046
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2019a). Shallow tunnel not aligned with the geostatic principal stress directions. In Geo-Congress 2019 (pp. 214-222). Philadelphia: ASCE. http://dx.doi.org/10.1061/9780784482155.023
  » http://dx.doi.org/10.1061/9780784482155.023
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2019b). Shallow tunnels misaligned with geostatic principal stress directions: analytical solution and 3D face effects. Tunnelling and Underground Space Technology, 89, 268-283. http://dx.doi.org/10.1016/j.tust.2019.04.006
  » http://dx.doi.org/10.1016/j.tust.2019.04.006
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2019c). Progressive failure due to tunnel misalignment with geostatic principal stresses. In ISRM 14th International Congress on Rock Mechanics (pp. 2292-2299). Foz do Iguaçu: ISRM.
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2019d). Buoyancy effect on shallow tunnels. International Journal of Rock Mechanics and Mining Sciences, 114, 1-6. http://dx.doi.org/10.1016/j.ijrmms.2018.12.012
  » http://dx.doi.org/10.1016/j.ijrmms.2018.12.012
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2020a). Analytical solution for a deep circular tunnel in anisotropic ground and anisotropic geostatic stresses. Rock Mechanics and Rock Engineering, 53(9), 3859-3884. http://dx.doi.org/10.1007/s00603-020-02157-5
  » http://dx.doi.org/10.1007/s00603-020-02157-5
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2020b). Tunnel misalignment with geostatic principal stress directions in anisotropic rock masses. Soils and Rocks, 43(1), 123-138. http://dx.doi.org/10.28927/SR.431123
  » http://dx.doi.org/10.28927/SR.431123
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2020c). Deformation patterns and 3D face effects of tunnels misaligned with the geostatic principal stresses in isotropic and anisotropic rock masses. In 54th US Rock Mechanics /Geomechanics Symposium (ARMA 2020) (pp. 1-9). Golden: ARMA.
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2021a). New modeling approach for tunnels under complex ground and loading conditions. Soils and Rocks, 44(1), 1-8. http://dx.doi.org/10.28927/SR.2021.052120
  » http://dx.doi.org/10.28927/SR.2021.052120
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2021b). Behavior of tunnels excavated with dip and against dip. Underground Space, 6(6), 709-717. http://dx.doi.org/10.1016/j.undsp.2021.04.001
  » http://dx.doi.org/10.1016/j.undsp.2021.04.001
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2021c). Construction strategies for a NATM tunnel in São Paulo, Brazil, in residual soil. Underground Space, 7(1), 1-18. http://dx.doi.org/10.1016/j.undsp.2021.04.002
  » http://dx.doi.org/10.1016/j.undsp.2021.04.002
• Vitali, O.P.M., Celestino, T.B., & Bobet, A. (2022). 3D face effects of tunnels misaligned with the principal directions of material and stress anisotropy. Tunnelling and Underground Space Technology, 122, 104347. http://dx.doi.org/10.1016/j.tust.2021.104347
  » http://dx.doi.org/10.1016/j.tust.2021.104347

Publication Dates

• Publication in this collection
  02 Feb 2024
• Date of issue
  2024

History

• Received
  21 Jan 2023
• Accepted
  18 Dec 2023