Three-dimensional numerical simulation of rock deformation in bolt-supported tunnels: A homogenization approach

Three-dimensional numerical simulation of rock deformation in bolt-supported tunnels: A homogenization approach

Tunnelling and Underground Space Technology 31 (2012) 68–79 Contents lists available at SciVerse ScienceDirect Tunnelling and Underground Space Tech...

1MB Sizes 0 Downloads 116 Views

Tunnelling and Underground Space Technology 31 (2012) 68–79

Contents lists available at SciVerse ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Three-dimensional numerical simulation of rock deformation in bolt-supported tunnels: A homogenization approach Samir Maghous ⇑, Denise Bernaud, Eduardo Couto Department of Civil Engineering, Federal University of Rio Grande do Sul, Porto Alegre-RS, Brazil

a r t i c l e

i n f o

Article history: Received 7 October 2011 Received in revised form 27 February 2012 Accepted 2 April 2012 Available online 28 April 2012 Keywords: Tunnel Bolted rock Bolt/rock interface Finite element method Homogenization Embedded model

a b s t r a c t Despite fully grouted bolts are nowadays widely used for tunnel support, bolting design is still based in many cases on empirical or semi-empirical considerations. This paper describes a three-dimensional theoretical and numerical model for the behavior of tunnels reinforced by bolts. From a theoretical viewpoint, the elastoplastic constitutive equations for the reinforced rock are derived in the framework of homogenization method. Emphasis is given to the formulation of the homogenized strength criterion with account for the bolt/rock interface properties. The anisotropic constitutive equations are then implemented in a 3D finite element computer code in which the processes of excavation, installation of bolts and lining placement are simulated by means of the ‘‘activation/deactivation’’ technique. Due to the multi-potential nature of the plastic flow, a specific iterative algorithm for plastic integration is proposed. The finite element model is applied to the analysis of the Kielder experimental tunnel for which in situ measured data are available. The accuracy of the numerical predictions based on homogenization method is also assessed by comparison with the results derived from implementation of the so-called embedded model. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The idea of installing metallic bolts within the ground for tunnel support has originally emerged in 1913 with the request of a technical patent to Germany authorities (Kovari, 2003). However, their effective use began only in the 1940s as roof support in the context of the American mining industry. From 1970s up today rock bolting has experienced a rapid increase and nowadays it is widely applied in both civil and mine tunnel projects. The success of such reinforcing system can be explained by its ease to carry out since it does not require heavy equipment for installation and by its relative flexibility in the sense that bolt density can be punctually adapted to account for local geotechnical properties. Furthermore, the use of bolts is an essential component of the New Austrian Tunneling Method (NATM) or the so-called ‘convergence-confinement’ approach for tunnel construction (Panet, 1995) and, for this reason, application of these approaches has accompanied the expansion of this reinforcement technique for tunnels excavated in soft rocks. Schematically rock bolts fall within the general category of rock reinforcement systems that act as structural inclusion installed within the ground with the objective to improve the behavior of the ground (Windsor and Thompson, 1993). Mechanically speaking, the mechanisms governing the interaction prevailing between ⇑ Corresponding author. Tel.: +55 51 33 08 35 88; fax: +55 51 33 08 39 99. E-mail address: [email protected] (S. Maghous). 0886-7798/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tust.2012.04.008

the rock bolts and the surrounding rock mass are complex and still not well understood. Hence, the modeling of such interaction is still a work in progress, and a comprehensive mechanical description is needed for reliable predictions regarding the performance of the rock bolting. The reader may refer for instance to the recent paper by Bobet and Einstein (2010) for a review on the technique of rock bolting. The geocomposite defined by association of rock bolts and the rock mass is strongly heterogeneous. Due to the large number of rock bolts involved in such a reinforcement technique, the implementation of direct numerical finite element analyses based on the explicit description by anchor element embedded within the rock mass, would require the use of sophisticated computational capabilities. Indeed, an accurate 3D geometrical discretization of the reinforced zones should involve rock mass finite elements whose typical size is close to the bolt diameter. This results in tedious mesh generation procedures with a high number of degrees of freedom. These difficulties have given rise to alternatives approaches to handle this kind of problem. On one hand, the homogenization approach is a theoretical-based alternative which offers an appropriate framework in the field of rock bolting. On the other hand, the so-called ‘‘embedded model’’ appears as a quite efficient numerical direct method to simulate bolt-supported tunnels in the case of perfect bonding assumed for the interface bolts/ surrounding rock. In terms of numerical implementation of both methods, the rock bolts are not treated as individual structural

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

elements to be geometrically discretized. This results in a simplified mesh generation that is similar to that used for a homogeneous rock. Basically, the reinforced medium is replaced in the homogenization method by an equivalent homogenized medium. The properties of the latter are derived from those of the individual constituents (i.e., rock and bolts) using the theoretical tools of homogenization for periodic media. This approach has been successively applied to model at a macroscopic scale the behavior of soil structures reinforced by stiff linear inclusions (see for instance de Buhan et al., 1989; de Buhan and Salençon, 1990). The pioneer works applying the homogenization approach may be attributed to Wullschläger and Natau (1983, 1987), Greuell (1993), and Greuell et al. (1994), who analyzed the simplified situation of a rock bolted cavity under axisymmetric plane strain condition. The formulation of an anisotropic elastoplastic behavior for the bolted rock regarded as a homogenized medium as well as the corresponding finite element implementation have been originally presented in Bernaud et al. (1995). Extension of these works to the situation of frictional rock material has been proposed within the context of a 2D setting by Bernaud et al. (2009). As far as design of bolt reinforced tunnels is concerned, several works have been dedicated, in the last decades, to develop analytical and/or numerical models aimed at providing useful methods for practical use in underground engineering. These works cover both aspects related to deformation and stability analyses in bolt-supported tunnels. Among the contributions related to the effect of face or radial bolting, the reader may refer for instance to Peila (1994), Peila and Oreste (1995), Oreste and Peila (1996), Peila et al. (1996) or Peila et al. (2000), to cite a few. Interestingly enough, Lunardi (2000) has suggested that understanding and controlling the behavior of the core ahead of the advancing tunnel face is crucial for successful tunneling in squeezing ground conditions. Based on this assumption, this author proposed a comprehensive approach for design and construction of tunnels that relies upon on the analysis of controlled deformation in rocks and soils. Conceived as an extension of the homogenization approach, the multiphase model has been applied to the prediction of deformation in a bolt-supported tunnel by Sudret and de Buhan (2001) and de Buhan et al. (2008). With respect to homogenization approach, the multiphase approach enables accounting for a possible slippage at the interface rock/bolt. It can thus be regarded as an improved homogenization procedure. However, first applications of the multiphase model in the context of reinforced soils have led to similar results than those obtained by homogenization (Sudret et al., 1998). The new aspects of the current work, which is an extension of that developed in Bernaud et al. (2009), are fundamentally:  the incorporation of a theoretical model for the bolt/rock interface within the homogenized elastoplastic constitutive equations;  the finite element implementation of the homogenized elastoplastic model within a three-dimensional setting;  the finite element implementation of the embedded model as a new alternative to homogenization method for the simulation of deformation in bolt reinforced tunnels. 2. Anisotropic elastoplastic constitutive equations for the bolted medium 2.1. Principle of the homogenization approach From a theoretical point of view, the homogenization of periodic media method requires the following assumptions to be fulfilled:

69

 the inclusions are arranged in a regular pattern;  the density of the bolt array is high enough, i.e., the spacing between the inclusions can be considered small when compared to a characteristic dimension of the problem. From a practical viewpoint, the above assumptions are met by the majority of bolt-supported tunnels: radial or frontal bolts are linear reinforcements which are regularly disposed within the ground and the bolt spacing and diameter are much smaller than the radius of the tunnel. Taking advantage of these specificities and that the bolt material is much stiffer than the rock material, it can be proved (Greuell, 1993; Bernaud et al., 1995; Bernaud et al., 2009) that the composite material obtained by association of rock and bolts behaves as a transversely isotropic medium around the direction of the direction of bolts. The tunnel considered in the following formulation is horizontal cylindrical opening parallel to direction ez and having a horseshoetype cross section as sketched in Fig. 1. Three kinds of rock bolts disposed perpendicular to the tunnel wall can be used separately or in association:  radial bolts, placed in the roof following the radial direction er;  side bolts, placed along the vertical wall sides of the tunnel parallel to direction ex;  face bolts, placed at the tunnel face parallel to direction ez. The reinforced zone defined above will be referred to in the sequel by (RZ), (SZ) and (FZ) respectively. The representative elementary volume is prismatic in (RZ) while it is parallelepipedic in regions (SZ) and (FZ). The radial bolts of length lr are regularly distributed following a horizontal spacing pf and angular spacingbf. The wall side bolts (resp. face bolts) have a length ls (resp. lf) and are distributed following a rectangular grid ps  bs (resp. ps  bs). The high heterogeneity of the composite material formed by rock and bolts, combined with the large numbers of inclusions generally involved, makes the homogenization method more attractive to deal with the problem of tunnel reinforced by bolts. This is mainly attributed to the fact that for numerical analyses, there is no need of a fine dicretization for capturing the strong heterogeneity of bolted zones. Clearly enough, the reinforcement patterns can be considered as dense in the zones (RZ), (SZ), and (FZ) if all the quantities br/p, pr/R, pf/R, ps/R, bf/R, bs/R are much small than unity. Length R refers to a characteristic dimension of the tunnel cross-section. The formulation of the constitutive equations for the bolted rock regarded as a homogenized medium requires defining the relevant reinforcement parameters. Firstly, the volume proportion of reinforcement is expressed as:

8 > gr ðrÞ ¼ prSbrr r > > < S g ¼ gs ¼ ps bs s > > > : g ¼ Sf f p b f f

for

ðZRÞ

for

ðZSÞ

for

ðZFÞ

ð1Þ

where Si, i 2 {r, s, f}, is the cross-section area of the bolts. This nondimensional parameter is constant in (ZS) and (ZF), while it is a decreasing function of the radial distance r in (ZR). Secondly, the density of reinforcement d is the ratio on the number of bolts per unit area on the tunnel wall or tunnel face:

8 dr ¼ p b1r R0 > > r < d ¼ df ¼ ps1bs > > :d ¼ 1 f p b f f

for

ðZRÞ

for

ðZSÞ

for

ðZFÞ

ð2Þ

where parameter R0 stands for the radius of the circular part of the tunnel cross-section.

70

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

A

pr B

lr

pf z

lf B'

A'

bf

ls

βr

pf

ps

R

R

Section AA'

Section BB'

Fig. 1. Geometry of the bolted tunnel and representative elementary volumes.

In the subsequent analysis, the rock material is modeled considered as a homogenous material. This assumption implies that, at the scale of modeling, the bolts are the only heterogeneity considered for the bolted rock. It is implicitly intended that this homogeneity results from a previous homogenization process that account for the possible micro-heterogeneities present within the rock. Hence, what will be referred to as the rock matrix is the equivalent medium whose behavior is obtained by homogenization. Several works have been devoted to this topic and the interested reader may refer for instance to de Buhan and Maghous (1997), Maghous et al. (2002), Deudé et al. (2002), Bernaud et al. (2002), Maghous et al. (2008) or Maghous et al. (in press). 2.2. Homogenized elastic moduli Within the framework of infinitesimal strains, the homogenized elastic relationship for the reinforced rock mass relates the macroscopic stress R and strain 2

R ¼ Ahom : e

ð3Þ

where Ahom is the fourth-order tensor of homogenized moduli. Variational principles in elasticity applied to the elastic concentration problem stated on the representative elementary volume of the bolted rock lead to the following result (Greuell, 1993):

Ahom ¼ Am þ A ei  ei  ei  ei with K ¼ gEb

ð4Þ

where Am stands for the elastic tensor of the rock and Eb is the Young’s modulus of the bolt disposed following direction ei. In the above equation, the second term in the right-hand side accounts for the contribution of bolts to the overall elastic stiffness. Similarly to the volume fraction of reinforcement, it remains constant in (ZS) and (ZF), while it decreases with distance r in (ZR). Eq. (4) emphasizes that the bolted rock behaves elastically as a transversally isotropic medium about direction ei of reinforcement. If elastic isotropy is adopted for the rock material, the longitudinal and transversal Young’s modulus are given by

71

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

Ehom ¼ Em þ K; l

Ehom ¼ Em t

Em þ K   Em þ 1  m2m K

Σhk

ð5Þ

Ghom

whereas the overall shear moduli are not affected by the presence of bolts, that is

lhom ¼ lhom ¼ lm l t

ð6Þ

Gm

Ghom

Σ ii

2.3. Macroscopic yield function Let n denotes the unit normal vector to the bolt/rock interface. The strength condition at the interface bolt/rock consists classically in a limitation of the values taken by the stress vector r  n acting upon this interface. In the context of the present analysis, the bolts are modeled as linear inclusions and, hence, n \ ei. Accordingly, the failure condition for interface shall be expressed in the following form

fint ðr  nÞ 6 0

ð7Þ

Starting from the strength properties of the individual constituents of bolted rock, the theoretical framework for the determination of the macroscopic failure criterionfhom (R) 6 0 is that of limit analysis homogenization of periodic media. Extending the results derived in de Buhan and Taliercio (1991) and Bernaud et al. (2009), the yield function fhom may be expressed as:

f

hom

ðRÞ ¼ maxð~f ðRÞ; f int ðR  nÞÞ 6 0

ð8Þ

The failure function ~f is defined by

X X  ~f ¼ minfm  rei  ei r2I

ð9Þ

where fm(r) 6 0 is the yield criterion of the rock material, i 2 {r, s, f} depending on the considered reinforced region (ZR), (ZS) or (ZF). The interval I is defined by

I ¼ ½kgr0 ; gr0 

ð10Þ

r0 represents the strength under uniaxial tension of the bolts and k is a non-dimensional number ranging between 0 and 1 that accounts for a reduction in the compressive strength of the bolts due to structural buckling. Clearly enough, the values of r0 and k may depend on the reinforced zone under consideration. The quantity gr0 represents the strength of the bolt per unit transverse area at the current point. The proof of (8) is based on the static approach of limit analysis, in which a statically admissible distribution of stress that is piecewise homogeneous on the representative elementary volume is considered. In view of further developments and sake of notation simplicity, the strength interface will be conveniently expressed through a function of the macroscopic stress tensor instead of the stress vector as follows:

~f ðRÞ ¼ f ðTÞ with T ¼ R  n int int

ð11Þ

T = R  n is the stress vector acting upon the bolt/rock interface. The strength macroscopic strength domain Ghom defined by the yield condition fhom(R) 6 0 is the intersection of two domains:

e\G e int Ghom ¼ G

ð12Þ

e (resp. G e int ) corresponds to yield criterion ~f ðRÞ 6 0 where convex G (resp. ~f int ðRÞ 6 0Þ. A schematic representation of the macroscopic strength domain e Ghom is given in Fig. 2. It may be drawn as follows: First, domain G is the convex envelope of two domains obtained by translating the rock strength domain Gm by algebraic distances  kgr0 and gr0

−k η σ 0

ησ0

Fig. 2. Geometrical interpretation of the macroscopic strength domain.

along Rii-axis. The convex Ghom is then obtained as the intersection e and G e int . of domains G The particular case of perfect bonding at the interface bolt/rock corresponds to Gint ¼ R6 , and therefore

e i:e:; f hom ðRÞ ¼ ~f ðRÞ Ghom ¼ G;

ð13Þ

2.4. Elastoplastic constitutive law The formulation of the elastoplastic behavior of the bolted rock regarded as a homogenized medium relies upon the heuristic idea that the macroscopic strength domain Ghom can be viewed as the elasticity domain of the material with the elastic properties defined in Section 2.2. This procedure simply amounts to neglecting the hardening induced at the macroscopic level by the homogenization process (Suquet, 1985). Accordingly, the rate-type state equation of the homogenized medium reads: p

R_ ¼ Ahom : ðe_  e_ Þ

ð14Þ

p where e_ is the total strain rate and e_ is the corresponding plastic part. As mentioned above, the macroscopic strength domain is

adopted as the elasticity domain for the homogenized material. This means that the yield criterion fhom(R) 6 0 actually constitutes p

the plasticity criterion. The direction of e_ is prescribed by means of the plastic flow rule:

e_ p ¼ v_

@g hom ðRÞ @R

ð15Þ

where v_ P 0 denotes the plastic multiplier which complies with the Kuhn–Tucker condition v_ f hom ¼ 0. From a theoretical viewpoint, the question of how the macroscopic plastic potential ghom can be derived from the plastic properties of the rock material, the bolt material and those of the interface bolt/rock reveals complex and still remain an open issue. In the particular case where the assumption of associated plastic rule can be adopted for the constituents of the homogenized medium, it can be proved (Suquet, 1985) that normality holds also at the macroscopic scale, i.e., p

g hom ¼ f hom ) e_ ¼ v_

@f hom ðRÞ @R

ð16Þ

In a more general situation, the plastic flow rule falls within the category of multi-potential flow rule situation, which is associated with the fact that the boundary @Ghomof domain Ghom is a

72

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

notion of subdifferential of fhom is classically introduced to handle

f int = 0

hom

the singularity: 0 @f@ R 0 belongs to the cone delimited by normal ~

~

vectors @@rf ðRÞ and @@f rint ðRÞ. The case of a rock material obeying a Drucker–Prager yield criterion has been specifically treated in Bernaud et al. (2009) considering perfect bonding between rock and bolts. The corresponding expressions shall be only recalled hereafter. Introducing the norm kak = (a:a)1/2, the Drucker–Prager yield function may be expressed in the following form:

∂f hom ∂Σ

f hom = 0

G hom

fm ðrÞ ¼ a tr

pffiffiffi

r þ bksk  2C k

ð22Þ

where s ¼ devðrÞ ¼ r  rÞ1 is the deviatoric part of r, C is the cohesion and parameter k is defined from the friction angle um of the rock material by 1 ðtr 3

f =0 ∂f ∂Σ

hom



Fig. 3. Normal vector to the boundary of strength convex domain Ghom.

non-smooth yield surface. More precisely, three distinct situations can occur regarding the plastic flow rule at the macroscopic level. They are detailed in Appendix A. 2.5. Further characterization of the constitutive behavior The numerical implementation of the anisotropic elastoplastic behavior presented previously requires making explicit some mathematical aspects of the formulation. The practical applicability of the approach relies in particular on the ability to compute for a P given stress state R the value taken by the yield function ~f ð Þ. The

convexity

of

function

fm

implies

that

function

r ? fm(R  rei  ei) admits a single minimum x0 = x0(R) 2 R defined by equation

@fm ðR  x0 ei  ei Þ ¼ 0 @ rii

ð17Þ

~f ðRÞ ¼ fm ðR  r ei  ei Þ

ð18Þ

where r = r (R) is the projection onto the interval I of x0(R):

if x0 ðRÞ 6 gr0 if x0 ðRÞ 2 I

ð19Þ

0 1 C 3B a 3 x0 ðRÞ ¼ @Sii þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S  2 Sii devðei  ei Þ A 2 2 b  32 a2

ð24Þ

In the above equality, S = dev (R) denotes the deviatoric part of the macroscopic stress tensor andSii its component following the bolt direction ei, i 2 {r,s,f}. The derivative of ~f reduces to

S  s @~f @fm ðRÞ ¼ ðR  r ei  ei Þ ¼ a1 þ b @r @r kS  s k

ð25Þ

where s⁄ = dev (r⁄ ei  ei). Remark. From a practical viewpoint, the plastic properties of the interface are classically defined in terms of the stress vector T = R  n acting upon the interface: fint(T) 6 0 (yield function) and of the rate of displacement jump). In the context of the present modeling, the following relationships prove useful: s @ g~int @g ðRÞ ¼ int ðTÞ  n @R @T

ð26Þ

s

where  is the symmetric part of the tensorial product operator.

if x0 ðRÞ P gr0

Additionally, it is convenient to specify, for numerical implementahom tion purposes, the expression of the derivative @f@R ðRÞ of the yield function. In view of definitions (8) and (11)

f hom

The non-dimensional parameters a and b are pffiffifficomputed from the value of k as a = (k  1)/3 and b ¼ ð2k þ 1Þ= 6. The expression of x0 = x0(R) that allows defining r⁄ = r⁄(R) by (19), which in turn is necessary for evaluating ~f ðRÞ, is given by

s @~f int @fint ð RÞ ¼ ðTÞ n; @R @T



8 > < gr0 r ðRÞ ¼ projI x0 ðRÞ ¼ x0 ðRÞ > : gr0

ð23Þ

gint(T) 6 0 (plastic potential prescribing the direction of plastic part

Hence,



1 þ sin um 1  sinum

 

     R ¼ max ~f R ; ~f int R 60

ð20Þ

and extending the result established by Bernaud et al. (2009) for the adherent case, it follows that

8 e > < @@rf ðRÞ ¼ @f@ rm ðR  r ei  ei Þ if ~f ðRÞ > ~f int ðRÞ @f hom ðRÞ ¼   > @R : @~f int R if ef ðRÞ < ~f int ðRÞ @r

ð21Þ

The situation ~f ðRÞ ¼ ~f int ðRÞ corresponds to a singular case and should be handled carefully. Since the considered stress tensor R is actually located at a singular point of @Ghom (Fig. 3), the derivative of fhom is not defined from a mathematical viewpoint. The

3. Finite element implementation This section deals with the three-dimensional numerical implementation of the anisotropic constitutive law formulated in Section 2 for the homogenized bolted rock. To assess the accuracy of the homogenization-based model, the so-called ‘‘embedded method’’ that is classically used for the numerical modeling of reinforced concrete, has been also implemented and the corresponding results compared to those obtained by homogenization. 3.1. Homogenization model The anisotropic elastoplastic behavior for the bolted material regarded as homogenized medium has been implemented within a 3D finite element code specifically devised for the analysis of tunnel deformation during the different stages of construction. The simulation of excavation/advancing face/reinforcement placement

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

is carried out by means of the ‘‘activation/deactivation’’ technique originally developed by Hanafy and Emery (1980) for the simulation of advancing face and lining placement. The method has been later extended within a 2D setting by Bernaud et al. (1995, 2009) to account for the phase of bolt installation. The iterative algorithm for plastic integration (return mapping) requires being able to compute the projection of any stress state onto the convex yield domain Ghom. However, a specific treatment is needed herein due to the multi-potential nature of the plastic flow (see Eq. (A.3) of Appendix A), which is associated with the fact that @Ghom is a non-smooth yield surface. The difficulty arises when yielding occurs in bolt or/and rock (i.e., ~f ðRÞ ¼ 0Þ together with yielding at the interface (i.e., ~f int ðRÞ ¼ 0Þ. In other words, the stress sate lies on a singular point of @Ghom. As a matter of fact, most of the numerical procedures dealing with this aspect that are available in the literature are devoted to usual isotropic yield surfaces (e.g., Lade and Nelson, 1984; Crisfield, 1991; Simo and Hughes, 1998). The procedure proposed in the subsequent analysis is an extension of that developed by Maghous et al. (2008) for jointed rocks. Let R be a given stress state (obtained for instance from an elastic computational analysis) lying outside the yield domain, i.e., such that fhom (R) > 0. It follows from definition (23) that two cases are to be distinguished:  Case A: only one of the two ‘constituents’ (rock + bolts) or (interface material) defining the homogenized material has yielded. The situation when yielding occurs in the (rock + bolt) constituent, i.e., ~f ðRÞ > 0, ^f ðRÞ < 0, has been addressed in Bernaud

linked to the unnecessary matrix mesh refinement. The embedded method has originally been developed by Elwi and Hrudey (1989) for the analysis of reinforced concrete structures. As regards the numerical implementation, this method presents similar specificities than the homogenization method in the sense that the bolts are not descretized individually in the mesh generation. By nature, the effectiveness of its implementation is however limited when slippage at the interface bolt/ground should be considered. Assuming a perfect bonding between bolts and the surrounding rock mass, the principle of the embedded model is outlined in the following. The initial step of the embedded model is the geometrical discretization of the rock mass into finite elements (brick elements). This mesh generation can be achieved independently of the bolts disposition within the ground. Detection of the intersecting points of the embedded reinforcement bolts in the brick elements is then performed. Consequently, the bolts are described by means of straight segments, the position of each one being completely defined in the global coordinate system by the corresponding intersecting points. Each reinforcement bar segment is endowed with the same kinematics that the embedding rock brick element. The presence of the reinforcement bar segments do not therefore induce supplementary degrees of freedom than those associated with the displacement of the rock mass, because the inclusions move jointly with the brick elements in which they are embedded. The contribution of bolts to the general stiffness matrix is accounted for by means of an increasing term in the expression of the crossed element stiffness. More precisely, the element stiffness matrix may be expressed as:

int

et al. (2009), where a fully analytical solution has been derived for the projection of R onto the convex yield domainGhom. The alternative situation, defined by the yielding occurring at the interface constituent of the homogenized material, i.e.,~f ðRÞ < 0, ^f int ðRÞ > 0, is classically addressed in computational plasticity handbooks (Tresca or Mohr–Coulomb interface).  Case B: both the (rock + bolts) and the (interface material) constituents undergo yielding, i.e., ~f ðRÞ > 0, ^f ðRÞ > 0. The pro-

73

K ¼ Km þ

ns X Kib

ð28Þ

i¼1

where subscript m stands for the rock matrix, b for bolts and ns is the number of reinforcement bolt segments embedded in the considered brick element. For the rock matrix stiffness, one has:

Km ¼

Z

T

Bm Dm Bm dV

ð29Þ

V

int

jection of the stress state R, denoted by R⁄, onto the convex Ghom is defined by

8   > < R ¼ R  Ahom : v_ 1 @ g~ ðR Þ þ v_ 2 @ g~int ðR Þ @R @R > :

v_ 1 P 0; v_ 2 P 0

where Dm is the rock secant constitutive of matrix and Bmthe wellknown rock strain–displacement matrix. For the reinforcement, the element stiffness is

Kb ¼ Sb Eb

Z

T

Bb Bb dS

ð30Þ

S

f hom ðR Þ ¼ maxð~f ðR Þ; ~f int ðR ÞÞ ¼ 0 ð27Þ

A two stage iterative procedure is proposed for the determination of R⁄. The procedure is summarized in Appendix B. 3.2. Embedded model In order to assess the accuracy of the homogenization-based model, the so-called ‘‘embedded method’’ has been also implemented for the three-dimensional numerical simulation of rock bolting. In the particular case of reinforced concrete and more generally for the simulation of a matrix material reinforced by linear inclusions, there are traditionally two three-dimensional numerical simulation approaches: the discrete model and the embedded model. In the framework of discrete model, referred to as direct approach throughout the present paper, reinforcement inclusions are often modeled as bar elements located along the matrix element nodes. As underlined in the Introduction section, the drawbacks of the discrete (direct) approach lay in its high computational cost

Sb is the transversal area of bolts, Eb is the Young’s modulus and Bbis the strain–displacement matrix for the bolts, which relates rock brick element displacements to the axial bolt strains. In the objective to perform 3D numerical simulations for comparison with the homogenization approach, the embedded model described above has been implemented and incorporated in a finite element code devised for simulating deformation of bolted-supported tunnels. The numerical formulation is based on a standard elastoplastic algorithm. 4. Application: the case of Kielder tunnel To evaluate the ability of the homogenization model to predict the convergence of a bolt reinforced tunnel, the case of Kielder experimental tunnel is analyzed herein. This underground structure in Northeast England was driven at100mdepth in mudstone. In the 1970s, Freeman (1978) performed pioneering work in studying the performance of fully grouted bolts in the Kielder experimental tunnel. He monitored both the loading process of the bolts and the stress distribution along the bolts.

74

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79 Table 3 Mechanical parameters of rock material. Young modulus Em Poisson coefficient vm Friction angle um Dilatancy angle wm Cohesion C Specific weight c Coefficient of horizontal thrust k0

Fig. 4. Transversal section of the Kielder experimental tunnel and reinforcement pattern.

Table 1 Geometrical and mechanical and parameters of bolts. Bolt diameter Bolt length lr Angular distance br Horizontal spacing pr Young modulus Eb Yield strength r0

25  103 m 1.8 m 25° 0.9 m 210 GPa 650 MPa

Table 2 Geometrical and mechanical and parameters of lining. Lining thickness Young modulus Poisson coefficient

0.14 m 28 GPa 0.2

For simulation purposes, it is assumed that the initial stress field is geostatic and takes the form:

r0 ¼ r0V ey  ey þ r0H ð1  ey  ey Þ

ð31Þ

ey is the upward vertical direction. Denoting by c the specific weight of rock, the vertical and horizontal stresses are given by:

r0V ¼ ch; r0H ¼ k0 ch

ð32Þ

where h is the depth at the current point, counted from the upper 0 surface of the ground. Parameter k ¼ r0H =r0V stands for the horizontal thrust coefficient.

1200 MPa 0.3 26.4–36.4° 10° 0.24 MPa 25 kN/m3 0.75–1.0

Description of geometry and most of mechanical properties have been provided in several publications (see e.g., Freeman, 1978; Ward et al., 1983; Cai et al., 2004; Oreste, 2008). The tunnel has a circular section of 3.3 m and the bolt reinforcement scheme is displayed in Fig. 4. In addition to rock bolts, a sprayed concrete lining has also been used. The model data are summarized in Tables 1–3. During the excavation processes, both lining and bolts were placed at a small distance from the facing that can be neglected in the modeling. In what follows, the behavior of the lining is considered as linear elastic and perfect bonding is assumed at the lining/rock and bolt/rock interfaces. Depending on the considered reference, the friction angle of the rock material ranges from um = 26.40 to um = 36.40. As regards the 0 value of parameter k ¼ r0H =r0V , the only reference providing comments on the question is the work by Ward et al. (1983) where it is said that the wide-open vertical joints in the stronger rocks suggest that the horizontal stresses are relatively low compared with the vertical stress. This means that k0 would be smaller than unity (k0 < 1). The mesh used for the finite element numerical simulation is showed in Fig. 5. The geometry domain is discretized into with 6 438 hexahedral 20-nodes elements. In the first simulations, a mean value of the rock friction angle um = 31.40 is adopted together with a ratio of initial horizontal stress to vertical stress smaller than unity k0 = 0.75. The predictions in terms of radial displacement versus distance to the facing are displayed in Figs. 6 and 7. The result corresponding to in situ readings from extensometers placed at the tunnel roof is also represented in Fig. 6. This value corresponds to the maximum displacement measure obtained after stabilization of the tunnel convergence. As it can be observed, both the homogenization method and embedded model lead to very close predictions, which slightly underestimate the measured value of radial displacement.

Fig. 5. Finite element mesh used for the numerical simulations.

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

75

Fig. 6. Radial displacement profile at the roof of the tunnel wall versus distance to the facing. Comparison between homogenized method, embedded method and in situ measurement (um = 31.40, k0 = 0.75).

Fig. 7. Longitudinal profile of radial displacement at the tunnel side wall. Comparison between homogenized method and embedded model (um = 31.40, k0 = 0.75).

Fig. 8. Radial displacements along a vertical axis located at plane z = 6 m and intersecting the tunnel axis. Comparison between the homogenized method and embedded method (um = 31.40, k0 = 0.75).

76

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

Fig. 9. Radial displacement at the roof of the tunnel for um = 26.40.

Fig. 10. Radial displacement of the roof of the tunnel for um = 31.40.

Fig. 11. Radial displacement of the roof of the tunnel for um = 36.40.

77

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

This discrepancy may be attributed to special stochastic variations in the rock material properties as well as to the uncertainty about the value of k0. Fig. 8 shows the profile of radial displacement along a vertical axis located at distance z = 6 m behind the tunnel face. Since there is uncertainty about the values of um and k0 in the available literature, further simulations have been carried out by varying the corresponding values. The predictions obtained by homogenization approach in terms of radial displacement versus the distance to the tunnel face are given in Figs. 9–11. Interestingly, the maximum radial displacement predicted by homogenization approach is an increasing function of the horizontal thrust coefficient. Furthermore, the numerical results suggest the assumption of isotropic confining stress field, i.e., k0 = 1, will always lead to predictions that underestimate the tunnel convergence. This is consistent with the analysis of Ward et al. (1983) that field observations suggest that k0 < 1. Depending the value adopted for um, the ‘‘best’’ value for k0, that is that fit better the measured displacement, lies between 0.5 and 0.75.

can be further extended to various kinds of geotechnical structures reinforced by linear inclusions. In contrast, the embedded model is a numerical technique that basically assumes the same kinematics for both the bolt element and the embedding rock brick element, hence avoiding the numerical difficulties inherent to discrete finite element discretization of bolts. The efficiency of such approach is restricted to the situation of perfect adherence between bolts and the surrounding rock mass. Unlike the embedded method, the homogenization approach can theoretically account for possible slippage at the bolt/rock interface. In connection with the latter aspect, the main issue which remains to be addressed consists in the practical numerical implementation and validation of the algorithm of plastic integration that has been proposed for dealing with the bond-slip between rock and bolts.

5. Conclusion

Regarding the plastic flow rule at the macroscopic level, three distinct situations can occur:

Previous works have already demonstrated the efficiency of the homogenization approach applied to the analysis of deformation in bolt-supported tunnels within a two dimensional setting. This contribution is conceived as an extension to three dimensional configurations of such structures. Regarding the theoretical aspect, the present formulation of the homogenized constitutive equations in elastoplasticity accounts for the bolt/rock interface properties. The multi-potential nature of the plastic flow rule associated with the fact that the yield surface is not smooth required the development of a specific iterative algorithm for plastic integration. The theoretical modeling and associated numerical implementation can readily be used to analyze stresses and strains induced by the tunnel excavation process and reinforcement placement. The numerical tool that has been elaborated allows dealing with either reinforcement by bolts or a lining support. The homogenizationbased approach has proved to be a useful way to capture the strong coupling between the rock mass and the reinforcement, avoiding tedious numerical difficulties that shall arise in a direct finite element approach. In addition, the so-called embedded model, classically used for the numerical simulation of reinforcement in concrete, was successfully implemented in order to assess the accuracy of the numerical predictions based on homogenization method. For illustrative purposes, the proposed approach has been applied to analyze the convergence of the Kielder experimental tunnel for which in situ measurement are available. The numerical predictions are compared to experimental data and show good agreement with readings extensometers. A parametric study that consists in varying some relevant parameters controlling the problem has been undertaken in the context of this analysis. Interestingly, both the homogenization approach and embedded model have led to very similar numerical predictions in the case of Kielder experimental tunnel. Despite the similarity of obtained results, there are however fundamental differences between these approaches. From a theoretical point of view, the constitutive equations formulated by means of homogenization are based upon a consistent mechanical framework, which allows discussing their relevancy in the light model assumptions. Taking advantage of the fully analytical character of the equations that express the homogenized behavior, the influence of relevant parameters on the overall behavior of the bolted rock mass (such as bolt reinforcement scheme and individual properties of the constituents) can readily be undertaken. The numerical tool can therefore be used as an efficient tool to elaborate guidelines for the optimal design of the reinforcement pattern. This clearly points the versatility of the theoretical and associated model whose range of applicability

Appendix A A.1. Some aspects of the macroscopic plastic flow rule

 If 0 ¼ ~f ðRÞ > ~f int ðRÞ then f hom ðRÞ ¼ ~f ðRÞ and consequently the p direction of e_ is controlled by the plastic flow rule of the ‘fictitious’ homogenized material that would result from association of rock and bolts along with perfect bonding at the interface (i.e., the material whose yield function would be ~f Þ. Accordingly, the plastic potential corresponds theoretically to that, say g~, of the ‘fictitious’ homogenized material and which is difficult to determine theoretically. However, if the rock material exhibits a plastic dilatancy angle wm that is different of the internal friction angle /m, the simplest way to defineg~ consists in observing first that the macroscopic plastic function f hom ¼ ~f depends on the friction angle of the rock material, that is f hom ðRÞ ¼ f hom ðR; um Þ ¼ ~f ðR; um Þ. The plastic potential g~ can therefore be defined by substituting in the expression of f hom ¼ ~f the friction angles by the corresponding plastic dilatancy angle:

g hom ðRÞ ¼ g~ðRÞ ¼ f hom ðR; wm Þ

ðA:1Þ

 If ~f ðRÞ < ~f int ðRÞ ¼ 0 then f hom ðRÞ ¼ ~f int ðRÞ and consequently the p direction of e_ is given by the plastic flow rule of the bolt/rock interface. Let us denote by g~int ðRÞ ¼ g int ðTÞ with T ¼ R  n the corresponding plastic potential expressed in terms of stress tensor instead of stress vector. Accordingly,

g hom ðRÞ ¼ g~int ðRÞ

ðA:2Þ

 The last case is that when yielding occurs in both rock and/or bolts (i.e., ~f ðRÞ ¼ 0Þ and at the bolt/rock interface bolts (i.e., ~f int ðRÞ ¼ 0Þ. In other words, the stress sate lies on the corner p or vertex singularity of @Ghom. The direction of e_ is not determined in a unique way. It belongs to the ‘‘cone’’ defined by the normal vectors @ g~=@ R and @ g~int =@ R

e_ p ¼ v_ 1

@ g~ @ g~int ðRÞ þ v_ 2 ðRÞ; @R @R

v_ 1 ; v_ 2 P 0

ðA:3Þ

Appendix B B.1. Iterative procedure for the stress projection onto convex Ghom The two stage iterative procedure aimed at evaluating the projection stress R⁄ onto the convex Ghom and which is defined by Eq. (27) is described below:

78

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79

Step 1: It is first assumed that the projection lies on the intersece \ @G e int . Conditions (30) that characterize the protion line of @ G jection R⁄ can therefore be rewritten as:

8   > < R ¼ R  Ahom : v_ 1 @ g~ ðR Þ þ v_ 2 @ g~int ðR Þ v_ 1 P 0; v_ 2 P 0 @R @R > :~  f ðR Þ ¼ 0; ~f int ðR Þ ¼ 0 ðB:1Þ The unknowns of the problem are the stress projection R⁄ and the Lagrange multipliers ðv_ 1 ; v_ 2 Þ. They are solution of the above system of non-linear equations. In step 2 the derivation of an approximated solution of this system is presented. Step 2: A first-order Taylor expansion is considered for both functions ~f ðR Þ and ~f int ðR Þ:

Cai, Y., Esaki, T., Jiang, Y., 2004. An analytical model to predict axial load in grouted rock bolt for soft rock tunneling. Tunnell. Undergr. Space Technol. 19, 607–618. Crisfield, M.A., 1991. Non-linear finite element analysis of solids and structures. Essentials, vol. 1. Wiley, New York. de Buhan, P., Maghous, S., 1997. Comportement élastique non linéaire macroscopique d’un matériau comportant un réseau de joints. C. R. Acad. Sc. Paris 324 (IIB), 209–218. de Buhan, P., Salençon, J., 1990. Yield strenght of reinforced soils as anisotropic media. In: Boehler, J.P. (Ed.), Yielding Damage and Failure of Anisotropic Media, pp. 791–803. de Buhan, P., Taliercio, A., 1991. A homogenization approach to the yield strength of composite materials. Euro. J. Mech. A/Solids 10, 129–154. de Buhan, P., Mangiavacchi, R., Nova, R., Pelligrini, G., Salençon, J., 1989. Yield design of reinforced earth walls by a homogenization method. Geotechnique 39 (2), 189–201. de Buhan, P., Bourgeois, E., Hassen, G., 2008. Numerical Simulation of boltsupported tunnels by means of a multiphase model conceived as an improved homogenization procedure. Int. J. Numer. Anal. Meth. Geomech. 13 (32), 1597– 1615.

8 ~ ~ hom @ g~ @~f > : @R ðR Þ  v_ 2 @@Rf ðRÞ : Ahom : @@gRint ðR Þ < 0 ¼ ~f ðR Þ ’ ~f ðRÞ  v_ 1 @ R ðRÞ : A ~ ~ > : 0 ¼ ~f int ðR Þ ’ ~f int ðRÞ  v_ 1 @@f Rint ðRÞ : Ahom : @@Rg~ ðR Þ  v_ 2 @@f Rint ðRÞ : Ahom : @@g~Rint ðR Þ

~

Additionally, the value of the derivatives terms @@gRint and @@Rg in (B.2) are evaluated at R instead of R⁄. Introducing the latter approximation in system (32) yields the linear system of equations

2 6 4

3

#

"~ f ð RÞ 7 v_ 1 ¼ 5 ~ ~f int ðRÞ : @gint ðRÞ v_ 2

@ ~f ð @R

RÞ : Ahom : @@g~Rint ðRÞ

@~f ð @R

RÞ : Ahom : @@g~Rint ðRÞ

@ ~f int @R

ðRÞ : Ahom : @@Rg ðRÞ

@~f int @R

ðRÞ : Ahom

~

~

@R

ðB:3Þ The Lagrange multipliers v_ 1 and v_ 2 are obtained by solving system (B.3). The evaluation of the Lagrange multipliers has been based upon the assumption that the projection R⁄ lies on the intere \ @G e int . This implicitly means that v_ 1 and v_ 2 are section line of @ G non-negative quantities. If the obtained values effectively meet this requirement, i.e., v_ 1 P 0 and v_ 2 P 0, an approximation expression for the projection is then given by

R ¼ R  Ahom : v_ 1

! @ g~ @ g~int ðRÞ þ v_ 2 ð RÞ @R @R

ðB:4Þ

In the case when the evaluated values of the Lagrange multipliers are such that v_ 1 > 0 and v_ 2 < 0 or v_ 1 < 0 and v_ 2 > 0, the negative multiplier should be set to zero and the other multiplier computed properly. If for instance v_ 1 > 0 and v_ 2 < 0, which means e of the boundary that the projection R⁄ actually lies on the part @ G of, the Lagrange multipliers are evaluated following a classical way

v_ 1 ’ @~f @R

~f ðRÞ ðRÞ : Ahom : @@Rg ðR Þ ~

;

v_ 2 ¼ 0

ðB:5Þ

and R⁄ is given by introducing (B.5) in (B.4). References Bernaud, D., de Buhan, P., Maghous, S., 1995. Numerical simulation of the convergence of a bolt-supported tunnel through a homogenization method. Int. J. Numer. Anal. Meth. Geomech. 19 (4), 267–288. Bernaud, D., Deudé, V., Dormieux, L., Maghous, S., Schmitt, D.P., 2002. Evolution of elastic properties in finite poroplasticity and finite element analysis. Int. J. Numer. Anal. Meth. Geomech. 26 (9), 845–871. Bernaud, D., Maghous, S., de Buhan, P., Couto, E., 2009. A numerical approach for design of bolt-supported tunnels regarded as homogenized structures. Tunnell. Undergr. Space Technol. 24 (5), 533–546. Bobet, A., Einstein, H.H., 2010. Tunnel reinforcement with rockbolts. Tunnell. Undergr. Space Technol. 26 (1), 100–123.

ðB:2Þ

Deudé, V., Dormieux, L., Kondo, D., Maghous, S., 2002. Micromechanics approach to nonlinear poroelasticity: application to cracked rocks. J. Eng. Mech. ASCE 128, 848–855. Elwi, A.E., Hrudey, T.M., 1989. Finite element model for curved embedded reinforcement. J. Eng. Mech. Div., ASCE 115 (4), 740–754. Freeman, T.J., 1978. The behaviour of fully-bonded rock bolts in the Kieder experimental tunnel. Tun. Tunnel. 5, 37–40. Greuell, E., 1993. Etude du soutènement des tunnels par boulons passifs dans les sols et les roches tendres, par une méthode d’homogénéisation. Ph.D. Thesis, Ecole Polytechnique, Palaiseau, France. Greuell, E., de Buhan, P., Panet, M., Salençon, J., 1994. Behaviour of tunnels reinforced by untensioned bolts. In: Proceedings of the 13th International Conference on Soil Mechanic and Foundation Engineering. New Delhi. Hanafy, E.A., Emery, J.J., 1980. Advancing face simulation for tunnel excavation and lining placement. In: Proceedings of the 13th Canadian Rock Mechanics Symposium on Underground Rock Engineering, Toronto, pp. 119–125. Kovari, K., 2003. History of the sprayed concrete liming method – part II: milestones up to the 1960s. Tunnell. Undergr. Space Technol. 18, 71–83. Lade, P.V., Nelson, R.B., 1984. Incrementalization procedure for elasto-plastic constitutive model with multiple, intersecting yield surfaces. Int. J. Rock Mech. Min. Sci. 8 (4), 311–323. Lunardi, P., 2000. Design and constructing tunnels: A.DE.CO.-RS approach. Tunnels and Tunnelling International, special supplement, May 2000. Maghous, S., de Buhan, P., Dormieux, L., 2002. Non-linear global elastic behaviour of a periodically jointed material. Mech. Res. Commun. 29, 45–51. Maghous, S., Bernaud, D., Fréard, J., Garnier, D., 2008. Elastoplastic behavior of jointed rock masses as homogenized media and finite element analysis. Int. J. Rock Mech. Min. Sci. 45, 1273–1286. Maghous, S., Dormieux, L., Kondo, D., Shao, J.F., in press. Micromechanics approach to poroelastic behavior of a jointed rock. Int. J. Numer. Anal. Meth. Geomech. http://dx.doi.org/10.1002/nag.1087. Oreste, P., 2008. Distinct analysis of fully grouted bolts around a circular tunnel considering the congruence of displacements between the bar and the rock. Int. J. Rock Mech. Min. Sci. 45, 1052–1067. Oreste, P., Peila, D., 1996. Radial passive rockbolting in tunnelling design with a new convergence-confinement model. Int. J. Rock Mech. Min. Sci. Geomech. Abst. 33, 443–454. Panet, M., 1995. Le calcul des tunnels par la méthode convergence-confinement. Presses de l’ENPC, Paris. Peila, D., 1994. A theoretical study of reinforcement influence on the stability of a tunnel face. Geotech. Geol. Eng. 12, 145–168. Peila, D., Oreste, P., 1995. Axisymmetric analysis of ground reinforcing in tunnelling design. Comput. Geotech. 17, 253–274. Peila, D., Oreste, P., Pelizza, S., Poma, A., 1996. Study of the influence of subhorizontal fiber-glass pipes on the stability of a tunnel face. In: North American Tunnelling’ 96, Balkema, ITA/AITES 1996 World Tunnel Congress, Washington, pp. 425–432. Peila, D., Oreste, P., Pelizza, S., 2000. Face reinforcement in deep tunnels. In: Proceedings of the International Underground Congress, Praga, pp. 232–243. Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Springer, Berlin. Sudret, B., de Buhan, P., 2001. Multiphase model for inclusion-reinforced geostructures. Application to rock-bolted tunnels and piled raft foundations. Int. J. Numer. Anal. Meth. Geomech. 25 (2), 155–182. Sudret, B., de Buhan, P., Maghous, S., Bernaud, D., 1998. Elastoplastic analysis of inclusion reinforced structures. Metal Mater. Int. 4 (3), 252–255.

S. Maghous et al. / Tunnelling and Underground Space Technology 31 (2012) 68–79 Suquet, P., 1985. Elements of homogenization for inelastic solid mechanics. In: CISM Lectures Notes. Homogenization Techniques for Composite Media, vol. 272. Springer-Verlag, pp. 93–278. Ward, W.H., TEDD, P., BERRY, N.S.M., 1983. The Kielder experimental tunnel: final results. Géotechnique 33 (3), 257–291. Windsor, C.R., Thompson, A.G., 1993. Rock reinforcement – technology, testing, design and evaluation. In: Hudson, J.A. (Ed.), Comprehensive Rock Engineering, Principles, Practices and Projects, vol. 4. Pergamon Press, Oxford, pp. 451–484.

79

Wullschläger, D., Natau, O., 1983. Studies of the composite system of rockmass and non-prestressed grouted rockbolts. In: International Symposium of Rock Bolting. Abisko, Sweden. Wullschläger, D., Natau, O., 1987. The bolted rockmass as an anisotropic continuum. Material behavior and design suggestions for rock cavities. In: Proceedings of the 6th International Congress of Rock Mechanical. Montreal.