Evaluation of the Implicit Gradient-Enhanced Regularization of a Damage-Plasticity Rock Model

In the present publication, the performance of an implicit gradient-enhanced damage-plasticity model is evaluated with special focus on the prediction of complex failure modes such as shear failure. Hence, it complements studies on predominant mode I failure frequently found in the literature. To this end, an implicit gradient-enhanced damage-plasticity rock model is presented and validated by means of 2D and 3D finite element simulations of both laboratory tests on intact rock specimens as well as a large-scale structural benchmark related to failure of rock mass. Thereby, a wide range of loading conditions comprising unconfined and/or confined, tensile and/or compressive stress states is considered. The capability of the gradient-enhanced rock model for representing the mechanical response objectively with respect to the finite element discretization and realistically compared to measurement data is assessed. It is shown that complex failure modes and the respective load–displacement curves are predicted in a mesh-insensitive manner.


Introduction
The appropriate representation of complex failure modes of cohesive-frictional materials in numerical simulations is of great interest for many engineering applications, which include, among others, geotechnical applications characterized by material failure playing a dominant role in the overall structural response [1,2].A prominent example of material failure under highly confined stress states is given in the context of tunnel construction: During the construction of deep tunnels, high geostatic stresses in the surrounding rock mass are redistributed in consequence of the excavation process [3,4].Depending on the quality of the surrounding rock mass, the installed supporting measures, and the type of excavation process, those stress redistributions can lead to damage in the rock mass, often accompanied by the transition from the rock mass as a continuum to rock blocks moving towards the tunnel center with localized shear bands indicating the sliding interfaces.It follows that an adequate representation of such failure phenomena by means of advanced constitutive models is of great importance and a prerequisite for the risk assessment of potential structural collapse.
Constitutive models based on the combination of the theory of plasticity and continuum damage mechanics, simply denoted as damage-plasticity models, provide a powerful framework for describing inelastic deformations, hardening and softening material behavior, as well as stiffness degradation due to damage.They are well suited for the description of cohesive-frictional materials such as concrete, rock, and soils, since different types of material failure, e.g., cracking in tension, crushing in compression, or failure under mixed stress states, can be described in a realistic manner [5].On the basis of damage-plasticity models, the material behavior is described in terms of mere continuum relations, i.e., constitutive relations describing the behavior of an infinitesimal material point.However, softening material behavior, described in terms of continuum models, exhibits several theoretical deficiencies, as reported and summarized in [6]: • the softening process zone is infinitesimally small; • at structural level, snapback behavior due to the infinitesimally small softening zone is observed; • the dissipated energy during the failure process is zero due to the infinitesimal zone in which energy is dissipated.
In a mathematical sense, those issues are related to the loss of ellipticity of the initial boundary value problem (IBVP) in static and quasi-static analyses.In finite element analyses, those deficiencies lead to mesh-dependent results, commonly referred to as a pathological mesh sensitivity.Failure patterns like cracks tend to localize into the smallest possible bandwidth, i.e., usually a single layer of finite elements, and upon mesh refinement the localization zone is decreasing to an arbitrarily small domain.Accordingly, the obtained results are not objective with respect to the finite element mesh, i.e., the results are sensitive with respect to the numerical discretization scheme.In the past decades, several techniques to overcome these deficiencies in numerical simulations have been proposed, for instance the crack band approach based on a mesh-adjusted softening modulus [7], nonlocal approaches of the integral type [8], models based on the Cosserat continuum [9], viscoplastic formulations [10], phase field models [11], and explicit and implicit gradient-enhanced formulations [12].
Among these approaches, models based on implicit gradient-enhanced formulations by now form a well-established branch in the literature due to their computational efficiency and numerical stability [12].Numerous implicit gradient-enhanced damage and plasticity models have been proposed in recent years [13][14][15][16][17][18][19][20][21], and their performance has been demonstrated based on different examples of material failure.In particular, many gradient-enhanced models have been developed explicitly for concrete [13,16,17,21], and special attention has been paid to the proper representation of cracking under predominantly tensile stress states [17,22].However, while cracking in tension (cf. Figure 1 left) is an important failure mode for concrete structures, considerably less attention was paid to shear failure (cf. Figure 1 right), i.e., pure mode II failure, or mixed failure under confined stress states.In fact, many of the available models have been validated based on examples of mode I failure, and often it is tacitly assumed that they perform equally well for more complex failure modes.The apparent gap in the literature on the assessment of gradient-enhanced models for describing such complex failure modes motivates a systematic investigation of a gradient-enhanced constitutive model applied to different types of material failure.To this end, in the present contribution, a gradient-enhanced damage-plasticity model for rock is proposed and evaluated.This model is considered as a representative for a wide class of gradient-enhanced damage-plasticity models.The model is based on the damage-plasticity model by Unteregger et al. [24], and its damage formulation is extended following the implicit gradient-enhanced approach proposed by Poh and Swaddiwudhipong [17].Based on numerical examples involving complex failure modes, the capability of the model to capture different types of material failure in a realistic and objective manner in finite element simulations will be demonstrated.
The remainder of the paper is organized as follows: In Section 2, the original damage-plasticity model for intact rock and rock mass, proposed in [24,25], is briefly summarized.In Section 3, the implicit gradient-enhancement of the rock model is presented, and the numerical implementation into a finite element framework is discussed.Section 4 covers numerical simulations of laboratory tests including wedge splitting tests as well as triaxial compression tests and triaxial extension tests with various levels of confining pressure.Additionally, a benchmark example of tunnel excavation will be presented, and the representation of complex failure modes in an objective manner with respect to the employed finite element mesh will be demonstrated.Finally, in Section 5, the paper is closed with a summary and a discussion on recommended future research activities.

Damage-Plasticity Model for Intact Rock and Rock Mass
The damage-plasticity model for intact rock and rock mass, denoted as RDP model in the following, was proposed originally for intact rock in [24] and was further extended to rock mass in [25] by incorporating empirical down-scaling factors to account for the influence of discontinuities according to [26][27][28].
The RDP model is based on the theory of plasticity formulated in the effective stress space combined with continuum damage mechanics.The stress-strain relation is expressed as in which σ describes the nominal Cauchy stress tensor (force per total area), ω the scalar isotropic damage parameter ranging from 0 (undamaged material) to 1 (fully damaged material), C the fourth order elastic stiffness tensor, ε the total strain tensor, and ε p the plastic strain tensor.The effective stress tensor σ (force per undamaged area) is linked to the nominal stress tensor by The elastic domain is bounded by the smooth Hoek-Brown yield criterion [29,30] formulated in the Haigh-Westergaard coordinates of the effective stress tensor, i.e., the mean stress σm , the deviatoric radius ρ, and the Lode angle in the deviatoric plane θ.In addition, a stress-like hardening variable q h (α p ) is incorporated leading to the definition of the yield function f p as Therein, f cu is the uniaxial compressive strength, m 0 is the friction parameter, r(θ, e) is the Willam-Warnke function to describe the shape of the yield surface in deviatoric planes, and parameters m b /m 0 and s are empirical down-scaling factors to account for discontinuities in rock mass, the latter depending on the geological strength index GSI and the disturbance factor D according to [26].A default value for the eccentricity e of 0.51 is proposed in [24].For representing material behavior of intact rock, m b /m 0 and s are equal to 1.
The flow rule for describing the evolution of plastic strains is defined in the effective stress space as with λ denoting the consistency parameter and g p ( σ, q h (α p )) the non-associated plastic potential function expressed as Therein, volumetric plastic flow is controlled by dilatancy parameters m g1,rm and m g2,rm .In the expression for m g1,rm = (m b /m 0 ) m g1 , m g1 is calibrated from experimental results (uniaxial tension, uniaxial compression, and triaxial compression tests) of intact rock specimens and m g2,rm is determined from m g2,rm = 2 m g1,rm − 6 f tu / f cu (6) such that the lateral plastic strain rate in uniaxial tension is zero.Uniaxial tensile strength f tu is calculated from (3) as Hardening material behavior is described by means of the stress-like internal variable which is conjugate to the strain-like hardening variable α p and contains the yield stress in uniaxial compression f cy .The evolution law for the strain-like hardening variable α p is given as in which ∂g p ( σ, q h (α p ))/∂ σ is the norm of the gradient of the plastic potential function with respect to effective stress, E rm /E i denotes the reduction of the Young's modulus of rock mass compared to intact rock [26] ranging from 0 (completely disintegrated) to 1 (intact rock), and x h ( σm ) is a measure for describing hardening ductility, defined as In (10), model parameters A h , B h , C h , D h , and G h control the hardening behavior.They are calibrated by experimental data from uniaxial tension, uniaxial compression, and triaxial compression tests.In [24,25], default values of B h = 10 −5 , D h = 10 −6 and G h = 0 are suggested in absence of respective experimental data.Damage is provoked when the hardening variable q h (α p ) attains its maximum value of 1.At this stage, the scalar isotropic damage parameter ω starts evolving dependent on the strain-like internal softening variable α d .This relation is described by means of an exponential softening law as with the softening modulus ε f controlling the slope of the softening curve.The rate of the strain-like internal softening variable α d is computed from the volumetric part of the plastic strain rate εp,vol = tr( εp Therein, the softening ductility measure accounts for the influence of multi-axial stress states on the softening behavior, with εp,vol = tr( − εp ) describing the compressive part of the volumetric plastic strain rate.Model parameters A s and B s are calibrated from experimental data of uniaxial compression and triaxial compression tests.Again, in absence of respective experimental data for parameters A s and B s , default values of A s = 15 and B s = 2 are proposed in [24].

Implicit Gradient-Enhancement of the Damage-Plasticity Model for Intact Rock and Rock Mass
Softening material behavior leads to an ill-posed initial boundary value problem and consequently to pathological mesh-sensitivity in finite element simulations.As a remedy, in [24], the crack band approach was employed for the RDP model.While the crack band approach is a rather simple regularization technique, it is also characterized by a number of shortcomings, which are addressed in [31].Motivated by those deficiencies, in the following, a more sophisticated regularization technique based on the implicit gradient-enhanced formulation [32] is presented.To this end, the approach by Poh and Swaddiwudhipong [17] proposed for a damage-plasticity model for concrete and based on the gradient of the internal softening variable, is adopted.By incorporating the gradient of an internal variable into the constitutive relations, nonlocality is introduced.Thus, the mechanical response of a material point does not exclusively depend on its local state, but is also influenced by the state in its neighborhood.
Nonlocality is incorporated by replacing the local softening variable α d by a weighted softening variable αd in the exponential damage law (11), which is expressed as Therein, the softening modulus ε f is a material parameter.The weighted softening variable αd is calculated from a combination of the local softening variable α d and its nonlocal counterpart ᾱd as αd = m ᾱd + (1 − m) α d (15) in which m denotes a weighting parameter.Choosing m larger than 1 yields the over-nonlocal formulation [33] to achieve full regularization of the problem, as proven in [34].Furthermore, by ensuring that αd can only increase, damage is considered as an irreversible process [35].
According to the implicit approach, the field of the nonlocal softening variable ᾱd , henceforth simply denoted as the nonlocal field, is defined implicitly as the solution of a higher-order partial differential equation.Adopting the formulation by Poh and Swaddiwudhipong [17], for the description of the nonlocal field a second order partial differential equation is employed as in which l denotes a length scale parameter defining the radius of nonlocal interaction, ∆ is the Laplace operator and α d represents the local softening variable of ( 12), and Ω is the spatial domain occupied by the body under consideration.Equation ( 16) is a screened-Poisson equation, commonly denoted as the Helmholtz equation in the context of implicit gradient-enhanced formulations [12,36].It is apparent that nonlocality affects only the damage part of the model and the plasticity part of the RDP model remains local.
As suggested in [32,37], homogeneous Neumann boundary conditions are assumed as ∇ᾱ d • n = 0 on the entire boundary Γ of the domain with the normal vector to the boundary n.This boundary condition was interpreted in [38] in the context of phase-field models enforcing cracks to occur perpendicular to the boundary.
The set of governing equations is completed by the equilibrium equation in which f is the vector of body forces.For the boundary conditions, the surface traction vector t = σ • n on Γ t and prescribed displacements u = ū on Γ u are assumed.Partial differential Equations ( 16) and ( 17) form a fully coupled system with the unknown displacement vector u and the nonlocal field ᾱd , which is solved by means of the finite element method.To this end, the weak form is formulated, which is subsequently discretized in space and in time.To obtain the weak form, the set of partial differential equations is multiplied by test functions δu for the displacement field and δ ᾱd for the nonlocal field.Integration over the domain, application of the divergence theorem, incorporation of the infinitesimal strain ε, and consideration of δu = 0 on Γ u yields the weak form of the IBVP expressed in Voigt notation as The displacement field u and the nonlocal field ᾱd are approximated over the domain using a Bubnov-Galerkin approach as in which N (•) contains the shape functions and q (•) are column vectors of the nodal unknown parameters, both expressed in the global form employing the standard assembly procedure, with (•) standing for the displacement field u and the nonlocal field ᾱd , respectively.The infinitesimal strain ε and the gradient of the nonlocal field ∇ᾱ d are discretized as with the strain-displacement matrix B u and the row vector of the spatial derivatives of the shape functions for the field of the nonlocal softening variable B ᾱd , again both expressed in the global form employing the standard assembly procedure.It follows that the shape functions for both fields must meet the requirement of C 0 -continuity.Due to the nonlinear and path-dependent character of the IBVP, an incremental solution procedure is employed.For the incremental solution procedure, a discrete (pseudo-)time interval [t (n−1) , t (n) ] is considered such that the body under consideration is in equilibrium at time t (n−1) .At this time, the nodal values q (n−1) u and q (n− 1)  ᾱd , the stress σ (n−1) and the internal variables α (n−1) p and α (n−1) d are known.An incremental load is applied such that the traction vector and the vector of body forces at time t The updated variables of the constitutive relations at time t (n) , i.e., σ (n) , α (n)  p , α (n) d , are evaluated by means of a stress-update algorithm, employing an implicit integration scheme following the return-mapping approach [39] for the plastic regime and a subsequent explicit evaluation of the damage part.Finally, the incremental discretized weak form is obtained as Ω δq (n)   ᾱd Since Equations ( 24) and ( 25) must hold for arbitrary kinematically admissible test functions δq (n)   u and δq (n)  ᾱd , they can be recast into the residual format as with f u int and f u,(n) ext denoting internal and external force vectors related to the displacement field, and R ᾱd is the residual vector for the nonlocal field.Due to the nonlinear dependence of the system of Equations ( 26) and ( 27) on the unknown nodal solution vector an iterative Newton-Raphson solution procedure is employed.The nodal unknowns at time t (n)  in the i-th iteration step are composed of q (n,i) = q (n−1) + ∆q (n,i) , where ∆q (n,i) has to be determined.Linearization of Equations ( 26) and ( 27) at the state q (n,i−1) with the initial guess of the nodal unknowns q (n,0) = q (n−1) yields the iterative procedure for the correction of the nodal unknowns ∆∆q (n,i) for time t (n) after the i-th iteration step ᾱd ᾱd with the submatrices of the system matrix given as in which ∂σ (n) /∂ε (n) , ∂σ (n) /∂ ᾱ(n) d , and ∂α (n) d /∂ε (n) are the consistent tangent stiffness submatrices.They represent the derivatives of the constitutive relations consistent with the numerical algorithm for integrating the path-dependent rate constitutive equations, which is essential for the full Newton-Raphson scheme in order to preserve a quadratic rate of convergence.Due to the non-associated plastic flow rule of the RDP model and the coupling of the displacement field and the nonlocal field, the system matrix is unsymmetric.From the computed correction of the nodal unknowns ∆∆q (n,i) , the updated solutions of the incremental nodal unknowns, ∆q (n,i) = ∆q (n,i−1) + ∆∆q (n,i) and the total nodal unknowns q (n,i) = q (n−1) + ∆q (n,i) are obtained.

Numerical Study
The aim of the present numerical study is to evaluate the performance of the implicit gradient-enhanced RDP model for predicting the mechanical response of structures, involving the softening behavior of rock in a realistic and mesh-insensitive manner for a wide range of loading conditions.To this end, a numerical study is presented, related to both laboratory tests and practical applications.It consists of the following parts: 1. 2D modeling of mode I failure in wedge splitting tests on Indiana limestone performed by Brühwiler and Saouma [40], 2. 3D modeling of shear failure in triaxial compression tests performed by Blümel [41] on specimens of Innsbruck quartz phyllite, considering the influence of confined stress states attaining the compressive meridian of the yield surface, 3. 3D numerical simulations of triaxial extension tests performed on the same type of specimens, investigating the influence of confined stress states attaining the tensile meridian of the yield surface, and 4. 2D numerical simulations of the excavation of a deep tunnel in Innsbruck quartz phyllite rock mass leading to the formation of shear bands in the vicinity of the tunnel for demonstrating the capability of the gradient-enhanced RDP model to predict failure of a complex structure.

Modeling of Mode I Failure, Demonstrated by Analyzing Wedge Splitting Tests on Indiana Limestone
In a first step, the capability of the gradient-enhanced RDP model for predicting mode I failure is assessed.To this end, the experimental study of wedge splitting tests on Indiana limestone performed by Brühwiler and Saouma [40] is considered.The investigated specimen is illustrated in Figure 2.During the experimental test, the splitting force was applied by a vertically driven wedge, exerting a pressure against roller bearings that were mounted on both sides of the groove in the specimen.The vertical (machine) force and the crack mouth opening displacement (CMOD) were recorded, and the latter was controlled during the experimental test.The energy conjugate force to the CMOD, the splitting force F sp , was calculated from the vertical force considering the geometry of the wedge and neglecting any frictional effects.In total, 5 tests were performed, but in [40] detailed load-displacement curves were presented only for Test 3. To investigate the degradation of the stiffness of the specimen during crack propagation, several loading/unloading cycles were performed.
The experimental test is simulated by means of a two-dimensional finite element model assuming plane stress conditions.According to Figure 2, the specimen is supported in the vertical direction at the bottom center over a width of 10 mm and in the horizontal direction at midpoint.The splitting force F sp , transmitted by the wedge, is approximated by the pressure p sp acting on the lateral groove faces.The simulation is performed in a displacement-driven manner by controlling the CMOD (i.e., the relative horizontal displacement between points (a) and (b) in Figure 2) during the loading and unloading cycles.The material parameters for representing the Indiana limestone were determined by a best fit with the recorded experimental results.In the present example of mode I failure, only few parameters have a significant influence on the results: E = 22000 MPa, ν = 0.15, f cu = 20 MPa, f cy = 2/3 f cu = 13.33 MPa, m 0 = 6.5, m g1 = 5, A h = 5 × 10 −3 , C h = 20, ε f = 8 × 10 −4 , m = 1.05, and l = 4 mm.For the weighting parameter m, any value larger than 1 is sufficient to ensure full regularization of the problem by avoiding spurious localization of plastic strain [33].A discussion on the influence of the weighting parameter m may be found in [42] for an integral-type nonlocal model, where, however, for parameter m = 2, an overestimation of the energy dissipation close to a notch was demonstrated.This non-physical effect was also observed by the authors with increasing influence for larger values of m.Thus, parameter m is chosen just slightly larger than 1, in accordance with proposed values from the literature [17,43].For the remaining model parameters e, B h , D h , G h , A s , and B s , the default values summarized in Section 2 are employed.From Equation ( 7), the uniaxial tensile strength is calculated as f tu = 3 MPa.
To investigate the influence of the finite element discretization on the predicted results, different meshes are employed: Three structured meshes with fully integrated 4-node quadrilateral elements and element sizes of 5 mm, 1 mm, and 0.5 mm in the vicinity of the expected crack path, as well as one unstructured mesh with the same element type and an element size of 1 mm along the crack path.Accordingly, the element size of 5 mm is slightly larger compared to the assumed length scale l of 4 mm for the coarse mesh and considerably smaller for the medium and fine mesh.The purpose of the unstructured mesh is to investigate a potential bias of the crack pattern by following the grid lines of the finite element mesh, since mesh-biased crack paths have been observed for smeared crack models based on the crack band approach [44].
In Figure 3, the resulting load-displacement curves, i.e., splitting force F sp versus CMOD, are shown for the considered experimental test and the numerical simulations.It can be concluded that the qualitative shape of the experimentally obtained curve is approximated quite well in the numerical simulations.Regarding the influence of the finite element mesh on the computed load-displacement curves, a slight difference between the predicted response based on the structured coarse mesh and the one based on the structured medium mesh is visible.This is explained by the rough approximation of the gradient of the nonlocal field by the coarse mesh.In contrast, almost identical load-displacement curves are obtained for the structured medium and the structured fine mesh, confirming that the gradient of the nonlocal field can be resolved sufficiently by those meshes.The unstructured mesh results in a similar load-displacement curve, demonstrating that the gradient-enhanced approach is not biased by the orientation of the finite element mesh.
The computed deformation of the specimen is depicted in Figure 4 for the three structured meshes.By increasing the splitting force during the numerical simulations, the tensile strength of the material is attained at first in the elements directly below the notch; consequently, damage is initiated.This leads to large, localized deformations in this region.Upon further loading, the increase of these deformations reflects the opening of a crack, propagating towards the bottom of the specimen.This is represented by the gradient-enhanced damage-plasticity approach in a smeared manner.The width of the damage zone is related to the length scale parameter l (cf.( 16)).Furthermore, it can be seen that an identical symmetrical response with respect to the vertical axis of symmetry is obtained for all three structured meshes.In Figure 5, the distribution of the damage variable ω computed for the three structured meshes is plotted at a CMOD = 0.3 mm.It is visible that damage is also accumulated slightly above the notch tip.This is a consequence of the diffusive character of the gradient-enhanced formulation.Furthermore, the present gradient-enhanced RDP model with the constant length scale l predicts a rather broad zone of complete damage (red region in Figure 5).In fact, this behavior does not represent damage localizing into a discrete macrocrack, and has been addressed in [45,46].However, comparison of the predicted damage distributions for the different meshes demonstrates mesh-insensitivity of the obtained failure patterns.Figures 3 and 6 show an identical load-displacement curve and identical deformation pattern and damage distribution computed by means of the unstructured mesh, which underlines the capability of the gradient approach to produce mesh-insensitive results.

Modeling of Shear Failure, Demonstrated by Analyzing Triaxial Compression Tests on Innsbruck Quartz Phyllite
In a second step, the performance of the gradient-enhanced RDP model for predicting shear failure under confined stress states attaining the compressive meridian of the yield surface is assessed.To this end, numerical simulations of triaxial compression tests on specimens of Innsbruck quartz phyllite performed by Blümel [41] are conducted.
For the experimental program, intact rock specimens with a diameter of 35 mm and a height of 70 mm were taken from drill cores sampled at the construction site of the Brenner Base Tunnel.A series of triaxial compression tests with different levels of confining pressures p (0) = 0 MPa (uniaxial compression), 12.5 MPa, 25.0 MPa, and 37.5 MPa was conducted.The experiments were performed in a sequential manner: Firstly, the confining pressure was applied, resulting in an initial hydrostatic stress state in the specimen, and subsequently, an axial pressure was applied displacement-driven up to failure.The mechanical response in the post-peak regime was also recorded.
Finite element simulations of triaxial compression tests are often performed in an approximate manner.Commonly, such tests are simplified as single element tests in which the non-homogeneous deformation of the specimen observed in the experiments cannot be captured, or by two-dimensional models in which three-dimensional effects are neglected, e.g., [24,[47][48][49].By contrast, in the present study, a three-dimensional finite element model is employed in order to capture the three-dimensional deformations in a realistic manner.The numerical model with the prescribed boundary conditions and loads is illustrated in Figure 7. Since the failure mode is expected to be symmetric with respect to a vertical plane through the center axis of the specimen, symmetry is exploited by considering only one half of the specimen.By analogy to the experimental tests, the numerical simulations are performed in two sequential steps: firstly, the specimen is supported vertically at the bottom face and the confining pressure is applied; secondly, the axial loading is applied by imposing a uniform vertical displacement at the top of the specimen.
The specimens are discretized with 20-node hexahedral elements employing reduced numerical integration.For both the displacement field and the nonlocal field, quadratic shape functions are used.To analyze potential mesh-sensitivity of the gradient-enhanced RDP model, three different structured finite element meshes are examined.A coarse mesh with 828 elements (an element size of 4 mm), a medium mesh with 5950 elements (an element size of 2 mm), and a fine mesh with 13,160 elements (element size of 1.5 mm) are employed.To trigger localized failure, at the center of the specimen, a small zone of slightly weakened elements is introduced, as indicated in Figure 7.In the numerical simulations, snapback behavior may occur, i.e., a simultaneous decrease of the load and the displacement after attaining the peak load.At this stage, displacement-controlled experiments become unstable.To overcome potential snapback behavior in the numerical simulations, the indirect displacement control technique [50,51] is employed.For the present example, this technique can be applied by enforcing a monotonic decrease of the vertical distance between two nodes, with each node located at one boundary of the expected shear band.The material parameters required for the elastic-plastic part of the RDP model were identified from single element simulations, as discussed in [52].The additional parameters for the softening regime A s , ε f , l, and m are calibrated from the present numerical simulations for a best fit with the experimental results for the confining pressure of p (0) = 37.5 MPa.The employed parameters for Innsbruck quartz phyllite are E = 56670 MPa, ν = 0.2, f cu = 42 MPa, f cy = 29.5 MPa, m 0 = 12.0, m g1 = 9.9, A h = 0.0045, C h = 8.8,A s = 4, ε f = 4 × 10 −4 , m = 1.05, l = 2 mm.For the remaining model parameters e, B h , D h , G h , and B s , the default values summarized in Section 2 are employed.The experimental results for confining pressures of p (0) = 0 MPa, 12.5 MPa, and 25 MPa, which have not been used for calibration, serve for validation of the numerical model.
Figure 8 shows the load-displacement curves obtained from the experiments and the numerical simulations for the different confining pressures.Note that the non-zero axial force at the beginning is the consequence of the initial hydrostatic stress state due to the applied confining pressure.Expectedly, for p (0) = 37.5 MPa and for the uniaxial compression test (p (0) = 0 MPa), the experimentally obtained peak load is represented very well since the test results were used for calibration.For p (0) = 12.5 MPa and 25 MPa, the peak loads are predicted satisfactorily, slightly underestimating the experimental results.The qualitative shape of the softening branch is also represented quite well.For the uniaxial compression test, no meaningful experimental results after the peak stress were recorded during the experiments.This unstable behavior is also manifested in the numerical simulations, for which strong snapback behavior is observed.The numerical results computed on the basis of the different finite element meshes reveal the capability of the gradient-enhanced RDP model to regularize the underlying IBVP: Once the mesh is sufficiently fine, mesh-insensitive load-displacement curves are obtained.
Figure 9a shows the deformed specimen with the confining pressure p (0) = 37.5 MPa in the final stage of the triaxial compression tests, computed by means of the three different meshes.While the displacements in the lower and the upper part of the specimen are almost uniform, the displacements localize into a single inclined shear band in the center part of the specimen.In the experiments, localization into a shear band was found to be the dominant failure mode, as shown in Figure 9b for a specimen tested with p (0) = 37.5 MPa.The formation of the shear band is also reflected by the distribution of the damage variable ω shown in Figure 10.Comparing the predicted damage distributions for the medium and the fine mesh confirms this distribution as insensitive with respect to the discretization.The influence of the level of confining pressure on the inclination of the shearing zone is demonstrated in Figure 11 based on the predicted distribution of the vertical displacement for p (0) = 0 MPa, 12.5 MPa, 25 MPa, and 37.5 MPa, respectively.For the uniaxial compression test (p (0) = 0 MPa), the zone of localized displacements is strongly inclined.With increasing confining pressure, the inclination angle of the shear band is decreasing, which is best visible by comparing the results for p (0) = 0 MPa and p (0) = 12.5 MPa.Upon further increase of the confining pressure, the inclination of the shear band becomes slightly smaller.The represented dependence of the inclination angle on the level of confining pressure is explained by the curvature of the compressive meridian of the yield surface and of the employed plastic potential function of the RDP model.

Modeling of Shear Failure, Demonstrated by Analyzing Triaxial Extension Tests on Innsbruck Quartz Phyllite
In a third step, the ability of the model for predicting shear failure under confined stress states attaining the tensile meridian of the yield surface is demonstrated.To this end, numerical simulations of triaxial extension tests on specimens with geometric and material properties identical to those of the previously presented triaxial compression tests are performed.Similar to the triaxial compression tests, confining pressures p (0) of 0 MPa, 12.5 MPa, 25 MPa, and 37.5 MPa are investigated.The same finite element meshes are employed for investigating the influence of the discretization.In contrast to the triaxial compression tests, subsequent to the application of the initial hydrostatic stress state, generated by the confining pressure, a displacement in positive vertical direction is applied at the top surface of the specimens.Since for triaxial extension tests experimental results are not available, the present study focuses on the assessment of the influence of the finite element mesh on the obtained results and serves as further verification of the gradient-enhanced RDP model.
The predicted load-displacement curves are shown in Figure 12.While nearly identical results are obtained for the medium mesh and the fine mesh, the load-displacement curves for the coarse mesh are somewhat more ductile.This discrepancy indicates the coarse mesh as insufficiently fine for the accurate resolution of the gradient of the nonlocal field.Compared to the results from the triaxial compression tests, a more brittle structural response is predicted, resulting in strong snapback behavior in the post-peak regime for all three confining pressures and, in particular, for the uniaxial tension test.In contrast to the triaxial compression tests, the structural response becomes increasingly brittle as confining pressure increases, and the peak load decreases gradually.This phenomenon was also observed in experiments on Berea sandstone in [53], and it is characteristic of brittle and quasi-brittle materials.
Figure 13 shows the computed deformation of the specimen for an applied axial displacement of 0.1 mm and a confining pressure of 37.5 MPa.Compared to the triaxial compression tests, a smaller inclination angle of the shear band is predicted for the medium mesh and the fine mesh, whereas for the coarse mesh a considerably steeper inclination angle of the shear band is obtained.Again, this discrepancy indicates the insufficient representation of the gradient of the nonlocal field by the coarse mesh.

Numerical Simulation of Localizing Deformations in Deep Tunnel Excavation
Finally, the performance of the gradient-enhanced RDP model for predicting the formation of multiple shear bands during the excavation of a deep tunnel is assessed.This benchmark is derived from a stretch of the Brenner Base Tunnel constructed by the drill, blast, and secure procedure, which has already been the subject of investigations in previous publications [52,54].An analogy to the present problem can be found in the context of petrol engineering, where the formation of shear bands has been observed and reported for borehole breakout [47,55,56].
Since in this contribution the major focus is on the assessment of the gradient-enhanced rock model, the tunnel excavation is approximated by means of a simplified two-dimensional model.Supporting measures like a shotcrete shell or rock anchors are neglected.Potential time-dependent effects of the mechanical behavior of rock mass due to the excavation procedure are not considered.For the excavation of the tunnel profile by means of drill and blast, either a full-face or a sequential excavation procedure can be employed.The chosen excavation sequence may have a considerable impact on the stability of the tunnel.In this numerical model, the worst case scenario is considered by assuming full-face excavation of the circular tunnel profile without any supporting measures, which results in the maximum loading of the rock mass in the vicinity of the tunnel.
The IBVP of tunnel excavation with its geometry, initial conditions, and boundary conditions is illustrated in Figure 14.Within the discretized domain of rock mass, a hydrostatic geostatic stress state characterized by a pressure of p (0) i = 25.7 MPa is assumed, corresponding to the overburden at the tunnel axis of 950 m.In the numerical simulations, initial equilibrium is established by applying the geostatic stress state together with the internal pressure p (0)  i acting on the excavation boundary.
The excavation procedure is simulated by gradually decreasing the internal pressure p (0) i to zero.Since the analyzed tunnel section is located in Innsbruck quartz phyllite rock mass, most of the material parameters of Section 4.2 are adopted.The material behavior of rock mass in contrast to intact rock is considered by empirical down-scaling factors based on the geological strength index GSI and the disturbance factor D proposed by Hoek and Brown [27] accounting for the influence of distributed discontinuities.They are taken from the geological survey, reported in [52].The additional material parameters for Innsbruck quartz phyllite rock mass are GSI = 40, D = 0, A s = 15, ε f = 7 × 10 −4 , and l = 50 mm.To trigger the formation of shear bands in spite of the axisymmetric problem, non-uniformly distributed rock mass properties are employed by introducing zones in which the strength of the rock mass is slightly weakened (indicated in Figure 14).It was verified that these weakened zones do not affect the predicted mechanical response before the onset of strain softening.For investigating the influence of the finite element mesh on the predicted results, three structured meshes with fully integrated and reduced integrated 8-node quadrilateral elements with element sizes of 300 mm (8078 elements), 140 mm ( 21 Upon decreasing the internal pressure, initially the rock mass behavior remains in the linear elastic regime, followed by the formation of plastic zones emerging from the excavation boundary.Once the strength of the rock mass is attained, strain softening is initiated.From the onset of strain softening, the strains are localizing into narrow zones of the rock mass, and, eventually, large displacements accumulate.Finally, at a certain release level of the internal pressure, equilibrium is lost.
Figure 15 shows the deformed rock mass at the level of 12% of the internal pressure for the three meshes with fully integrated elements.For the medium and the fine mesh, the localization of displacements into narrow zones has already reached a very progressed stage close to failure, and the formation of shear bands is clearly visible.The non-axisymmetric displacement field indicates the transition from an initially quasi-continuous rock mass to quasi-discontinuous rock mass.The latter is characterized by blocks of rock mass moving towards the tunnel center, so potential failure of the tunnel is imminent.Concerning the influence of the finite element discretization, for the medium and the fine mesh, an almost identical displacement field is obtained.Slight differences between those two solutions can be observed only for the upper right quadrant.Figure 16 shows the load-displacement curves, i.e., the normalized internal pressure versus the mean displacement magnitude along the excavation boundary computed for each mesh.The load-displacement curves predicted by the three meshes are very close to each other, and, in particular, the results for the medium and the fine mesh are almost identical.The results for the coarse mesh are slightly different due to the already discussed required mesh size for a sufficient resolution of the gradient of the nonlocal field.Figure 17 depicts the damaged zones in the rock mass at the level of 12% of the initial internal pressure for the three meshes employing full numerical integration.The highly damaged zones correspond to large shear deformations, which were shown previously in Figure 15.A similar shape of failure zones was obtained by Addis et al. [57] in laboratory tests on a bore hole in weak sandstone.Hence, the potential of the gradient-enhanced approach to predict the onset of failure of a structure was demonstrated in spite of the rather complex failure mode.

Conclusions
The present contribution addressed the prediction of complex failure modes in numerical simulations by means of a new implicit gradient-enhanced damage-plasticity model for intact rock and rock mass.It was derived from the damage-plasticity rock model presented by Unteregger et al. [24,25] and extended by adopting the implicit gradient-enhancement proposed in [17].For the assessment of the model, a comprehensive numerical study was presented.In particular, numerical simulations of wedge splitting tests for evaluating the representation of mode I failure, simulations of triaxial compression and extension tests for evaluating the representation of shear failure under confined stress states, and finally simulations of deep tunneling for examining the prediction of failure mechanisms of a complex structure were performed.From the obtained results of the numerical study, the following conclusions can be drawn:

•
In numerical simulations of wedge splitting tests, the formation of a crack propagating from the notch is modeled by the implicit gradient-enhancement in a smeared manner over a width related to the assumed length scale parameter.The capability of the gradient-enhanced damage-plasticity rock model of representing the experimentally observed material behavior was realistically demonstrated.

•
Finite element analyses with both structured and unstructured meshes confirmed the regularizing effect of the implicit gradient-enhancement in mode I failure and thus revealed mesh-insensitive results.In particular, it was demonstrated that the crack direction is not biased by the orientation of the finite element mesh.

•
Regarding the simulations of triaxial compression tests, a proper representation of shear failure under confined stress states, attaining the compressive meridian of the yield surface, was shown.After attaining peak strength upon initiation of damage, localization into a distinct shear band was observed.Furthermore, reasonable agreement with experimental results, and in particular the experimentally observed increasingly ductile material behavior with increasing confining pressure was obtained.

•
The simulations of triaxial extension tests demonstrated that failure under confined stress states attaining the tensile meridian of the yield surface is represented reasonably well, i.e., the localization of damage into a distinct shear band is predicted.In contrast to the triaxial compression tests, a more brittle structural response was observed.

•
For both the triaxial compression and extension tests, a mesh study confirmed mesh-insensitive results for shear failure under confined stress states.

•
In the numerical simulations of deep tunnel excavation, the rock mass was subjected to softening material behavior due to the excavation procedure, which leads to localization of strains into multiple distinct shear bands.In spite of the complex structural failure mechanism involving the formation of multiple shear bands, mesh insensitive load-displacement curves were obtained.By analyzing the convergence of the damage patterns upon mesh refinement, only slight differences were recognized.

•
For the adopted formulation of the implicit gradient-enhanced rock model, the length scale l is assumed as a constant parameter.Thus, a rather broad zone of completely damaged material is predicted by the model.Possible remedies were proposed in the literature, for instance, based on variable length scale parameters [58,59] or the so-called concept of decreasing interactions [45], which is motivated by a physically based micromechanical homogenization theory [60], cf.[46] for a discussion of these approaches.For a more realistic representation of sharp cracks within the present gradient-enhanced rock model, further investigations of these concepts are recommended.

•
Regarding the identification of the parameters controlling the gradient-enhanced softening part of the model, i.e., softening modulus ε f , length scale parameter l, and weighting parameter m, an ad hoc approach was employed: While the length scale parameter l was treated as a model parameter conforming to the size of the employed finite element mesh to ensure a sufficient resolution of the nonlocal gradient, ε f was identified by a best fit with experimental results for a chosen l and m.
A more systematic approach for identifying these parameters based on experimental results is an open issue.In particular, a scheme to compute ε f based on a prescribed value of l and a typical material parameter, such as the specific mode I fracture energy, is desirable.
Summarizing, it can be concluded that the presented damage-plasticity model is capable of representing a wide range of failure mechanisms in numerical simulations in a realistic and objective manner.

Figure 2 .
Figure 2. Geometry and boundary conditions of the numerical model of the specimen for the wedge splitting test.

Figure 3 .
Figure3.Splitting force vs. crack mouth opening displacement (CMOD) for the wedge splitting test: experimental results (data taken from[40]) and numerical results.

Figure 7 .
Figure 7. Geometry and boundary conditions of the specimen for the triaxial compression test.

Figure 11 .
Figure 11.Triaxial compression tests: distribution of the vertical displacement u z in the symmetry plane for the four different levels of confining pressure.

Figure 12 .Figure 13 .
Figure 12.Computed load-displacement curves for triaxial extension tests for different levels of confining pressure p (0) .

Figure 14 .
Figure 14.2D initial boundary value problem of deep tunnel excavation: full model (left) and the detail center view (right).

Figure 15 .
Figure 15.Distribution of the magnitude of the displacement vector in the vicinity of the tunnel surface at the level of 12 % of the initial internal pressure with a displacement scale factor of 10 for the three finite element meshes with full numerical integration: coarse (left), medium (center), and fine (right).

Figure 16 .
Figure 16.Load-displacement curves: normalized internal pressure versus the mean displacement magnitude at the tunnel surface for the three finite element meshes with full and reduced numerical integration: total view (left), detailed view (right).

Figure 17 .
Figure 17.Distribution of the damage variable ω in the rock mass in the vicinity of the tunnel surface at the level of 12 % of the internal pressure for the three finite element meshes employing full numerical integration: coarse (left), medium (center), and fine (right).