New modeling approach for tunnels under complex ground and loading conditions

Rock masses may present pronounced stress and material anisotropy. Stress measurements in rock masses show that stress anisotropy may be highly pronounced, as shown by McGarr & Gay (1978). Also, the data compiled by Brown & Hoek (1978) show that large horizontal stresses are common at shallow depths. According to Brady & Brown (2006), the major horizontal stress (σH) and the minor horizontal stress (σh) rarely have the same magnitude. Large and highly anisotropic horizontal stresses were reported by Haimson et al. (2003) and Park et al. (2014) in gneissic and granitic rock masses in South Korea. Those stresses were associated with the intense tectonic activity in the area. Rock masses may present pronounced fabric structure; thus, anisotropic mechanical behavior is expected. The data compiled by Worotnicki (1993) in metamorphic rocks showed that the ratio between the Young modulus perpendicular to the rock structure and parallel to the rock structure was larger than 2 for more than 50% of the rocks tested (e.g. schists, slates, quartzites, mudstones and phyllites), and the largest stiffness ratio was 6. This is relevant because anisotropic rock properties strongly affect the behavior of tunnels and should be considered in tunnel design (Fortsakis et al., 2012; Wittke, 1990; Armand et al., 2013; Bobet, 2011, 2016; Bobet & Yu, 2016; Vitali et al., 2020a, b, c; Vitali, 2020; Goricki et al., 2005; Schubert & Mendez, 2017, Klopčič & Logar, 2014). In anisotropic rock masses, the tunnel alignment with one of the principal directions of stress and material anisotropy is unlikely. In this case, asymmetric displacements are induced near the face and anti-symmetric axial displacements occur far-behind the face of the tunnel (Vitali et al. 2019b, 2020a, c; Vitali, 2020). The asymmetric displacements near the face affect the performance of the support and rock surrounding the excavation and may produce asymmetric plastic deformations around the tunnel (Vitali et al., 2019b, c, 2020a). Further, asymmetric displacements and asymmetric failure at the tunnel walls are commonly observed (Schubert & Budil, 1995; Goricki et al., 2005; Schubert et al., 2005; Schubert & Moritz, 2011; Klopčič & Logar, 2014; Lenz et al., 2017). 2D analyses cannot capture the 3D face effects that occur in tunnels during construction and, in particular, when the tunnel axis is not one of the principal directions of material anisotropy or a principal far-field stress; thus, 3D analyses are required. Because of recent advances in hardware and software, 3D FEM modeling is nowadays possible in the practice of engineering. However, the numerical modeling of tunnels not aligned with one of the principal directions of material anisotropy may be cumbersome and time consuming. Abstract The behavior of tunnels in anisotropic rock masses is highly complex and heavily dependent on the orientation of the tunnel axis with respect to the geostatic principal stress directions and to the rock structural planes. 2D solutions cannot capture the 3D face effects of such complex scenario; thus, 3D numerical modeling is required. The modeling of such tunnels using conventional boundary conditions may be cumbersome since the tunnel may not be parallel to the boundaries. The issue is further complicated if the principal far-field stresses are not parallel to the principal axes of material anisotropy. In this case, the use of conventional boundary conditions may be problematic. In this paper, a new approach is presented to impose the boundary conditions and the far-field stresses on 3D numerical models of tunnels under complex ground and loading conditions. With the proposed approach, it is possible to easily simulate any orientation of the tunnel with respect to the principal directions of stress and material anisotropy. The numerical results obtained with the proposed approach were validated with an analytical solution and with numerical results using traditional boundary conditions.


Introduction
Rock masses may present pronounced stress and material anisotropy. Stress measurements in rock masses show that stress anisotropy may be highly pronounced, as shown by McGarr & Gay (1978). Also, the data compiled by Brown & Hoek (1978) show that large horizontal stresses are common at shallow depths. According to Brady & Brown (2006), the major horizontal stress (σ H ) and the minor horizontal stress (σ h ) rarely have the same magnitude. Large and highly anisotropic horizontal stresses were reported by Haimson et al. (2003) and Park et al. (2014) in gneissic and granitic rock masses in South Korea. Those stresses were associated with the intense tectonic activity in the area. Rock masses may present pronounced fabric structure; thus, anisotropic mechanical behavior is expected. The data compiled by Worotnicki (1993) in metamorphic rocks showed that the ratio between the Young modulus perpendicular to the rock structure and parallel to the rock structure was larger than 2 for more than 50% of the rocks tested (e.g. schists, slates, quartzites, mudstones and phyllites), and the largest stiffness ratio was 6. This is relevant because anisotropic rock properties strongly affect the behavior of tunnels and should be considered in tunnel design (Fortsakis et al., 2012;Wittke, 1990;Armand et al., 2013;Bobet, 2011Bobet, , 2016Bobet & Yu, 2016;Vitali et al., 2020a, b, c;Vitali, 2020;Goricki et al., 2005;Schubert & Mendez, 2017, Klopčič & Logar, 2014. In anisotropic rock masses, the tunnel alignment with one of the principal directions of stress and material anisotropy is unlikely. In this case, asymmetric displacements are induced near the face and anti-symmetric axial displacements occur far-behind the face of the tunnel (Vitali et al. 2019b(Vitali et al. , 2020aVitali, 2020). The asymmetric displacements near the face affect the performance of the support and rock surrounding the excavation and may produce asymmetric plastic deformations around the tunnel (Vitali et al., 2019b(Vitali et al., , c, 2020a. Further, asymmetric displacements and asymmetric failure at the tunnel walls are commonly observed (Schubert & Budil, 1995;Goricki et al., 2005;Schubert et al., 2005;Schubert & Moritz, 2011;Klopčič & Logar, 2014;Lenz et al., 2017).
2D analyses cannot capture the 3D face effects that occur in tunnels during construction and, in particular, when the tunnel axis is not one of the principal directions of material anisotropy or a principal far-field stress; thus, 3D analyses are required. Because of recent advances in hardware and software, 3D FEM modeling is nowadays possible in the practice of engineering. However, the numerical modeling of tunnels not aligned with one of the principal directions of material anisotropy may be cumbersome and time consuming.
Traditionally, for deep tunnels, the geostatic stress field is generated by applying a uniform pressure perpendicular to the boundaries of the model. With this approach, different 3D FEM meshes need to be created for each attempted orientation of the tunnel with respect to the principal stress directions (Vitali et al., 2018b). Also, if the principal stress directions are not aligned with the principal material directions, displacements parallel to the boundaries are induced, which may be problematic.
In this paper, a new approach for general numerical modeling of tunnels under complex anisotropic conditions is presented. The basic idea is to impose body forces to all the FEM elements to generate the geostatic stress field and to constrain the displacements at the boundaries. Because no displacements are expected far from the tunnel, fixing the nodes at the boundaries is acceptable, with the assumption that the boundaries are sufficiently far from the tunnel. The paper shows that the numerical results obtained with the 3D FEM model imposing the proposed boundary conditions match the analytical results (Vitali et al., 2020b) and the results of 3D FEM models with traditional boundary conditions (Vitali et al., 2020a). Figure 1 shows the 3D FEM mesh used with the proposed boundary conditions. The tunnel is assumed deep. The geostatic stress field is generated by imposing appropriate body forces in the 3D FEM elements that discretize the rock mass; thus, any initial stress state can be easily created. Midas GTS NX, which is the FEM code used in this paper, has a feature that allows the user to impose body forces into 3D finite elements by providing the components of the Cauchy stress tensor (σ xx , σ yy , σ zz , τ yz , τ xz , τ xy ). A similar feature exists in other FEM codes. The nodes at the boundaries are fixed, as illustrated in Figures 1c and 1d. Consequently, external forces are generated at these nodes (i.e. the reaction forces) that ensure equilibrium of the imposed geostatic stress field. Also, this boundary condition is reasonable since no displacements are expected far from the tunnel. Obviously, to achieve accurate numerical results, the model should be large enough, and the mesh properly refined.

3D FEM mesh
The tunnel investigated numerically was circular with 5m radius (r 0 ). The 3D FEM mesh had a cylindrical shape with 100r 0 diameter and 120r 0 length. The adopted size of the FEM mesh ensured accurate results even for highly nonlinear material (Vitali et al., 2018a). The mesh was refined near the tunnel face and was gradually coarsened towards the boundaries. The mesh refinement adopted follows the recommendations by Vitali et al. (2018a). 2 nd order hexahedron elements were adopted. The FEM mesh shown in Figure 1 has around 300,000 nodes and 76,000 elements. The results from the simulations were extracted from a refined mesh at the core of the model (Figure 1b). The length of the hexahedron elements at the core was 0.2r 0 in the axial direction, as recommended by Vitali et al. (2018a).
To ensure the accuracy of the numerical results, the refined mesh was extended to a distance of 6r 0 ahead the face and 12r 0 behind the face, as illustrated in Figure 1b. The results far-behind the face presented in this paper were extracted at a distance of 8r 0 behind the face, which is far enough from the face such that the 3D face effects are negligible and the results can be compared with the analytical solution.
The first phase of the analyses imposed the body forces in the 3D finite elements that discretize the rock mass, to generate the far-field stresses, while the displacements at the boundaries of the model were constrained. After this stage, the elements inside the tunnel were deactivated to simulate the tunnel excavation. Although the results presented in this paper were obtained with the FEM code Midas GTS NX, this modeling approach is general; thus, any FEM code that allows the user to impose body forces in the elements may be used. Note that the presented new approach is valid for elastoplastic rock masses, tunnels with support systems and for any geometry and construction sequence.

Tunnel in anisotropic rock and complex geostatic stress field
To verify the accuracy of the 3D FEM mesh with the proposed boundary conditions, as shown in Figure 1, the displacements and stresses at the tunnel perimeter were compared with those obtained with the analytical solution proposed by Vitali et al. (2020b). The transversely anisotropic elastic model was selected to represent the rock mass. The anisotropic rock properties are: Young modulus parallel to the rock structure, 2.67GPa, and perpendicular to the rock structure, 1.33GPa; shear modulus parallel to the rock structure, 1GPa, and perpendicular to the rock structure, 0.76GPa; Poisson's ratio parallel to the rock structure, 0.33, and perpendicular to the rock structure, 0.25. The dip angle was 64° and the strike direction, 37°. The tunnel was assumed aligned with the North (i.e. it is assumed that the z-axis is parallel to the North); thus, the axis of the tunnel is not parallel to any of the principal directions of material anisotropy. The strike direction is measured from the positive z-axis towards the positive x-axis (coordinate system shown in Figure 1). A highly complex geostatic stress field was selected. The farfield stresses with respect to the tunnel coordinate system ( Figure 1b) were σ xx,ff =7.5MPa; σ yy,ff =5MPa; σ zz,ff =7.5MPa; τ yz,ff =2.5MPa; τ xz,ff =2.5MPa and; τ xy,ff =-1.25MPa, where positive normal stresses denote compression. Using the proposed approach, such complex geostatic stress field was easily generated by imposing the body forces into the 3D elements. The nodes at the boundaries of the model were fixed; that is, the displacements at the boundaries were zero and the reaction forces balanced the (complex) geostatic stress field. Note that, in this scenario, the geostatic principal stress directions and the principal directions of material anisotropy are not aligned; thus, the use of traditional boundary conditions (Vitali et al., 2018a(Vitali et al., , b, 2020a could be problematic. Also, if traditional boundary conditions were used, the tunnel would not be aligned with the boundaries. Each tunnel orientation attempted would thus require a different mesh (Vitali et al., 2018b), which is a time-consuming task. Figure 2 presents the normalized stresses and displacements at the tunnel perimeter obtained with the 3D FEM model Figure 2. Comparison between numerical and analytical results: (a) tangential stress (σ θθ ) and tangential axial shear stress (τ θz ), normalized with respect to the vertical stress (σ v ), at the tunnel perimeter; and (b) radial (u r ) and axial (u z ) displacements, normalized with the tunnel radius (r 0 ), at the tunnel perimeter.
with the proposed boundary conditions ( Figure 1) and with the analytical solution (Vitali et al., 2020b). Positive radial displacements are towards the inside of the tunnel. The axial displacements sign convention is consistent with the z-axis direction (i.e. negative axial displacements are towards the excavated tunnel and positive towards the rock mass ahead of the face of the tunnel). As one can see, analytical and numerical results match, which shows that the proposed boundary conditions work well. The stresses and the displacements at the tunnel perimeter are not symmetric with respect to the horizonal and vertical axes as a consequence of the stress and material anisotropy. Anti-symmetric axial displacements and anti-symmetric tangential axial shear stresses (τ θz ) are induced around the tunnel perimeter. The tunnel cross-section is distorted in the axial direction about the axis of anti-symmetry at θ=126 o , as illustrated in Figure 2b. Note that the tangential axial shear stresses are larger at the locations where the axial displacements are zero and are smaller where the axial displacements are maximum (i.e. maximum axial displacement refers to the magnitude regardless of the direction). Also, the tangential stresses are larger where the radial displacements are minimum and are smaller where the radial displacements are maximum. The deformed tunnel cross-sections are shown in Figure 2b, that illustrates the axial distortion of the tunnel cross sections and the ellipsoidal shape of the deformed cross section.

Tunnel in vertically-foliated rock mass
The cases investigated by Vitali et al. (2020a) of tunnels in vertically-foliated rock masses were selected for further verification of the proposed method to impose initial/geostatic stress conditions. The cases were assessed using a 3D FEM model and the proposed method and boundary conditions ( Figure 1). The results were compared with those from the analytical solution presented by Vitali et al. (2020b) and with the 3D FEM model using traditional boundary conditions (Figure 3), from Vitali et al. (2020a). The same geostatic stress field selected by Vitali et al. (2020a) was adopted (i.e. major horizontal stress, σ H , of 10MPa; minor horizontal stress, σ h , of 5MPa and; vertical stress, σ v , of 5MPa). The tunnel was at an angle of 45° with the major horizontal stress. Two orientations of the major horizontal stress with respect to the rock structure were considered: major horizontal stress, σ H , perpendicular to the strike, and major horizontal stress, σ H , parallel to the strike. In these cases, the principal directions of stress and material anisotropy are aligned, which is ideal for the use of traditional boundary conditions (Vitali et al., 2020a); note that, when those directions are not aligned, the use of traditional boundary conditions may be difficult.
The far-field stresses and the rock structure are shown in Figure 4. The rock properties selected were the same as the previous case. With the tunnel angle at 45°, the far-field stresses σ xx,ff , σ yy,ff , and σ zz,ff , in the coordinate system attached to the tunnel, were: σ xx,ff = σ yy,ff =5MPa, σ zz,ff =7.5MPa, and τ xz,ff = ±2.5MPa, with the sign of the far-field shear stress being the only difference between the cases, as one can see in Figure 4. The figure also shows that the tunnel orientation, with respect to the rock structure, was the same in both cases. The mesh shown in Figure 1 was used for all the cases. This is an advantage of the proposed technique: the same 3D FEM mesh can be used to analyze the tunnel excavation under any geostatic stress state in any full anisotropic rock mass, which is not the case when using traditional boundary conditions. The FEM model with traditional boundary conditions is presented in Figure 3. This was the same 3D FEM model used by Vitali et al. (2020a). As one can see, the mesh was rather large, to prevent effects from the boundaries. All elements used were 2 nd order hexahedron elements, and the mesh refinement and the size of the model (Figure 3) were selected to ensure the accuracy of the results, following the recommendations provided by Vitali et al. (2018a). Figure 3b shows the refined mesh at the core of the model where the results were extracted. Figure 3c illustrates the plan view of the mesh with the boundary conditions. The rock structure was aligned with the sides of the discretization. Note that the tunnel was not aligned with the mesh, but at an angle Ψ=45°. The far-field stresses, as shown in Figure 3c, were applied to the boundaries of the discretization. At the faces of the mesh, opposite to where the stresses were applied, rollers were used. This was the first stage of the simulation, where far-field stresses were imposed, and all displacements were zeroed. That is, the geostatic stress conditions were imposed. In the second stage of the simulation, the elements of the tunnel were deactivated, without changing the boundary conditions imposed in the first stage.
The normalized stresses and displacements along the tunnel perimeter, far-behind the face, are presented in Figure 5. As one can see, the results obtained with different methods are the same (i.e. 3D FEM model with the proposed boundary conditions, Figure 1; analytical solution, Vitali et al. (2020b); and 3D FEM model with traditional boundary conditions, Figure 3). The consistency of the results with all three different methods indicates that the new approach is essentially correct. Anti-symmetric tangential axial shear stresses (τ θz ) and anti-symmetric axial displacements (u z ) were induced far-behind the face, as illustrated by the axially deformed tunnel cross-sections shown in Figures 5a.2 and 5b.2 (i.e. the tunnel cross-section is distorted about the vertical axis). Note that the direction of the axial distortion is not the same for the two cases investigated. The axial and the radial displacements are larger for the case where the major horizontal stress is perpendicular to the strike (Figure 5a.1) than for the case where it is parallel (Figure 5b.2). As discussed by Vitali et al. (2020a, b, c) and Vitali (2020), when the major horizontal stress is perpendicular to the strike, the axial distortion produced by the material anisotropy and by the stress anisotropy have the same tendency; thus, axial and radial displacements are increased. The opposite occurs when the major horizontal stress is parallel to the strike direction. Note that, as shown by Vitali et al. (2018bVitali et al. ( , 2019a, if the rock mass is isotropic and elastic and the tunnel is unsupported, the far-field shear stress does not affect the radial and the tangential displacements far-behind the face because the in-plane displacements do not depend on the axial stresses. However, if the rock mass is anisotropic or the tunnel is not aligned with the principal axes of material anisotropy, the far-field axial shear stresses (τ xz,ff and τ yz,ff ) induce displacements on the plane of the tunnel cross section (Vitali et al., 2020b). The far-field axial shear stress did not affect the tangential stresses (σ θθ ) at the tunnel perimeter, as one can see by comparing Figures 5a.1 and 5b.1. In contrast, the tangential axial shear stresses (τ θz ) were affected by the far-field shear stress. Figure 6 shows the normalized radial and axial displacements at the face of the tunnel for the two cases investigated. As one can see, the results with both 3D FEM models are the same (i.e. the 3D FEM model with the proposed boundary conditions, Figure 1, and the 3D FEM model with traditional boundary conditions, Figure 3). This is further evidence that the proposed method provides accurate results. As one can see in Figure 6, the displacements at the face are highly asymmetric. For the case where the major horizontal stress is perpendicular to the strike (Figure 6a), the tunnel cross section translates to the right (i.e. towards the positive x-axis) and, for the case where it is parallel (Figure 6b), to the left. The asymmetric displacements are more pronounced, and the axial displacements larger, when the major horizontal stress is perpendicular to the strike. Also, the location where the radial and axial displacements are maximum is the same in both cases analyzed, as well as the location where radial and axial displacements are minimum. For instance, for the case where the major horizontal stress is perpendicular to the strike (Figure 6a), the maximum axial and radial displacements occur at the  (1) tangential stresses (σ θθ ) and tangential axial shear stresses (τ θz ) normalized with respect to the vertical stress (σ v ); and (2) radial (u r ) and axial (u z ) displacements normalized with the tunnel radius (r 0 ) along the tunnel perimeter. (a) major horizontal stress perpendicular to the rock structure; and (b) major horizontal stress parallel to the rock structure. Figure 6. Normalized radial (u r ) and axial (u z ) displacements with respect to the tunnel radius at the face of the tunnel, for: (a) major horizontal stress perpendicular to the strike; and (b) major horizontal stress parallel to the strike. right springline, and the minimum, at the left springline. The opposite is observed when the major horizontal is parallel to the strike (Figure 6b). Thus, as one can see in Figure 6, the tunnel cross section translates towards the location where the axial displacement is smaller. Note that negative axial displacements are towards the excavated tunnel and positive, towards the rock mass ahead of the face of the tunnel. The anti-symmetric axial displacements are partially constrained at the face, which may explain the asymmetric deformations near the face. A detailed discussion on the influence of the stress and rock anisotropy on tunnel behavior is provided by Vitali et al. (2020a).

Conclusions
A new approach for 3D numerical modeling of tunnels in complex conditions is proposed. The geostatic stress field is generated by imposing body forces to the elements, while the boundaries of the model are fixed. The proposed approach is validated by comparing its results with those of a 3D FEM model where conventional boundary conditions are used (Vitali et al., 2020a), and with results from an analytical solution (Vitali et al., 2020b). The results from all three different methods are the same; thus, when the mesh is properly refined and the model sufficiently large, the numerical results obtained should be correct. The approach is well-suited for the design of tunnels under complex loading and/or ground properties.