Large strain correction for tunnel analyses considering hydromechanical coupling and ground anisotropy

The deformations resulting from tunnel analyses for heavily squeezing ground may be very large, necessitating numerical formulations that consider geometric nonlinearity. Alternatively, for a certain class of problems, routine small strain analyses can be performed, and their results can be corrected to account for large strains by means of a simple hyperbolic expression proposed a few years ago. The present paper shows that this correction equation is sufficiently accurate for practical purposes even in the case of anisotropic material behaviour and for hydromechanically coupled, steady state or transient analyses of tunnels. The accuracy of the equation prediction varies amongst these cases but is satisfactory overall for the purpose of preliminary calculations, thus broadening its value and usefulness as a preliminary design tool.


Large strain correction for tunnel analyses considering hydromechanical coupling and ground anisotropy Matteo Natale , Alexandros N. Nordas * , Seraina Kopp & Georgios Anagnostou
The deformations resulting from tunnel analyses for heavily squeezing ground may be very large, necessitating numerical formulations that consider geometric nonlinearity.Alternatively, for a certain class of problems, routine small strain analyses can be performed, and their results can be corrected to account for large strains by means of a simple hyperbolic expression proposed a few years ago.The present paper shows that this correction equation is sufficiently accurate for practical purposes even in the case of anisotropic material behaviour and for hydromechanically coupled, steady state or transient analyses of tunnels.The accuracy of the equation prediction varies amongst these cases but is satisfactory overall for the purpose of preliminary calculations, thus broadening its value and usefulness as a preliminary design tool.

List of symbols E o
Young's modulus normal to the anisotropy plane E p Young's modulus parallel to the anisotropy plane G op Shear modulus on planes orthogonal to the anisotropy plane H w Depth of the tunnel underneath the water table R Outer radius of the region affected by tunnel drainage U a,ss Normalised small-strain radial displacement of the tunnel boundary U a,ls Normalised  Dilatancy angle in the anisotropy plane In the analysis of tunnelling through heavily squeezing rock, where deformations can be very large, the application of small-strain theory, which evaluates stiffness and equilibrium in the undeformed initial configuration, can significantly overestimate tunnel convergences and lead even to nonsensical predictions.Such problems must be solved instead in the framework of the large-strain theory which considers geometric nonlinearity by evaluating stiffness and equilibrium in the deformed configuration.
Vrakas and Anagnostou 1 showed that the differences between small-and large-strain analyses become considerable when convergences exceed 20% of the tunnel radius, demonstrating the limitations of small-strain analyses on prominent examples of tunnelling under heavily squeezing conditions: In a critical 6 km stretch of the Gotthard base tunnel in Switzerland, which crosses low stiffness and strength kakiritic rocks at depths up to 900 m, small-strain analysis was shown to overestimate the over-excavation required to accommodate the designated profile clearance by more than 50%.In the Yacambú-Quibor tunnel in Venezuela, where complete closure of the cross-section (convergence > 80%) occurred in a 23.3 km stretch crossing weak graphitic phyllites at depths up to 1270 m, a back-calculation based on small-strain analysis was shown to result in a fourfold overestimate of the ground cohesion.In a critical zone of the planned subaqueous Gibraltar Strait tunnel, which crosses breccia over 4 km at a depth of 200 m below the seabed, preliminary small-strain calculations of the undrained ground response were shown to yield physically unrealistic predictions (convergences > 100% vs. 48% from large-strain analyses).
While the above examples underscore the significance and necessity of adopting a large-strain formulation in the analysis of severe squeezing conditions, geometrically nonlinear finite element (FE) analyses can be excessively demanding in respect of computing time and resources.In view of this, Vrakas and Anagnostou 1 proposed the following simple, accurate and theoretically well-founded hyperbolic equation for obtaining large-strain solutions from small-strain analyses: where U a,ss = u a,ss /a and U a,ls = u a,ls /a, u a,ss and u a,ls are the displacements from small-and large-strain analyses, respectively, and a denotes the tunnel radius in the initial undeformed configuration.Equation (1) constitutes a valuable tool for preliminary design, since it entails performance only of routine small-strain analyses with a subsequent "self-correction" of the results.Although it has been derived considering the plane-strain and rotationally symmetric tunnel problem, it has been shown to be reasonably accurate also for other, general 2D problems (cylindrical and spherical openings, dilatant and hardening material behaviours, non-cylindrical tunnels and non-hydrostatic in-situ stress fields, undrained response of saturated low-permeability ground), as well as for 3D problems of an advancing tunnel heading.
Building upon the work of Vrakas and Anagnostou 1 , the present paper investigates the applicability of Eq. (1) (hereafter referred also as "correction equation") to anisotropic problems (Section Anisotropic stress analyses), and to hydromechanically coupled, steady state or transient problems with isotropic (Section Coupled isotropic analyses) or anisotropic (Section Coupled anisotropic analyses) material, which have not been considered previously.Throughout the paper, either isotropic or transversely isotropic rocks with a vertical axis of symmetry and horizontal planes of anisotropy (foliation, stratification) are considered.The constitutive behaviour of the rock is assumed as in Vrakas and Anagnostou 1 , i.e. linear-elastic and perfectly plastic, following the Mohr-Coulomb yield criterion and a non-associated flow rule, which is the most commonly used model in rock and tunnel engineering practice.

Problem layout and computational assumptions
Consideration is given herein to the plane-strain problem of a deep, cylindrical tunnel of radius a, crossing homogeneous, transversally isotropic rock and subjected to a homogeneous and hydrostatic in-situ stress field (Fig. 1).The computational domain considers one quarter of the tunnel cross-section, on account of symmetry about the horizontal and vertical planes, and a far field boundary (where the radial in-situ stress σ 0 prevailing at the elevation of the tunnel axis is applied) at a sufficiently large distance from the tunnel (b = 100 a; Fig. 1b).The ground response to tunnelling is simulated via an unloading of the tunnel boundary from the initial in-situ stress σ 0 to zero, which provides the maximum convergence.The numerical computations have been performed in Abaqus® 2 , employing a structured FE mesh of 1′120 4-noded, linear, quadrilateral, plane strain elements.
A linear elastic and perfectly plastic constitutive model with a Mohr-Coulomb yield criterion and a nonassociated flow rule is adopted, which considers transversal isotropy and reduced strength for shearing along an anisotropy plane compared to shearing through the rock matrix.The model behaviour is defined by 11 constants 3 : two Young's moduli E p and E o ≤ E p parallel and orthogonal to the anisotropy plane, respectively; two Poisson's ratios of the anisotropy plane ν pp , ν op for loading parallel and orthogonal to it, respectively; the cohesion c, angle of internal friction φ and angle of dilation ψ for shearing through the matrix; the three corresponding plasticity constants (c b , φ b , ψ b ) for shearing along an anisotropy plane; and the shear modulus G op on planes orthogonal to www.nature.com/scientificreports/ the anisotropy plane.The latter is an independent material constant of the transversely isotropic material and is taken equal to 0.5 E o /(1 + ν op ), as suggested by Wittke 4 .

Accuracy of the correction equation
Dimensional analysis makes it possible to express the normalised displacement of any point of the tunnel boundary in the following form: where f c is the uniaxial compressive strength disregarding strength anisotropy (f c = 2 c cosφ /(1-sinφ)), and θ denotes the angular position of a point on the tunnel boundary.The parameter sets considered are given in Table 1 and have been selected such that the normalised convergences from small-strain analyses (U a,ss ) range between 15 and 100%.The sets with c b /c = ∞ disregard strength anisotropy and are intended to consider the isolated effect of stiffness anisotropy.Anisotropy results in general in a non-uniform deformation of the tunnel boundary.For the assumed horizontal orientation of the anisotropy planes, stiffness anisotropy results in higher convergences at the crown (θ = 90°), where unloading takes place in the "softer" direction orthogonal to the planes, compared to the sidewall (θ = 0°), where unloading takes place in the "stiffer" direction parallel to the planes (see inset in Fig. 2).Strength anisotropy results in a local increase of convergences over a narrow region close to the crown, where the lower strength of shearing along the anisotropy planes prevails; however, this effect is in general small in comparison with the absolute magnitude of convergence considered.The crown and sidewall can thus be considered as reference points for evaluating the results, since their convergences are indicative of the maximum and minimum over the boundary, and their average value is also approximately equal to the average convergence over the entire boundary.
(2) www.nature.com/scientificreports/ Figure 2 shows the numerically determined ratio U a,ls /U a,ss as a function of the predictions U a,ss of small-strain analyses at the tunnel sidewall (triangular markers) and crown (square markers), along with the approximation of this relationship based on Eq. (1) (dashed line).One can readily verify that the correction equation approximates the numerical predictions with very high accuracy in the absence of strength anisotropy (red markers).In cases considering strength anisotropy (black markers), the correction equation is less accurate overall, mostly slightly underestimating the convergences at the crown and overestimating those at the sidewall; however, it still provides reasonably accurate prediction of the average convergence throughout the considered range, even for excessive convergences U a,ss close to 100%.

Problem layout and computational assumptions
Consideration is given in the sequel to a tunnel deep under the water table (H w = 100 a; Fig. 1a).The seepage flow field is approximated as rotationally symmetric, assuming that the pore pressure is equal to the in-situ hydrostatic pressure p 0 prevailing at the elevation of the tunnel axis at a radius R = H w (Fig. 1b); this simplification has been shown to be sufficiently accurate 5,6 .Identical computational specifications to those outlined in Section Anisotropic stress analyses are otherwise adopted, whilst also considering the linear interpolation of the pore pressure within each FE.
The rock is considered as a porous, fully saturated medium obeying Terzaghi's principle of effective stresses and isotropic Darcy's law.Two borderline cases are investigated for the permeability: a high permeability (k = 10 -6 m/s) where drained conditions prevail continuously; and a low permeability (k = 10 -13 m/s) where the conditions are undrained during tunnel excavation, and become drained after a sufficiently long transient consolidation period.The parameters adopted in the computations are given in Table 2.
The high-permeability case assumes that the hydraulic head field readily achieves steady state prior to excavation, as in the case of perfect advance drainage 7 .This is achieved by initially prescribing atmospheric pore pressure over the fixed tunnel boundary (p a = 0, where p a denotes the pore pressure at the tunnel boundary).Subsequently, the boundary is completely unloaded while maintaining the atmospheric pressure boundary condition.
The low-permeability case initially considers a complete unloading of the tunnel boundary assuming a noflow hydraulic boundary condition (q a = 0, where q a denotes the flux at the tunnel boundary), which provides the instantaneous undrained convergence.Subsequently, a transient consolidation analysis is performed until steady state is achieved, considering a mixed hydraulic boundary condition over the tunnel boundary.This ensures water seepage from the ground into the tunnel, but not vice versa (q a = 0 if ∂p/∂r < 0, p a = 0 if ∂p/∂r > 0); 8 , Table 1.Parameter sets of the anisotropic stress analyses (parameters common to all sets:  www.nature.com/scientificreports/if an atmospheric pressure condition (p a = 0) was specified instead, watering of the ground from the tunnel could occur in the early stages of consolidation, when negative pore pressures may develop.The effect of the adopted hydraulic boundary condition is illustrated qualitatively for two time-instances (t 1 , t 2 ) in Fig. 3.

Accuracy of the correction equation
Dimensional analysis enables expression of the normalised convergence of the tunnel boundary in the following form: The parameters adopted in the computations are given in Table 2.
Figure 4 shows the numerically determined ratio U a,ls /U a,ss as a function of the predictions U a,ss of small-strain analyses (round markers) for the high-permeability case, along with the approximation of this relationship based on Eq. (1) (dashed line).The correction equation achieves a very good agreement with the numerical predictions, becoming only slightly conservative towards the range of excessively large U a,ss .
The results for the case of low permeability are shown in Fig. 5, where it is readily seen that the correction equation achieves a virtually perfect match with the numerical predictions.This is expected for the instantaneous convergences, since these lie in the range where geometric nonlinearity effects are limited (U a,ss < 15%), but holds also for the steady-state convergences, even in the range of excessively large, nonsensical values (U a,ss > 100%). (3) Anisotropic stress analyses: ratio of the large-to the small-strain displacements versus the smallstrain displacements (parameters: s.Table 1; labels besides the markers: data-set numbers after Table 1).
Figure 5. Isotropic coupled analyses for low-permeability ground: ratio of the large-to the small-strain displacements versus the small-strain displacements (parameters: s.

Coupled anisotropic analyses
The applicability of the correction equation to coupled anisotropic problems is examined considering the computational model discussed in Section Coupled isotropic analyses, with the additional specification of an anisotropic material behaviour analogously to Section Anisotropic stress analyses.The parameters adopted in the computations are given in Tables 3 and 4 for the high-and low-permeability cases, respectively.In addition to the strength and stiffness anisotropy, the case of a pronouncedly anisotropic seepage flow was considered in certain cases, assuming that the permeability parallel to the bedding or schistosity plane is 10 times higher than perpendicularly to it.Figure 6a assumes isotropic permeability and shows the numerically determined ratio U a,ls /U a,ss as a function of the predictions U a,ss of small-strain analyses, along with the approximation of this relationship based on Eq. (1) (dashed line) for the high-permeability case.The correction equation is very accurate in the absence of strength anisotropy (red markers) and reasonably accurate otherwise (black markers), slightly underestimating the displacements at the crown and overestimating those at the sidewall in some cases; overall, the accuracy achieved is very good in terms of average convergences.This also applies in the case of anisotropic permeability (Fig. 6b).
Figure 7a shows the results for the low-permeability case, considering isotropic permeability.The instantaneous (undrained) convergences U a,ss are anyway small (< 15%) and thus accurately predicted even by the infinitesimal-strain analysis.As for the steady-state convergences, the accuracy of the correction equation is high in the absence of strength anisotropy (red markers) but inferior otherwise (black markers), particularly for the crown, where the correction equation may underestimate displacement by up to 40% (data set 21).There are also some outliers (tunnel wall; data sets 15, 16 and 22) for which even the correction equation overestimates displacement.Notwithstanding the inaccuracies observed in the case of strength anisotropy, the correction equation provides a reasonable indication of the anticipated magnitude of convergences; therefore, it is still valuable for preliminary calculation purposes, particularly in cases where small-strain analyses provide unusable results, i.e. predicting excessive and nonsensical convergences (U a,ss > 100%).The same conclusions can be drawn when considering permeability anisotropy (Fig. 7b).

Concluding remarks
By utilising parametric plane strain numerical simulations, the present paper demonstrated that the equation proposed by Vrakas and Anagnostou 1 for obtaining large-strain solutions from small-strain analyses can be applied also to anisotropic and hydromechanically coupled, steady state or transient analyses.The accuracy of the equation has been shown to be: (i) sufficient for estimating average convergences in anisotropic ground (Fig. 2); (ii) very high for estimating instantaneous undrained convergences in low-permeability isotropic (Fig. 5) or Figure 6.Anisotropic coupled analyses for high-permeability ground: ratio of the large-to the small-strain displacements versus the small-strain displacements assuming, (a) isotropic permeability or, (b) isotropic and anisotropic permeability (parameters: s.Table 3; labels besides the markers: data-set numbers after Table 3).anisotropic (Fig. 7) ground; (iii) very high for estimating steady state convergences in high-permeability (Fig. 4) and low-permeability (Fig. 5) isotropic ground; (iv) reasonable for estimating steady state convergences in highpermeability (Fig. 6) and low-permeability (Fig. 7) anisotropic ground.Complementing the work of Vrakas and Anagnostou 1 , the present findings broaden even further the already wide applicability of the correction equation as a preliminary design tool.

Figure 1 .
Figure 1.(a) Problem layout and, (b), model (data in brackets: applies to the coupled analyses only).

Figure 4 .
Figure 4. Isotropic coupled analyses for high-permeability ground: ratio of the large-to the small-strain displacements versus the small-strain displacements (parameters: s.Table2).
large-strain radial displacement of the tunnel boundary b Angle of internal friction in the anisotropy plane ν op Poisson ratio expressing the effect of the normal stresses orthogonal to the anisotropy plane to the normal strains parallel to the anisotropy plane OPEN ETH Zurich, Stefano-Franscini-Platz 5, 8093 Zurich, Switzerland.* email: anordas@ethz.chwww.nature.com/scientificreports/

Table 3 .
Parameter sets of the anisotropic coupled analyses for high-permeability ground (parameters common to all sets:

σ o f c /σ o E p /E o c b /c φ b /φ p 0 /σ 0
100). a For this data-set the case of permeability anisotropy with k h /k v = 10 has been additionally considered.

Table 4 .
Parameter sets of the anisotropic coupled analyses for low-permeability ground (parameters common to all sets: ν pp = ν op = 0.2, φ = 25°, ψ = ψ b = 5°, b/a = R/a = 100).a For this data-set the case of permeability anisotropy with k h /k v = 10 has been additionally considered.