Gradient-Enhancement of a Rock Model Applied to Numerical Simulations of Tunnel Advance

Rock is classiﬁed as a frictional cohesive material characterized by nonlinear stress-strain relations including softening leading to structural failure. To describe the mechanical behavior of rock comprising the accumulation of irreversible strains, strain hardening, and the degradation of stiffness and strength, an isotropic damage plasticity model for intact rock and rock mass was proposed in [1]. In the context of numerical simulations of tunneling, considering strain softening of rock allows to predict the formation of shear bands emanating from the tunnel surface due to the excavation process, which is crucial for engineers to forecast potentially dangerous situations related to failure. To ensure objective results in the softening regime upon mesh reﬁnement, in the aforementioned rock model the concept of the crack band approach by Bažant and Oh was adopted. As a remedy of known deﬁciencies related to this regularization approach, the rock model is extended by the over-nonlocal implicit gradient-enhancement proposed in [2]. Accordingly, a nonlocal ﬁeld of the damage-driving variable, implicitly deﬁned as the solution of a Helmholtz-like partial differential equation, is introduced yielding a fully coupled system of second-order PDEs. In the present contribution, the over-nonlocal gradient enhancement of the rock model is presented and applied together with a linear-elastic material model for shotcrete to ﬁnite element simulations of a deep tunnel constructed by the New Austrian Tunneling Method. Based on this real example, objectivity of the results with respect to the ﬁnite element discretization is discussed.

Rock is classified as a frictional cohesive material characterized by nonlinear stress-strain relations including softening leading to structural failure. To describe the mechanical behavior of rock comprising the accumulation of irreversible strains, strain hardening, and the degradation of stiffness and strength, an isotropic damage plasticity model for intact rock and rock mass was proposed in [1]. In the context of numerical simulations of tunneling, considering strain softening of rock allows to predict the formation of shear bands emanating from the tunnel surface due to the excavation process, which is crucial for engineers to forecast potentially dangerous situations related to failure. To ensure objective results in the softening regime upon mesh refinement, in the aforementioned rock model the concept of the crack band approach by Bažant and Oh was adopted. As a remedy of known deficiencies related to this regularization approach, the rock model is extended by the over-nonlocal implicit gradient-enhancement proposed in [2]. Accordingly, a nonlocal field of the damage-driving variable, implicitly defined as the solution of a Helmholtz-like partial differential equation, is introduced yielding a fully coupled system of second-order PDEs. In the present contribution, the over-nonlocal gradient enhancement of the rock model is presented and applied together with a linear-elastic material model for shotcrete to finite element simulations of a deep tunnel constructed by the New Austrian Tunneling Method. Based on this real example, objectivity of the results with respect to the finite element discretization is discussed.

Implicit Gradient-Enhancement of a Damage-Plasticity Rock Model
For modeling the highly nonlinear mechanical behavior of rock, an isotropic damage plasticity model for intact rock was proposed in [1], denoted henceforth as RDP model. It is formulated based on the combination of the plasticity theory and the continuum damage mechanics theory. It represents linear elastic material behavior, non-associated plastic flow, nonlinear isotropic hardening and softening material behavior, and nonlinear behavior in predominantly hydrostatic compression. For representing the material behavior of quasi-isotropic quasi-homogeneous rock mass, the RDP model incorporates empirical reduction factors to account for the influence of distributed discontinuities within the rock mass. Its stress-strain relation is expressed as σ = (1 − ω) C : (ε − ε p ), in which σ denotes the nominal Cauchy stress (force per total area), ω the scalar isotropic damage parameter ranging from 0 (undamaged material) to 1 (fully damaged material), ε the total strain and ε p the plastic strain. The nominal stress is related to the effective stress (force per undamaged area)σ by σ = (1 − ω)σ. The evolution of the damage parameter is defined by an exponential softening law ω(α d ) = 1 − exp(−α d /ε f ) with the softening modulus ε f and the strain-like softening variable α d . The latter depends on the evolution of the plastic volumetric strainε p,vol byα d =ε p,vol /x s (ε p,vol ), in which x s denotes a softening ductility measure accounting for the influence of multi-axial stress states on the softening behavior.
By considering softening material behavior in terms of a mere continuum formulation, in finite element simulations sensitive results with respect to the discretization are obtained. As a remedy, in the original version of the RDP model the mesh-adjusted softening modulus according to the crack-band theory was employed. It is a rather straightforward technique to overcome the pathological mesh sensitivity, however it has some severe deficiencies [4]. To eliminate some of these deficiences, in the following the implicit gradient-enhancement of the damage formulation is adopted for the RDP model. This approach forms a well-established regularization framework due to its computational efficiency and numerical stability. So far, several implicit gradient-enhanced damage and plasticity models or a combination of the latter were proposed, whereas in the present case the implicit gradient-enhanced approach proposed by Poh and Swaddiwudhipong [2] is pursued. Therein, nonlocality is introduced by replacing the strain-like softening variable α d in the evolution law of the damage parameter by a weighted softening variableα d leading to the modified softening law ω(α d ) = 1 − exp(−α d /ε f ). The weighted softening variable depends on the local softening variable α d and its nonlocal counterpartᾱ d byα d (α d ,ᾱ d ) = mᾱ d + (1 − m) α d with the weighting parameter m. This parameter should be chosen larger than 1 to achieve full regularization of the problem and hence this approach is often pronounced as the over-nonlocal formulation. Following the implicit definition of the nonlocal softening variableᾱ d , the latter is defined from a second order partial differential equationᾱ d − l 2 ∆ᾱ d = α d in Ω, with the length scale parameter l related to the radius of nonlocal interaction, the Laplace operator ∆ and Ω is the spatial domain occupied by the body of interest. As a consequence, the mechanical behavior of a point in the continuum does not only depend on its local state, but it is influenced by its neighborhood. However, nonlocality is introduced only for the damage regime and the plasticity part of the RDP model remains local. Together with the standard equilibrium equation a fully coupled set of equations is obtained describing the unknown displacement field u and the nonlocal softening variableᾱ d . Along the domain boundary, for the nonlocal softening variable the homogeneous Neumann boundary condition is assumed. For the solution of the system of equations, the finite element method is employed.

Numerical Simulation of Deep Tunnel Advance
The implicit gradient-enhanced RDP model is applied to numerical simulations of deep tunnel advance. Thereby, its performance to predict the formation of multiple shear bands emanating from the tunnel surface due to the construction process in a mesh-insensitive manner is assessed. The analyzed tunnel section is part of the Brenner Base Tunnel, which is located at a depth of 950 m below ground, is surrounded by Innsbruck quartz phyllite rock mass and is characterized by a circular profile with a diameter of 8.5 m. It was constructed by employing the drill, blast and secure procedure.
For the finite element simulations, the actually three-dimensional tunneling process is approximated by an equivalent twodimensional model assuming plane-strain conditions. Therein, full-face excavation of the tunnel profile and the installation of a shotcrete shell (thickness 0.2 m) are considered. Within the discretized domain of rock mass of 200 × 200 m, a hydrostatic geostatic stress state with a constant value of p (0) i =−25.7 MPa is assumed, according to the overburden. The influence of the time-dependent tunneling process on the analyzed tunnel profile is considered based on the convergence-confinement method: Initial equilibrium is established by applying the geostatic stress state and the internal pressure p i =p (0) i acting on the tunnel surface; Subsequently, the excavation procedure is simulated by gradually decreasing the internal pressure p i to zero following a step-wise reduction scheme, whereas the shotcrete shell is placed at a pressure reduction by 80 %. The employed material parameters for the surrounding rock mass of Innsbruck quartz phyllite are taken from [3]. The additional parameters required for the implicit gradient-enhancement are specified as ε f =7 × 10 −4 , l=50 mm and m=1.05. For the shotcrete shell, linear elastic material behavior is assumed with Young's modulus of 7000 MPa. The spatial domain of the rock mass and the shotcrete shell is discretized with 8-node quadrilateral finite elements employing full numerical integration. For interpolation of the displacement field and the nonlocal field, quadratic shape functions are employed. To trigger localization of strains into multiple shear bands around the tunnel profile, a small share of randomly distributed slightly weakened zones are introduced. The influence of the finite element discretization on the numerical results is studied by employing three different meshes, i.e., a coarse (300 mm element size in the vicinity of the tunnel), a medium (140 mm) and a fine mesh (70 mm).
The predicted evolution of the magnitudes of the displacement vectors at the rock-shotcrete interface over time are illustrated in Fig. 1, in terms of the mean values along the tunnel perimeter. Regarding the influence of the employed finite element mesh, a noticeable difference is recognized between the coarse and the medium mesh, however almost identical results are predicted by the medium and the fine mesh. In this case, the sufficiently resolved gradient of the nonlocal field is apparent. Fig. 2 shows the deformed tunnel structure after 72 h. The localization of displacements into shear bands is clearly visible, which indicates the transition from an initially quasi-continuous rock mass to a quasi-discontinuous rock mass. Again, for the medium and the fine mesh almost identical deformations are obtained, only the distribution in the upper left corner slightly deviates. Despite the complex failure mode accompanied by the formation of multiple shear bands, the capability of the regularization by means of the gradient-enhanced approach is demonstrated.