Bolt-Grout Interactions in Elastoplastic Rock Mass Using Coupled FEM-FDM Techniques

Abstract

Numerical procedure based on finite element method (FEM) and finite difference method (FDM) for the analysis of bolt-grout interactions are introduced in this paper. The finite element procedure incorporates elasto-plastic concepts with Hoek and Brown yield criterion and has been applied for rock mass. Bolt-grout interactions are evaluated based on finite difference method and are embedded in the elasto-plastic procedures of FEM. The experimental validation of the proposed FEM-FDM procedures and numerical examples of a bolted tunnel are provided to demonstrate the efficacy of the proposed method for practical applications.

1. Introduction

Rock bolts have been widely used as a primary support system to stabilize the rock masses around tunnel, underground mine galleries, slopes, and others structure made in rock masses. In general, rock bolts reinforce rock masses through restraining the deformation within rock masses [1] and reduces the yield region around the excavation boundary. During the last four decades, different types of rock bolts have been in practice, out of which fully grouted active/passive bolts were the most common types. For a fully grouted passive rock bolt installed in deformable rock masses, a neutral point exists on the bolt rod, where shear stress at the interface between the bolt and grout material vanishes. The pickup length is then defined as the length between free end from the tunnel boundary to the neutral point, and the shear stress along this bolt length drags the bolt towards the tunnel [1]. The bolt length between the neutral point and the other free end of the bolt (inside the rock mass) is designated as anchor length and the developed shear stress drags the bolt towards the rock mass or, in other words, anchors the bolt into the rock mass. Based on these concepts, shear stresses and axial loads developed in a bolt rod are analytically formulated by many researchers who consider the effect of faceplate at the free end of the bolt [1, 2]. Bolt grout interactions around a circular tunnel in Hoek and Brown medium have been formulated analytically considering a bolt density factor [3]. In their study, Indraratna and Kaiser [3] demonstrated that due to the installation of bolts, the Hoek and Brown parameters 𝑚_𝑏, 𝑠, and uniaxial compressive strength of intact rock, 𝜎_𝑐𝑖, would change to new values 𝑚^∗_𝑏, 𝑠^∗, and 𝜎^∗_𝑐𝑖 based on the bolt density factor which consists of bolt spacing, shear stress on the bolt grout interactions, tunnel radius, and bolt diameter. Considering different approaches to bolt performance, Stille [4] presented a closed form elasto-plastic analytical solutions of grouted bolts. Based on shear lag model, Cai et al. [5] derived an analytical solution of rock bolts for describing the interaction behaviours of rock bolt, grout material and rock mass. In addition, recently Torres [6] presented closed-form solution of elastic rock bolt in elasto-plastic rock mass for grouted and anchored rock bolt. In his work, stiffness and strength of grout material are assumed to be infinitely high. Hence, there is virtually no shear displacement in rock bolt interface. Brady and Lorig [7] numerically analyzed the interactions of bolt grout in Mohr Columb medium using finite difference method (FDM) technique and showed that radial displacement, peak tangential stress, and yielded region would reduce due to the installation of grouted bolts around a circular tunnel. FDM technique has also been used for solving differential equation, which governs the interactions between rock bolt and grout material [2].

In addition, numerous studies have been published on the analytical solution of stresses and displacements around a circular tunnel considering elasto-plastic rock mass with Hoek and Brown yield criterion [8–11]. Finite Element Method (FEM) and/or FDM-based procedures are also developed for the analysis of the said problem and have been presented in many references for solving geotechnical problems [12–15].

In this study, FDM technique is applied to evaluate the interactions between bolt and grout considering neutral point, pickup and anchor length concepts. Displacements and stress distributions in the rock mass are computed using FEM technique, considering material nonlinearity based on Hoek and Brown yield criterion. The analysis procedure of bolt grout interactions are evaluated based on FDM and embedded into the elasto-plastic procedure of FEM. The proposed FEM-FDM procedure is validated with the experimental results of pull-out tests conducted by Farmer [16]. It is found that the observed shear strains along the bolt rod are close to those obtained from numerical models. The coupled FEM-FDM technique is then applied to evaluate rock mass behaviour and bolt performance for a circular bolted tunnel. Results are obtained by varying bolt densities and shear stiffness at the bolt grout interface considering the faceplate option.

2. Implementation of Hoek-Brown Yield Criterion in FEM

The generalized Hoek and Brown rock mass failure criterion for drained condition is expressed as follows [17],𝐹(𝜎₁,𝜎₃)=𝜎₁−𝜎₃−𝜎_𝑐𝑖(𝑚_𝑏𝜎₃𝜎_𝑐𝑖+𝑠)^𝑎=0, (1) where 𝜎_1,3 are principal stresses and should be arranged as 𝜎₁≤𝜎₂≤𝜎₃ (considering compressive stress as negative), 𝜎_𝑐𝑖 = uniaxial compressive strength of intact rock and negative value is considered, 𝑚_𝑏 is reduced material constant based on intact property of rock, 𝑚_𝑖, and Geological Strength Index (GSI), and 𝑠 and 𝑎 are rock mass properties based on GSI [18].

For numerical computations, it is convenient to express yield surface using mean hydrostatic stress, σ_m, and deviatoric stress invariants, namely, second invariant of deviatoric stress, 𝐽₂ and Lode angle 𝜃 [13, 18]. This is particularly helpful since computations of principal stresses and their directions in three dimensions can be avoided. In this study, (1) has been reformulated using invariant of stresses as 𝐹(𝜎_𝑚,𝐽₂,𝜃)=0 and is given as 𝐹=𝛼(𝐽₂)^1/𝑎+𝛽𝐽₂+𝑚_𝑏𝜎_𝑚+𝑠𝜎_𝑐𝑖=0, (2) where, the expression for parameters 𝛼=(2cos𝜃)^1/𝑎(−𝜎_𝑐𝑖)^{(𝑎−1)/𝑎},   𝛽=𝑚_𝑏(cos𝜃+(sin𝜃/√3)), and 𝜎_𝑚 = (1/3)𝜎_𝑖𝑖,   𝐽₂ = (1/2)𝑠_𝑖𝑗𝑠_𝑖𝑗,   𝐽₃ = (1/3)𝑠_𝑖𝑗𝑠_𝑗𝑘𝑠_𝑘𝑖,sin3𝜃 = (3√3/2)(𝐽₃/𝐽₂^3/2),   𝑠_𝑖𝑗 = 𝜎_𝑖𝑗−𝜎_𝑚𝛿_𝑖𝑗.

The angle 𝜃 is called Lode angle ranging between ±30° and 𝛿_𝑖𝑗 is the Kronecker delta, whose value is unity when 𝑖=𝑗 and zero when 𝑖≠𝑗. From (2), it is clear that the generalized Hoek and Brown yield function is dependent on all three invariants of stresses. This yield function is plotted in three mutually perpendicular principal stress axes (𝜎₁≤𝜎₂≤𝜎₃) as shown in Figure 1. In principal stress space, the yield locus is conical whose normal section at any point is an irregular hexagon, just as Mohr-Coulomb yield surface [11]. The conical shape of this yield surface is the evidence that hydrostatic stress does influence the yielding. In addition, the apex of the hexagonal pyramid is zero (0) when mean hydrostatic stress is 𝜎_𝑚=−𝑠𝜎_𝑐𝑖/𝑚_𝑏.

Figure 1 

Generalized Hoek and Brown yield surface in principal stress space.

2.1. Properties of the Yield Surface

Figure 2 shows the two dimensional deviatoric or 𝜋 plane (𝜎₁+𝜎₂+𝜎₃=0) represent of the generalized Hoek and Brown yield surface. The yield surface is convex and has six sharp corners or points of singularity where tangent of yield locus is indeterminate. In all other points, the yield surface is regular. The singular points occur when 𝜃 approaches to ±30 degrees. As a result, the direction of plastic straining is also indeterminate and thus elasto-plastic constitutive matrix cannot be determined at those points on the yield surface.

Figure 2 

Yield surfaces in deviatoric plane.

2.2. Flow Rule

A nonassociative flow rule assumes that the normality condition, the direction of plastic strain increment is normal to the yield surface, is not satisfied. Hence an alternative plastic potential function, 𝑄=𝑄(𝜎_𝑚,𝐽₂,𝜃), must be assumed, such that 𝑄≠𝐹. On the contrary, associative flow rule assumes that 𝑄=𝐹. The linear plastic potential function in 𝜎₁-𝜎₃ plane as given in (3a) has been considered in conjunction with Hoek and Brown yield criterion for the analysis of ground reaction curve and elasto-plastic analysis of circular tunnel [8, 9]. Equation (3a) is modified to obtain the expression for plastic potential function using mean stress and deviatoric stress invariants, as given in (3b). In deviatoric plane, this plastic potential surface is convex and resembles an irregular hexagon, except when zero dilation angle, it is regular, having six points of singularity when 𝜃 approaches ±30 degrees,𝑄=−𝜎₁+𝜎₃𝑁_𝜓,𝑄=𝐽₂[(1+𝑁_𝜓)cos𝜃+(𝑁_𝜓−1)sin𝜃√3]+(𝑁_𝜓−1)𝜎_𝑚=const, (3a)(3b)where 𝑁_𝜓=(1+sin𝜓)/(1−sin𝜓) and 𝜓 = dilation angle.

2.3. Elastic-Plastic Analysis in Plane Strain Condition

The elasto-plastic constitutive matrix for elastic-perfectly plastic (strain hardening parameter taken to be zero) is represented by the following equation [12]:𝐃_𝑒𝑝=𝐃−𝐃{𝜕𝑄}𝜕𝝈{𝜕𝐹/𝜕𝝈}^𝑇𝐃{𝜕𝐹/𝜕𝝈}^𝑇𝐃{𝜕𝑄/𝜕𝝈}, (4) where for plane strain condition, stresses and strains are expressed in Cartesian coordinate system as {𝝈}^𝑇={𝜎_𝑥𝜎_𝑦𝜏_𝑥𝑦𝜎_𝑧}, {𝜺}^𝑇={𝜀_𝑥𝜀_𝑦𝛾_𝑥𝑦0}. The matrix D is the elastic constitutive matrix estimated based on elastic modulus, 𝐸, and Poisson’s ratio, 𝜈, which gives𝐸𝐃=(1+𝜈)(1−2𝜈)1−𝜈𝜈0𝜈𝜈1−𝜈0𝜈001−2𝜈20𝜈𝜈01−𝜈, (5) For elastic stress increment, first three rows and columns are considered. Representing {𝜕𝐹/𝜕𝝈}=𝐶₁𝐚₁+𝐶₂𝐚₂+𝐶₃𝐚₃ and {𝜕𝑄/𝜕𝝈}=𝐷₁𝐚₁+𝐷₂𝐚₂+𝐷₃𝐚₃, for plane strain condition, a_i (i = 1,…, 3) can be obtained in Cartesian coordinate system as 𝐚^𝐓₁=𝜕𝜎_𝑚=1𝜕𝜎3,𝐚1101_𝐓2=𝜕𝐽₂=1𝜕𝜎2𝐽₂𝑠_𝑥𝑠_𝑦2𝜏_𝑥𝑦𝑠_𝑧,𝐚^𝐓₃=𝜕𝐽₃𝜕𝜎={(𝑠_𝑦𝑠_𝑧+𝐽₂3)(𝑠_𝑥𝑠_𝑧+𝐽₂3)−2𝑠_𝑧𝜏_𝑥𝑦(𝑠_𝑥𝑠_𝑦−𝜏²_𝑥𝑦+𝐽₂3)}, (6) where 𝐽₂=(1/2)(𝑠²_𝑥+𝑠²_𝑦+𝑠²_𝑧)+𝜏²_𝑥𝑦 and 𝐽₃=𝑠_𝑧(𝑠²_𝑧−𝐽₂).

The constants 𝐶_𝑖 and 𝐷_𝑖 are the derivatives of yield and plastic potential functions, respectively, with respect to 𝜎_𝑚, 𝐽₂ and 𝐽₃ and given in following equations: 𝐶₁=𝑚_𝑏,𝐶₂=𝑅₂−tan3𝜃𝐽₂𝑅₁,𝐶₃=√3𝑅₁2cos3𝜃(𝐽₂)³,𝑅₁=−(2sin𝜃𝑎)(2cos𝜃)^{(1−𝑎)/𝑎}(−𝜎_𝑐𝑖)^{(1−𝑎)/𝑎}(𝐽₂)^1/𝑎−𝑚_𝑏(−sin𝜃+cos𝜃√3)𝐽₂,𝑅₂=𝛼𝑎(𝐽₂)^{(1−𝑎)/𝑎}𝐷+𝛽,₁=𝑁_𝜓−1,𝐷₂=𝑇₂−tan3𝜃𝐽₂𝑇₁,𝐷₃=√3𝑇₁2cos3𝜃(𝐽₂)³,𝑇₁=𝐽₂[−(1+𝑁_𝜓)sin𝜃+(𝑁_𝜓−1)cos𝜃√3𝑇],₂=[(1+𝑁_𝜓)cos𝜃+(𝑁_𝜓−1)sin𝜃√3]. (7) Finite element implementation of the above yield functions is performed, based on the procedure outlined in various literatures [12, 14, 19]. In this case, Newton-Raphson solution method is applied with initial stiffness conditions. Hence, large numbers of load steps are taken for satisfying convergence criterion. The increments of plastic strain and stress are estimated using the following equations:{Δ𝜀}^𝑝=𝜆{𝜕𝑄𝜕𝝈},Δ𝝈^𝑡_𝑖=𝐃_𝑒𝑝{Δ𝜀}, (8a)(8b)where 𝜆= plastic multiplier

3. Analysis of Bolt-Grout Interactions Using FDM

A bolt installed in deformable rock mass is subjected to axial loading and it provides resistance to the movement of rock mass through shear stresses developed axially in the bolt-grout interfaces. The bolt is represented by a one-dimensional bar and hence the relationship between the shear force and the axial load can be established considering a small (differential) section of the bolt as shown in Figure 3. The force equilibrium in the axial direction leads to the following expression: 𝐴_𝑏𝜕𝜎_𝑥𝑥𝜕𝑥𝑑𝑥=−𝐹_𝑠𝑑𝑥 (9) or 𝑑𝜎_𝑥𝑥𝐹𝑑𝑥=−^𝑠𝐴_𝑏, (10) where 𝐹_𝑠= shear force due to bond per unit length, 𝐴_𝑏= is the area of the cross-section of the bolt =𝜋𝑟²_𝑏, 𝐿 = length of the bolt, 𝜎_𝑥𝑥= axial stress in the bolt, and 𝑟_𝑏= radius of the bolt rod.

Figure 3 

Stress component along a fully grouted bolt.

Now, for an elastic bolt with elastic modulus of 𝐸_𝑏, 𝜎_𝑥𝑥=𝐸_𝑏𝑑𝑢_𝑥𝑑𝑥. (11) Hence, substituting (11) into (10) gives 𝑑²𝑢_𝑥𝑑𝑥²𝐹=−^𝑠𝐴𝐸_𝑏. (12) Farmer [16] proposed and justified that shear force developed due to bond between bolt and grout for a unit length of bolt. This relationship is a linear function (presumed) of the relative movement between rock and bolt [2, 20]. The relationship is expressed as 𝐹_𝑠=2𝜋𝑟_𝑏𝑘(𝑢_𝑟−𝑢_𝑥)=𝑘^∗(𝑢_𝑟−𝑢_𝑥), (13) where 𝑘 is the shear stiffness of bolt grout interface and can be determined from stress-displacement response during laboratory pullout tests. In general, 𝑘^∗ depends on the bolt and rock displacements due to debonding between bolt rod and grout. Experimental studies of pull-out tests conducted by Farmer [16] showed that load-displacement curve is non-linear, which is evident from de-bonding between bolt rod and grout. However, in this present study, de-bonding or elasto-plastic behaviour of grout material has not been considered and hence shear stiffness is kept constant in (13). For yielding of bolt steel material, Tresca criterion has been implemented. In the example problem, given in the later part of the paper, the axial stresses are found to be far less than the yield strength of bolt steel and hence bolts do not yield for any condition. As a result, theoretical background of plasticity of bolt steel is omitted.

Displacement function 𝑢_𝑟(𝑥) describes the movement of rock mass over bolt length (0<𝑥<𝐿). In general, 𝑢_𝑟 will be a monotonically decreasing function with radial distance r, measured from a circular tunnel boundary. The form and rate of decrease of 𝑢_𝑟 will depend on, namely. the size and shape of the opening, the strength and structure of the rock mass, stress redistribution occurring around excavation due to the yield of rock mass, and number of bolts installed at the excavation boundary.

Now, combining (12) and (13), the distribution of displacement along bolt length is described by the second order inhomogeneous linear ordinary differential equation as 𝑑²𝑢_𝑥𝑑𝑥²−𝑘^∗𝐴𝐸_𝑏𝑢_𝑥𝑘=−^∗𝐴𝐸_𝑏𝑢_𝑟. (14) Equation (14) suggests that the solution of 𝑢_𝑥 can be obtained if 𝑢_𝑟 is known beforehand. In this study, the solution 𝑢_𝑟 is obtained from the nodal displacements of a FEM model with bolts installed at the excavation boundary.

For this purpose, bolt is discretized along its length for implementation of finite difference scheme as shown in Figure 4. Hyett et al. [2] provided the basic formulation of (14) using the finite difference method. In the following, the significant part of that concept is mentioned for implementation of (14) in FEM procedure. In discretized form, (14) can be written with central difference form as (𝑢_𝑥^𝑖+1−𝑢^𝑖_𝑥)/ℎ^𝑖+1−(𝑢^𝑖_𝑥−𝑢_𝑥^𝑖−1)/ℎ^𝑖(ℎ^𝑖+1+ℎ^𝑖−𝑘)/2𝐴𝐸^∗_𝑏𝑢^𝑖_𝑥𝑘=−^∗𝐴𝐸_𝑏𝑢^𝑖_𝑟, (15) where the superscript i refers to the 𝑖th discretized point on the bolt and ℎ^𝑖 denotes step length as ℎ^𝑖=𝑥^𝑖−𝑥^𝑖−1. Equation (15) can be rewritten as−𝐴𝐸_𝑏ℎ^𝑖𝑢_𝑥^𝑖−1+(𝐴𝐸_𝑏ℎ^𝑖+𝐴𝐸_𝑏ℎ^𝑖+1+𝑘^∗(ℎ^𝑖+1+ℎ^𝑖2))𝑢^𝑖_𝑥−𝐴𝐸_𝑏ℎ^𝑖+1𝑢_𝑥^𝑖+1=𝑘^∗(ℎ^𝑖+1+ℎ^𝑖2)𝑢^𝑖_𝑟. (16) Equation (16) can be rearranged in more compact form as −𝐵^𝑖−1𝑢_𝑥^𝑖−1+(𝐵^𝑖−1+𝐵^𝑖+𝑘^∗(ℎ^𝑖+1+ℎ^𝑖2))𝑢^𝑖_𝑥−𝐵^𝑖𝑢_𝑥^𝑖+1=𝑘^∗(ℎ^𝑖+1+ℎ^𝑖2)𝑢^𝑖_𝑟, (17) where 𝐵^𝑖=𝐴𝐸_𝑏/ℎ^𝑖+1.

Figure 4 

Finite difference scheme and notation.

Now, using force-free boundary condition for the free end 𝑖=1 (without faceplate bolt), the first equation transforms:(𝐵₁+𝑘^∗(ℎ²2))𝑢¹_𝑥−𝐵₁𝑢²_𝑥=𝑘^∗𝑢_𝑟0(ℎ²2). (18) For a faceplated bolt, the displacements of rock mass and bolt are same at 𝑖=1. Thus, reaction force is developed at that end. Then the unknown reaction force 𝐹¹_𝑥 is determined as−𝐹¹_𝑥−𝐵¹𝑢²_𝑥=−𝐵¹𝑢¹_𝑟. (19) It may be noted that installation of faceplate ensures that the entire length of bolt acts as anchor length and the maximum tensile stress occurs at the tunnel boundary end.

For a passive grouted bolt, the other end of bolt, for both with or without faceplated, 𝑛th nodal point is free. Hence, the nth equation becomes −𝐵^𝑛−1𝑢_𝑥^𝑛−1+(𝐵^𝑛−1+𝑘^∗(ℎ^𝑛2))𝑢^𝑛_𝑥=𝑘^∗(ℎ^𝑛2)𝑢^𝑛_𝑟. (20) Equations (16), (19), and (20) corporate 𝑛 linear simultaneous equations, which become a tri-diagonal matrix with n unknown values of 𝑢_𝑥. The matrix can be solved directly by Thomas algorithm or Gauss elimination method. In addition, axial tensile load on bolts is calculated from the equation 𝐹=𝐴𝐸_𝑏𝑑𝑢_𝑥𝑑𝑥 (21a) or in discretized form of above equation becomes 𝐹^𝑖=𝐴𝐸_𝑏(𝑢_𝑥^𝑖+1−𝑢^𝑖_𝑥𝑥^𝑖+1−𝑥^𝑖). (21b)The shear force at bolt grout interface is obtained from (13).

4. Implementation of Bolt-Grout Interactions into Elasto-Plastic FEM Procedures

It is assumed that the shear force developed at the bolt-grout interface, resists the movement of rock mass along bolt axis direction. In order to implement this concept numerically, shear force at discretized points (FDM) is transformed to nodal reaction forces (FEM) of a rock element and is added to the reaction forces generated by that element. These reaction forces are applied into the next iteration in Newton-Raphson solution scheme (see (22)). In the following section, this concept is discussed in detail. Given𝐊^𝑡_𝑖Δ𝐪^𝑡_𝑖=𝐅^𝑡−[𝐑^𝑡_𝑟𝑖+𝐑^𝑡_𝑏𝑖], (22) where 𝐊^𝑡_𝑖 global stiffness matrix at load step 𝑡 and iteration i = ∑^𝑛𝑒_𝑒=1∫_𝑉𝑒𝐁^𝐓𝐃_𝑒𝑝𝐁𝑑𝑣, 𝐃_𝑒𝑝=𝐃 if yielding of rock does not occur at the Gauss point. The above integration is performed based on Gauss-Legendre quadrature rule.𝑉_𝑒:volume of the element e;𝑛_𝑒:number of elements;▵𝐪^𝑡_𝑖:incremental displacement vector;𝐅^𝑡:applied load at load step 𝑡;𝐑^𝑡_𝑟𝑖:reaction force vector (rock element) at load step 𝑡 at iteration 𝑖−1;𝐑^𝑡_𝑏𝑖:reaction force vector (bolt element) at load step 𝑡 at iteration 𝑖−1.

Let us consider a 4-noded rock element is bisected by a bolt element as shown in Figure 5. The discretized points of the bolt element are i and j, respectively, and numerals 1 to 4 denote the node numbers of the rock element. It may be noted that points i and j are automatically calculated as the intersection points between bolt and edges of solid elements. These points are taken as nodal points for FDM procedure. From (13), the shear force per unit length at points i and j can be estimated after the bolt displacements are computed using FDM procedures. Let 𝐹^𝑆_𝑖 and 𝐹^𝑆_𝑗 denote the shear forces per unit length at points i and j respectively. If the values of shape functions at the midpoint (denoted by m) between the points i and j are 𝑁_𝑘(𝑘=1,…,4), then shear force acting on the bolt segment 𝑙_𝑖𝑗 can be transferred to nodal reaction forces as𝑅_𝑘=−𝑁_𝑘(𝐹^𝑆_𝑖+𝐹^𝑆_𝑗)2𝑙_𝑖𝑗for𝑘=1,…,4. (23)

Figure 5 

Discretization rock and bolt interaction form.

Nodal reaction force is then transformed into global coordinate system as 𝑅_𝑘𝑥=𝑅_𝑘cos𝜃 and 𝑅_𝑘𝑦=𝑅_𝑘sin𝜃. In (23), the negative sign denotes shear force acting on rock mass. The Newton-Raphson method is applied to obtain solution for displacements. In this procedure, for each load step and for every iteration within a load step, Hoek-Brown yield criterion (see (2)) is evaluated at each Gauss point of a rock element with the updated trial stress, as follows.𝝈^𝑡_𝑖=𝝈^𝑡_𝑖−1+Δ𝝈^𝑡_𝑖. (24) If the criterion is satisfied, plastic correction is applied to estimate the increment of plastic strain and corrected stress level (8b). If the previous state of stress was in elastic domain then transition from elastic to plastic regime is estimated as elasto-plastic ratio (a value between 0 and 1.0) for correction of total strain increment in plastic regime. With the updated stress (24), nodal reaction force for each element is computed using 𝑅^𝑡_𝑟𝑖=_𝑉𝑒𝐁^𝑇𝝈^𝑡_𝑖𝑑𝑉. (25) For each bolt element, bolt displacements are computed at each discretized points using (17). Then, shear force per unit length is obtained at each point using (13). Shear force is transformed into nodal reaction forces using (23) and added to the above nodal reaction forces. The global internal out-of-balance force vector, as given in right-hand side of (22) is applied as input parameter to N-R procedure for the next iteration. In this way, iterations are continued within a load step until convergence criterion is satisfied. Then load step is incremented and the same process continues until the total load is applied. The algorithm for coupled FEM-FDM procedure for bolt-grout interactions in elasto-plastic rock mass is given in Algorithm 1.

Thus, for 𝑡th load step  𝑖th iteration, do the following:    For each Gauss (sampling) point of an rock element       Initialize 𝐃𝑒𝑝=𝐃       Find the trial stress Δ𝝈𝑡𝑖=𝝈𝑡𝑖−1+𝐃𝑒𝑝Δ𝜺𝑡𝑖       Using trial stresses       Check if 𝐹>0       (see (3a) and (3b))       If “No” trial stress is the updated stress       If “Yes” plastic correction is required (see (5) and (9))       Determine nodal reaction forces (see (26a) and (26b))    End For loop    For each bolt element       Use FDM procedure to compute displacements       at each discretized points on the bolt (𝑢𝑥) (see (18))       Determine shear force at each discretized point (see (14))       Transform the shear forces into nodal reaction       forces for the rock elements having a bolt (see (24)) Add this reaction force to the above reaction force    End For Loop  End loop iteration End loop load step

Algorithm 1 

5. Experimental Validations

The proposed FEM-FDM procedure has been validated using the results obtained from “pull-out” tests of rock bolts conducted by Farmer [16]. In general, the “pull-out” test is intended to measure the short-term strength of rock bolt anchor [21]. The strength is measured by a pull test in which bolt-head displacement is measured as a function of the applied load. In his pioneering work, Farmer [16] conducted pull-out tests of resin grouted rock bolts anchored in elastic concrete, limestone, and chalk. In his study, grouted lengths of bolts were 350 and 500 mm in concrete. The tests were conducted after 24 hours from the installation of the bolts. Five strain gauges were installed on the rod at regular interval along the length of the bolt. During testing, the anchors were loaded against surface bearing plate of 300 mm², at a force rate of 5 kN/min considering conditions of rigid boundary. Farmer [16] also derived an ordinary differential equation of bolt displacement, similar to (14), ignoring the displacement of rock or concrete.

Here, a concrete block of 1.0 m × 1.5 m has been discretized with 4 noded quadrilateral elements of 100 mm size and grouted rock bolt of length 500 mm is loaded axially, as shown in Figure 6. In pull-out tests, since axial load is applied only on the bolt rod, axial force of the first discretized point (i = 1) on bolt can be rewritten using (21a) and (21b) as 𝐹=𝐴𝐸_𝑏(𝑢²_𝑥−𝑢¹_𝑥ℎ²), (26a) or𝐵¹𝑢¹_𝑥−𝐵¹𝑢²_𝑥=−𝐹¹_𝑥. (26b)

Figure 6 

Discretized model for pull-out test.

In (26b), the axial force 𝐹¹_𝑥 (pull-out load) is a known parameter. It is incorporated into (17) for i =1 keeping all other conditions same. shear strain along the bolt is determined using (27) and compared with those of experimental observations recorded at gauge stations,𝛾_𝑥=𝑑(𝑢_𝑟−𝑢_𝑥)𝑑𝑥. (27) The initial shear stiffness of bolt grout interface 𝑘^∗ is calculated from load-displacement curve recorded by Farmer [16] and it is found to be 6 GN/m per unit length. In numerical models, the elastic moduli of concrete, and bolt rod are considered to be 20 GPa and 180 GPa, respectively, as given by Farmer [16]. The bolt and hole diameters are 20 mm and 28 mm, respectively.

Farmer [16] had recorded the gauge strains for different pull-out loads ranging between 20 kN and the breaking load. Shear strains along bolt length of the pull-out test are shown in Figure 7. Farmer [16] mentioned that experimental shear strain/stress distribution for pull force of 20 kN is close to the theoretical distribution which he had proposed. However, he observed substantial differences in shear stress at 40 kN pull force and this difference increased with increased force. He concluded that de-bonding between bolt rod and grout at higher level of pull force might be the reason for the differences and it was also evident from the non-linear load-displacement curve.

Figure 7 

Distributions of strain along the bolt length.

Numerical study conducted on the same elastic rock medium suggests that during pull-out load, the displacement of rock medium (𝑢_𝑟) is negligible as compared to bolt displacement (𝑢_𝑥). Hence, ignoring 𝑢_𝑟, (27) provides the axial strains which are similar to the measured values. In this study, shear strains obtained from numerical models are compared with those of experimental results of different pull-out loads. Shear strains are estimated at discretized points along bolt length. It is found that up to 80 kN pull-out force, shear strain distribution along the bolt rod is very close to the values obtained from the experimental results (Figure 7). However, numerical results differ from those of experimental observations if the pull-out force exceeds 80 kN. The difference in results are contributed due to de-bonding between bolt rod and grout which is not considered in FEM-FDM procedure (𝑘^∗ being constant). From this study, it is clear that FEM-FDM procedure proposed in this paper simulates pull-out test results and hence can be applied to other applications, one of which is elaborated in the next section.

6. Numerical Examples of a Circular Tunnel Having Rock Bolts

A circular tunnel made in rock mass, having 4 m in diameter, is considered for numerical studies using the above FEM-FDM procedures. The finite element and finite difference codes are developed by the authors and are used for the numerical computations of the above problem. The rock mass is subjected to yield, based on Hoek and Brown criterion, as given in (1). Bolts are installed around the tunnel boundary (tunnel and bolt-grout parameters are given in Table 2). The geometry, finite element mesh and boundary conditions are shown in Figure 8(a). Finer mesh is developed near the tunnel boundary for improving the accuracy of displacement and stress distributions. Since the problem is axisymmetric, it is sufficient to analyze only one quarter of the tunnel. The models are made up with 352 quadrilateral elements having 405 nodes and are analyzed in plane strain condition. It is also considered that faceplate may be installed at the free end of the bolt (tunnel boundary side). Bolt density and shear stiffness at bolt grout interface are also varied to evaluate the performance of bolts and response of the rock mass for a uniform far field stress of 8.0 MPa.

Figure 8 

(a) Meshed tunnel with uniform far field stress, (b) zoomed view near tunnel boundary (bold lines signify 2 m bolts).

Figure 8(b) shows the enlarged view near the tunnel boundary depicting the 9 installed bolts. It may be noted that, bolts can be installed at any part of rock element. However, in this study, they are installed along the boundary of the rock elements for symmetry. Material properties of intact rock and Hoek-Brown parameters are listed in Table 1. It may be seen that residual properties of Hoek and Brown parameters are kept same as the peak properties. Results of bolt models are compared with those of no-support model. Again, results of no-support model are compared with analytical solutions having elastic perfectly plastic conditions in Hoek-Brown media, Hence, in this study, all numerical models are analyzed considering elastic perfectly plastic conditions. The initial stiffness scheme is adopted for solving displacements using Newton-Raphson method.

Table 1 

Rock mass properties used in following discussion.

Table 2 

6.1. Results and Discussions

In the following subsections, stress and displacement distributions around the tunnel boundary, and load and shear response of bolts are discussed in details. Influence of faceplate on the bolt performance, significance of different shear stiffness of bolt grout interface, and effect of bolt density on rock mass are mainly focused in the discussion. For comparison purpose, a model without bolts is also analyzed and the corresponding displacements and stresses in the rock mass are plotted with those bolt models. From here on, the model without support will be termed as “no-support” model. The models having and not having faceplated bolts will be called as “with faceplate” and “without faceplate” models, respectively.

6.1.1. Influence of Faceplated Bolts in Rock Mass Behaviour

Figure 9 shows the radial displacement of rock mass with radial distance from tunnel boundary. It can be seen that the peak radial displacement is reduced by about 15% and 25% for without and with faceplate models, respectively, as compared to no-support model. In other words, faceplated bolts reduce the peak radial displacements by about 11% as compared to bolts without faceplate. The influence of faceplate can also be observed in the stress distribution plots shown in Figure 10. It is clear that the extent of yield region around the tunnel has decreased by at least 0.25 m from no-support model due to the installation of bolts. The peak tangential stress is lower and the overall radial (confining) stress is higher in the rock mass, if faceplate is installed in the bolts. These results suggest that bolts especially with faceplates can provide effective support to the yielding rock mass around a tunnel boundary, provided sufficient shear stiffness exists between bolt and grout interface.

Figure 9 

Distributions of radial displacements in rock mass around the tunnel boundary.

Figure 10 

Distributions of stresses in rock mass around the tunnel boundary.

6.1.2. Effect of Faceplate on Bolt Performance

Figure 11 depicts the bolt load and shear load per unit length obtained from with and without faceplate models. The maximum tensile load carried by a faceplated bolt is almost twice than that of without faceplated bolt. The maximum shear load in the anchored portion of the bolt length is found to be 3.5 times in faceplated bolt as compared to nonfaceplated bolts. Moreover, for non-faceplated bolts, the anchor length is about 1.9 m as compared to 2.0 m for faceplated bolts. These results clearly demonstrate the need for faceplated bolts for enhancing the performance of a 2 m rock bolt.

Figure 11 

Distributions of load and shear stress on bolt rod (with and without faceplate).

6.1.3. Significance of Shear Stiffness of Bolt-Grout Interface on Rock Mass Behaviour

As the shear stiffness at the bolt grout interface increases, overall stresses and displacements around the tunnel boundary reduce, due to increase in shear resistance along the bolt length. In this study, three shear stiffness values are considered, namely. 𝑘=1, 10, and 100 GPa/m. Shear stiffness of 1 GPa/m may represent a weak grout mix whereas 100 GPa/m may represent very strong grout (probably resin mix). If shear stiffness is increased from 1 GPa/m to 100 GPa/m, the radial displacement at the boundary of the tunnel is reduced by 22% for models with faceplate bolted tunnel (Figure 12). Similar trend is observed in models without faceplate. It is also found that the faceplated bolts have significantly reduced the rock movement all along the bolt length and thus the stability of the rock mass have increased around the tunnel boundary.

Figure 12 

Effect of shear stiffness on distributions of displacements in rock mass along the bolt length (with faceplate).

Stress distributions around the tunnel boundary are shown in Figure 13 for with faceplate models. As the shear stiffness increases, the radius of elasto-plastic zone decreases. For these models having 9 bolts installed equidistant from each other, it is found that the radius of elasto-plastic zone (measured from the centre of the tunnel) is reduced from 2.438 m to 2.313 m once 𝑘 is increased from 1 GPa/m to 100 GPa/m. These results emphasize upon bolts as being primary support system would be effective if the bond (shear stiffness) between the bolt and grout is sufficiently stronger.

Figure 13 

Effect of shear stiffness on stress distributions in rock mass (with faceplate).

6.1.4. Effect of Shear Stiffness on Bolt Performance

As before, bolt performance is analyzed in terms of bolt load and shear load developed in a rock bolt. As expected, grout having stronger bond strength or higher shear stiffness produces higher axial load and shear load in the bolt grout interface providing better support against the rock movement. The peak shear load is 0.042 MN/m for 𝑘=1 GPa/m and 0.9217 MN/m for 𝑘=100 GPa/m. Hence, the peak anchoring shear load is about 22 times higher for 𝑘=100 GPa/m as compared to 𝑘=1 GPa/m for faceplated bolts (Figure 14). In the same token, the peak axial tensile load increases by 290 kN for the same conditions. These results are the manifestation for the need of higher bond strength (shear stiffness) grout material, for better performance of bolts.

Figure 14 

Performance of bolt with varying shear stiffness (with faceplate).

6.1.5. Influence of Bolt Density on Rock Mass Behaviour

As mentioned earlier, bolt density is defined as the number of bolts installed per unit area of the tunnel boundary. Mathematically, for a circular tunnel of radius 𝑎, bolt density, 𝛽, is expressed as𝑛𝛽=𝑆_𝑇𝑆_𝐿 (28a) or 𝑛𝛽=2𝜋𝑎𝑆_𝐿, (28b)where 𝑛 is number of bolts, 𝑆_𝑇=2𝜋𝑎 = Transverse distance (m), 𝑆_𝐿= Longitudinal distance between two bolts (m).

In this case, since quarter tunnel is considered, the above equation can be expressed as 𝛽=(𝑛−1)/𝜋 for 𝑎=2 m and 𝑆_𝐿=1 m. For maintaining symmetry, 2, 3, 5, and 9 bolts are included in the numerical models corresponding to 4, 8, 16, and 32 bolts for the entire tunnel boundary, respectively. It is obvious that, as the number of bolt increases, value of 𝛽 also increases. For no-support condition, 𝛽 is considered to be 0.

In general, if the frequency of bolting (bolt density) is increased then the displacements, the peak tangential stress, and the radius of elasto-plastic boundary may reduce. On the contrary, the overall radial stress will increase with increasing values of 𝛽. Results obtained from numerical models support the above concepts as mentioned below. Figure 15 depicts the variation in the peak radial displacement at tunnel boundary with bolt density. From this figure, it is clear that (i) the peak radial displacement decreases with bolt density or number of bolt installed, and (ii) bolts with faceplate performs better as compared to bolts without faceplate as far as reduction in radial displacement is concerned. Figures 16 and 17 show the distribution of radial displacements for different bolt densities for faceplated and without faceplated model. It is found that in faceplated bolt model, the peak normalized radial displacement in percentage is reduced from 0.225 to 0.175 for 𝑘=100 GPa/m, if 𝛽 is increased from 0.0 to 2.55. Similar reduction in radial displacement is observed for without faceplated bolt model with increased bolt densities (Figure 16).

Figure 15 

Variations in radial displacement in rock mass at tunnel boundary with bolts density.

Figure 16 

Distributions of radial displacements in rock mass (without faceplate) with bolt density.

Figure 17 

Distributions of radial displacements in rock mass (with faceplate) with bolt density.

6.1.6. Effect of Bolt Density on Bolt Performance

Axial tensile loads on the bolt are plotted with different bolt densities for without and with faceplate models (Figures 18 and 19, resp.). Load on the bolts has marginally decreased with increasing 𝛽. This suggests that as the number of bolts increases, loads are distributed equally among all installed bolts. It may also be noticed that with faceplate, bolts carry almost twice the load as compared to without faceplate bolts, and hence the performance of the former is superior than the later.

Figure 18 

Distributions of load on bolt along bolt length (without faceplate).

Figure 19 

Distributions of load on bolt along bolt length (with faceplate).

7. Conclusions

This paper proposes a methodology for the analysis of interaction between rock bolt and grout material and behaviour of rock mass using coupled FEM and FDM procedures. The FDM procedures are applied for evaluation of displacements, axial load of a bolt, and shear load per unit length at the bolt grout interface. The developed shear load at the interface is transformed as nodal reaction forces of rock element in FEM analysis. The FEM-FDM procedure has been validated with results obtained from pull-out tests. These procedures have been implemented in elasto-plastic finite element analysis using Hoek-Brown yield criterion. The developed technique is then applied for analyzing a bolted circular tunnel with variations in shear stiffness in the bolt-grout interface, bolt density, and faceplate options. This study shows the efficacy of the proposed method for practical applications. The specific conclusions of this study are summarized as follows.(i)The performance of faceplated bolts in terms of axial tensile load and shear load in the interface is superior to non-faceplated bolts. The study finds that for faceplate bolts the peak radial displacements in rock mass have reduced by 11% and 25% from in place of non-faceplated bolts and no-support conditions, respectively.(ii)The peak tangential stress is lower and the overall radial stress (confining stress) is higher if faceplates are attached in bolts.(iii)The radius of elasto-plastic zone reduces by 0.25 m from no-support conditions if bolting is considered along the 2 m radius tunnel boundary.(iv)The shear stiffness (bond strength grout material) plays an important role in bolt performance. The study finds that the peak axial tensile load increases from 110 to 400 kN if shear stiffness increases from 1 to 100 GPa/m. Similarly, shear load per unit length of the bolt also increases 22 times if the same condition exists.(v)Bolt density (frequency of bolting) around the tunnel also has an effect on the rock mass movements. It is found that the peak radial displacement in rock mass reduces by about 1.08 mm (with faceplated) and by 1.0 mm (without faceplated) from no-support condition if the number of bolts in a quarter tunnel of 2 m radius is increased from 0 to 9.