Study on Mechanical Characteristics of Fully Grouted Rock Bolts for Underground Caverns under Seismic Loads

This study establishes an analytical model for the interaction between the bolt and surrounding rock based on the bearing mechanism of fully grouted rock bolts. The corresponding controlled differential equation for load transfer is deduced. The stress distributions of the anchorage body are obtained by solving the equations. A dynamic algorithm for the bolt considering shear damage on the anchoring interface is proposed based on the dynamic finite element method. The rationality of the algorithm is verified by a pull-out test and excavation simulation of a rounded tunnel. Then, a case study on the mechanical characteristics of the bolts in underground caverns under seismic loads is conducted. The results indicate that the seismic load may lead to stress originating from the bolts and damage on the anchoring interface. The key positions of the antiseismic support can be determined using the numerical simulation. The calculated results can serve as a reference for the antiseismic optimal design of bolts in underground caverns.


Introduction
Southwestern China is a region characterized by high mountains, rolling hills, and abundant water resources.Many hydropower stations have been established or are under construction in this region.Restricted by geological conditions, these hydropower stations are mainly built across underground rocks, and many large-scale underground caverns have been formed.The region is located at plate boundaries with a high probability of earthquakes, and its seismic fortification intensity is generally above VII.The antiseismic performance of underground caverns is directly related to the safe operation of hydropower stations.With the application of the New Austrian Tunneling Method, bolt support has been the most widely used as a flexible reinforcement measure.Bolting static design theory has been increasingly maturing.However, still there is not a full understanding of the joint mechanisms between rock bolts and the surrounding rock and the dynamic response characteristics of the bolts under seismic loads.Therefore, establishing a reasonable dynamic analytical model and developing an efficient numerical calculation platform for the bolt are of great significance for analyzing the antiseismic stability of underground caverns.
As one of the main approaches for bolt support in underground caverns, fully grouted rock bolts have shown excellent performance in reinforcing the rock.Currently, researchers mainly focus on the force of the bolt.Recent research on basic mechanical transform mechanisms between the bolt and medium is mainly based on pull-out tests, including establishing theoretical models [1][2][3][4][5][6][7][8] and performing tests [9][10][11].However, the mechanical characteristics of the bolt in practical engineering projects are considerably different from those in observed pull-out tests.Freeman [12] and Wang et al. [13] first proposed the neutral point theory, which has been widely used for underground caverns.Many scholars have studied the action mechanism between the bolt and rock in tunnels based on this theory and obtained theoretical solutions [14][15][16][17][18] for the stress distribution of the bolt.However, these methods are mainly suitable for single cavern featuring a simple structure.Numerical methods should be studied for large-scale groups of underground caverns under complex geological conditions.Few studies have presented dynamic calculations for bolts to date.Xue et al. [19] analyzed the influence of seismic loading on bolts using FLAC and obtained the variation of the axial force of the bolt with dynamic loading times.Zhang and Xiao [20] proposed an algorithm that was suitable for antiseismic analysis to simulate bolts based on the numerical algorithm adopted by FLAC3D.However, these algorithms neglect the inhomogeneity of the stress distribution of the bolt and do not focus on damage evolution properties of the anchoring interface during earthquakes.
This paper proposes a new analytical algorithm for the interaction between the bolt and rock considering shear damage on the anchoring interface under a seismic load based on the bearing mechanism of fully grouted rock bolts in underground caverns.The rationality of the algorithm is verified; then, the dynamic response of an engineering case is evaluated.Several meaningful conclusions are obtained.

Analytical Model for the Interaction between the Bolt and Rock
Engineering practice shows that slip separation of the anchorage body is a common failure mode of bolts.Therefore, the following basic assumptions were made to study the interaction between the bolt and rock.(1) The bolt is made of a linear-elastic material, and only the axial deformation is considered.(2) The perfect combination between the bolt and mortar causes them to deform together.(3) The slip failure of the bolt occurs only on the contact interface between the anchorage body and rock.A one-dimensional local coordinate system was established along the axial direction of the bolt, whose positive direction points to the deep rock from the anchor head.A small piece of the anchorage body was taken for study, whose coordinate is  and length is d, as shown in Figure 1.The basic mechanical equations can be expressed as follows.
The balanced differential equation is The constitutive equation is The geometric equation is where () and () are the axial force and normal stress of the anchorage body at , respectively.() is the shear stress of the anchoring interface. a () and  a () are the axial strain and axial displacement of the anchorage body, respectively. a is the radius of the anchorage body. a is the composite elastic modulus of the anchorage body.
The controlled equation for the interaction between the bolt and rock can be obtained by combining (1)-(3): Supposing that the axial displacement of the rock at  is  r (), the shearing relative displacement between the rock and anchorage body is () =  r () −  a ().The trilinear shear slip model [21] was adopted to describe the relationship between () and (), as shown in Figure 2. The unified expression is  () =  () + . ( The values of  and  are determined according to three different stages. (I) Elastic Stage.The shear stress increases proportionately with the relative displacement, and there is no damage on the anchoring interface.
where  1 is the shear stiffness of the anchoring interface during the elastic stage, which is calculated by the following equation: where  r is the shear stiffness of the rock. r is 5-10 GPa/m for hard rock and 1.5-3 GPa/m for soft rock. m is the shear stiffness of the mortar and is determined by [22] where  m is the shear modulus of the mortar,  b is the radius of the bolt, and  is the thickness of the mortar.
(II) Plastic Softening Stage.The shear stress proceeds in the opposite direction when the relative displacement increases, and partial damage exists on the anchoring interface.
where  1 and  2 are the relative displacements corresponding to the peak shear strength and the end of the softening stage, respectively, and  2 is the shear stiffness of the anchoring interface during the softening stage. p is the peak shear strength of the anchoring interface, which meets the Mohr-Coulomb yield condition: where   and   are the cohesive force and internal frictional angle of the anchoring interface, respectively, and  r is the hydrostatic confining pressure of the rock perpendicular to the anchorage body.
(III) Slip Stage.The shear stress remains the same when the relative displacement increases, and the anchoring interface suffers extensive damage.
where   is the residual shear strength of the anchoring interface, which is written as Substituting ( 5) into (4) yields Equation ( 13) is the controlled equation for the interaction between the bolt and rock considering shear damage on the anchoring interface.
The damage coefficient, , is defined to describe the shear damage degree at different parts of the anchoring interface and is calculated by The value of  ranges from 0 to 1.  = 0, 0 <  < 1, and  = 1 correspond to states (I), (II), and (III), respectively.

Solution for the Stress of the Bolt Based on the Finite Difference Method
Equation ( 13) is a second-order linear differential equation with constant coefficients.It is difficult to obtain the functional expression of  a () using analytical methods.The finite difference method was adopted in this paper to acquire the numerical solution of  a ().Supposing that the bolt length is , the anchorage body is divided into  equal segments, and each segment length is Δ = /.The segmented points from  = 0 to  =  are numbered in order as 1, 2, . . .,  + 1.Then, the values of  r () and  a () at the th segmented point can be expressed as  r  and  a  , respectively.The variables   and   represent the material constants of the anchoring interface at the th segmented point.According to the finite difference method, (13) can be discretized as follows.
When  = 1, Equations ( 15)-( 17) can be rewritten in the matrix form as where Mathematical Problems in Engineering When [ r ] is known, (18) becomes a tridiagonal linear equations system for the unknown variable [ a ], which can be solved using the speedup method.
The values of (), (), and () at different segmented points are calculated by Based on the same deformation of the bolt and mortar, the normal stress,  b (), and axial force,  b (), of the bolt are obtained as

Dynamic Finite Element Calculation under Bolt Support
Regarding the dynamic finite element calculation for underground caverns under bolt support, the central difference method was mainly used to solve the differential equation of motion as follows [20]: where  is the lumped mass of the nodes of the finite element model, ä is the acceleration of the nodes,  ext and  int are the external and internal forces of the nodes, respectively,  damp is the damping force of the nodes, and  mg is the support reaction force offered by the bolts, which can be calculated by the following method.
Having learned the action results for the seismic load under bolt support at  −1 , the free deformation calculation of the rock without support was conducted first at   .This process can be described by solving the following differential equation of motion: The rock displacement at all segmented points of the bolt in the global coordinate system is obtained according to the interpolation theory of the finite element shape function.Then, the axial displacement of the rock at segmented points can be acquired by projecting the displacement in the vector direction of the bolt.Substituting [ r ] into (18) yields [ a ].
Having acquired [ a ] and [ r ], the functional relationship between () and () can be redefined and the matrices [], [], and [] can be renewed.
The bolts can restrict the rock deformation by providing a support reaction force.The shear stress of the anchoring interface can be equivalently transformed into the concentrated loads of the segmented points.Then, the concentrated loads are exerted on the rock elements in a contrary direction, which can simulate the anchoring effect of the bolts.Supposing that the shear stress of each segment of the bolt is linearly distributed, the equivalent load of each segmented point is expressed as where   is the equivalent load of the th segmented point:  = 1, 2, . . ., .
The three components of   in the global coordinate system are given as where , , and  are the directional cosines of the bolt.
According to the interpolation theory of the shape function, the support reaction force exerted on the rock by the th segmented point is deduced by where    is the support reaction force at the th node of the rock element in the , , and  directions,   is the value of the shape function at the th node of the element at the segmented point, and    represents the three components of   in the global coordinate system.
Iterating this process for all of the bolts yields the total support reaction force of the rock system.The rock displacement is updated by solving (24).Several iterations were repeated according to the steps above, and the convergent solution of the rock displacement can be obtained, which will help to obtain the internal force of the anchorage body.The detailed calculation procedure is shown in Figure 3.

Verification of the Numerical Model
The rationality of the analytical model for the bolt was verified by assessing two static calculation examples before the dynamic calculation was conducted.

Example 1.
During the pull-out test, the bolt is subjected to a tensile load, and the rock or concrete is approximately fixed; that is, [ r ] = 0.Then, (15) is modified as where  0 is the tensile load.Therefore, (18) can be rewritten as where the meanings of [], [ a ], and [] are the same as above and  1 =  1 /( a  a ) −  0 /( a 2  a Δ).
In this example, the bolt was buried directly in the concrete.The physical and mechanical parameters of the bolt and anchoring interface are provided in Table 1.Supposing that there is no damage on the interface of the bolt,   =  1 ,   = 0, and  = 1, 2, . . .,  + 1.The matrix [ a ] in (30) can be calculated directly without iterations.
The normal stress distributions of the bolt for  0 = 100 kN and  0 = 150 kN were calculated by using the bolt algorithm  in this paper and referring to the tests from Rong et al. [9], as shown in Figure 4.
The numerical values of the normal stress are nearly consistent with those of the experimental results, except for a slight difference in the quantity.The difference is smaller when  0 = 100 kN because there is only minor damage on the anchoring interface under smaller loads, and the interfacial shear stiffness was considered to remain constant during the experiment.For  0 = 150 kN, the difference is slightly larger because more damage occurs on the anchoring interface under larger loads.The interfacial shear stiffness was reduced during the experimental process but is thought to remain constant in the numerical simulation, which leads to  the smaller values in the simulation calculation than those of the experiment.

Example 2.
Example 2 was derived from a deeply buried rounded tunnel with an excavation radius of 5 m [16].The physical and mechanical parameters of the materials are provided in Table 2.A bolt at the top arch was monitored to observe its force state.The stress distributions of the bolt when the tunnel was excavated were obtained, as illustrated in Figure 5.
The shear stress is characterized by positive and negative values.Negative values indicate that the direction of the shear stress points to the free face of the cavern from 0 to 0.7 m, whereas positive values indicate that the direction of the shear stress points to the deep rock from 0.7 to 2.5 m.The location where the shear stress equals 0 is the location where the axial force is maximal, which corresponds with the neutral point theory.
Three monitoring points were set at 0.25, 1.25, and 2.25 m along the bolt length.The measured values and calculated values of the axial force are compared in Figure 5.The calculated values of the axial force derived from the three monitoring points show a closer match with the measured values, indicating the high accuracy of the bolt algorithm.
Therefore, the bolt algorithm proposed in this paper can reasonably reflect the factual bearing characteristics of the bolt during the pull-out test and in underground caverns.That is, the algorithm can satisfy the basic requirements for dynamic analysis.A three-dimensional finite element model consisting of 119,028 nodes and 112,905 elements of an 8-node hexahedron was built, as shown in Figure 6, which includes the diversion tunnel, main powerhouse, busbar chamber, main transformer room, tailrace tunnel, and surge shaft.The maximal characteristic mesh size was limited to less than 10 m, which    [23] is adopted to analyze the seismic response of underground caverns.The initial three-dimensional geostress field is captured by stress inversion of the measured points.The disturbing stress field after the excavation of the caverns is taken as the initial condition for the dynamic calculation.The bolts in the main powerhouse are arranged alternately according to @1.5 m × 1.5 m,  = 6 m/9 m.The bolts in the  main transformer room are arranged according to @1.5 m × 1.5 m,  = 6 m.The physical and mechanical parameters of the materials are provided in Table 3.

Engineering Case Study
The engineering project is in a frequently earthquakestricken area whose basic seismic intensity is VII.The seismic wave was imported from the bottom of the model.The  direction acceleration time-history is shown in Figure 8, and the  direction acceleration time-history is 2/3 that of the  direction.
The viscous-elastic artificial boundary was applied at the bottom and top of the model, and the free field artificial boundary was applied at the other sides.The partial damping was adopted as 0.157 for the rock.

Analysis of the Calculation Results
. Five monitoring bolts were arranged at five different positions of the mid-section of the 3 numbered unit sections to observe the force changes of the bolts during the earthquake, as shown in Figure 9.

Force Analysis of the
Bolts before the Earthquake.The stress distributions of the monitoring bolts after all of the excavation stages were completed are plotted in Figure 10.
The stress distributions conform to the neutral point theory, as shown in Figure 10.The maximal normal stress appeared at the neutral point, whereas the shear stress at both ends of the bolts reached the maximum value.The stresses of the bolts at the top arches were significantly smaller than those at the sidewalls, whereas the maximal normal stress and shear stress of the bolts at the sidewalls did not exceed the yield strength and peak shear strength, respectively.That is, the bolts were in good contact with the rock.

Force Time-History of the Bolts under Seismic Loads.
The time-history curves of the maximal normal stress and shear stress at the anchor head of the monitoring bolts during the earthquake are plotted in Figure 11.
The maximal normal stress grew consistently over time when under the seismic load.The maximal normal stresses of the bolts at the sidewalls in the main powerhouse changed considerably as their values increased by 102 MPa and 125 MPa after the earthquake.The seismic load had a considerable influence on the force of the bolts.The main interpretation of this result is that the rock maintained a long state of stress adjustment during the earthquake, and the unrecoverable plastic deformation accumulated continuously over time.This led to the continuous growing strain of the bolts, whereas the majority of the bolts were in the plastic zone of the rock.
The shear stress at the anchor head initially increased gradually over time.Then, damage on the anchoring interface at the sidewalls in the main powerhouse occurred when the shear stress reached the peak shear strength.As the damage accumulated over time, the shear stress decreased slowly.

Force Analysis of the
Bolts after the Earthquake.The stress distributions of the monitoring bolts when the seismic load was completed are plotted in Figure 12.
The stress distribution laws along the bolt length remained similar to those in Figure 10, except that the stress   values were larger.The neutral point positions of the bolts moved to the deep rock by 0-0.6 m, which shows that the influence scope of the earthquake expanded gradually to the anchoring depth.
normal bolts at sidewalls in the main powerhouse was larger, and the length of the sections with normal stress values exceeding 200 MPa reached 40-60% of the full length of the bolts.The damage length of the anchoring interface near the anchor head of the bolts at the sidewalls in the main powerhouse reached 0-0.3 m, which shows that the earthquake disturbed the high sidewalls of the main powerhouse more considerably.
Figure 13 illustrates the maximal normal stress distribution of the bolts in the main powerhouse after the earthquake.At the upstream sidewall near the diversion tunnel and downstream sidewall near the tailrace tunnel, the normal stress    values of a few bolts reached the yield strength (300 MPa), and the yield length reached 1.5-3.0m, whereas the normal stress of the bolts at other positions was smaller.
Figure 14 shows the damage coefficient distribution of the anchoring interface at the anchor head in the main powerhouse after the earthquake.Many of the anchoring interfaces in the main powerhouse suffered damage.At the upstream sidewall near the diversion tunnel and downstream sidewall near the tailrace tunnel, the damage coefficient of the anchoring interface of a few bolts reached 1.0, indicating that Mathematical Problems in Engineering the anchoring interface was in a slip state.These positions should be reinforced to guarantee the antiseismic safety of the bolt support system.Once slip occurs at the anchor head, the slip position will gradually expand to the anchoring depth under the subsequent seismic load, which will negatively impact the force of the bolts.

Conclusions
This paper proposed a dynamic algorithm for bolts using the trilinear shear slip model based on the basic mechanical equations of the anchorage body.The following conclusions can be drawn from this study: (1) It is feasible to simulate the force of the bolts in underground caverns under seismic loads using the bolt algorithm.The calculated stress distributions of the bolts corresponded well with the neutral point theory.
(2) It is effective to consider the shear damage on the anchoring interface because it can help reflect the damage evolution properties of the anchoring interface under seismic loads and can also provide new avenues for revealing the antiseismic failure mechanisms of the bolts in underground caverns.
(3) The normal stress of the bolts increases continuously during an earthquake, which may cause the bolts to yield.The shear stress of the anchoring interface at the anchor head increases rapidly; hence, this part is apt to be damaged or slip.
(4) The antiseismic safety evaluation of the bolt support system in underground caverns can be undertaken as a quantitative analysis of the bonding state of the anchoring interface and the material yield of the bolts.For this case study, a few bolts yielded, and the anchoring interface at the anchor head experienced slip at the upstream sidewall near the diversion tunnel and downstream sidewall near the tailrace tunnel in the main powerhouse.These positions should be the key zones for antiseismic support.
(5) This proposed algorithm mainly focused on the continuous contact state of the anchoring interface and ignored the discontinuity between the mortar and bolt, rock, or within mortar itself, which is helpful to simplify the dynamic finite element analysis.

Figure 1 :
Figure 1: Interaction between the bolt and rock.

Figure 3 :
Figure 3: Flowchart of the dynamic finite element calculation for the anchoring system at   .

Figure 4 :
Figure 4: Normal stress distributions of the bolt.

6. 1 .
Engineering Profile and Calculation Model.The Wudongde hydropower station is characterized by powerhouses located separately on the left and right banks of the Jinsha River.The underground powerhouse on the left bank includes 6 hydroelectric generating sets.The size of the main powerhouse is 310.0 m × 30.0 m × 84.2 m (length × width × height).
Shear stress distribution of the anchoring interface

Figure 5 :
Figure 5: Stress distributions of the bolt.

Figure 8 :
Figure 8: Time-history curve of the seismic acceleration.

Figure 9 :
Figure 9: Layout of the monitoring bolts.
Shear stress distributions of the anchoring interface

Figure 10 :
Figure 10: Stress distributions of the monitoring bolts (before the earthquake).

Figure 11 :
Figure 11: Stress time-history of the monitoring bolts.
Shear stress distributions of the anchoring interface

Figure 12 :
Figure 12: Stress distributions of the monitoring bolts (after the earthquake).

Figure 13 :
Figure 13: Maximal normal stress distribution of the bolts in the main powerhouse after the earthquake (unit: MPa).

Figure 14 :
Figure 14: Damage coefficient distribution of the anchoring interface at the anchor head in the main powerhouse after the earthquake.

Table 1 :
Physical and mechanical parameters.

Table 2 :
Physical and mechanical parameters.

Table 3 :
Physical and mechanical parameters.Note. b is the yield strength of the bolt.