Elastic-viscoplastic solution of a circular tunnel in Hoek-Brown rock mass

E-mail: linfei19891017@126.com *Corresponding author https://doi.org/10.1590/1679-78256328 ABSTRACT A closed-form analytical solution to rheological phenomena is obtained based on Hoek-Brown (HB) yield criterion and the generalized Bingham model. Once the circular tunnel is excavated, an initial plastic region is formed. When the total stresses fulfil the yield condition, the initial plastic region steps into the viscoplastic region, in which the displacement develops with time. The one-dimensional constitutive equation of the generalized Bingham model is transformed into a three-dimensional expression in polar coordinates by combining the plane strain condition. With the help of equilibrium differential equation, geometric equation and boundary condition, the stress and displacement in the initial plastic region can be expressed. At the moment of t , the stresses and displacements are described similar to the solutions of the elasto-plastic state. The stresses, strains and displacements after transformation can be expressed in the viscoplastic region. The proposed solution was verified by the closed-form solution of circular


List of Symbols
a Excavation radius σ r , σ θ Radial and tangential stresses ε r , ε θ Radial and tangential strains σ r vp , σ ϑ vp Radial and tangential viscoplastic stresses Derivative of variable to time

Introduction
Because of the urgent need of mineral resources, the mining dramatically deep into the earth, many coal and metal mines have been mined to a depth of more than 1 km. In order to facilitate the passage, deep buried tunnel construction has been common. However, deep buried roadway is located in the harsh environment of high confining pressure, high temperature and high permeability, which is easy to produce engineering geological disasters such as rheology, rock burst and large deformation of surrounding rock. Especially in the soft rock roadway, rock lithology metamorphism and surrounding rock failure occur frequently, the construction process is full of challenges [1]. This will seriously affect the stability of roadway surrounding rock, cause safety hazards to construction and traffic personnel, and cause huge economic losses to the roadway repeatedly renovated and maintained. Correctly predicting the deformation, stress and failure mode of surrounding rock can provide reliable reference for the stability evaluation of underground engineering.
The expansion of surrounding rock fracture range is time dependent. This leads to the nonlinear continuous deformation of surrounding rock, which makes it more difficult to maintain the stability of surrounding rock. For soft rock mass, the creep displacement is significant, and the deformation of soft rock mass can reach tens of millimeters every day, and when the support strength is weak, the deformation can reach tens of centimeters every day. In the 1960s, Lecian's concept of aging was first introduced into the analysis of stress and strain around tunnels [2]. An increasing number of underground engineering, due to high stress, the rheologic phenomena of the rock mass are more obvious, and research on the mechanical state of the rheological rock mass is of great value to engineering stability and safety design. Therefore, it is necessary to derive analytical solutions considering rheological factors.
Axisymmetric analysis of a circular tunnel is a classical problem. Many studies in existing research reports have been performed on the analytical solutions after excavation, and many engineering problems have been solved. In 1983, Brown et al. [3] proposed a method for analyzing the distribution of the stress and displacement in the plastic region into surrounding rock by using the Hoek-Brown (HB) criterion. Sharan, S. K. [4][5] presented stress and displacement field analysis and analytical solutions of the rock surrounding tunnels under the HB criterion. Carranza et al. [6][7]

Problem description
Although the rock mass is actually discontinuous medium, Fazelabdolabadi et al. [23] develops a Bayesian framework to quantify the absolute permeability of water in a porous structure from the geometry and clustering parameters of its underlying pore-throat network. Herein, A circular tunnel with a radius a is excavated in a homogeneous infinite isotropic rock mass subjected to a hydrostatic stress (p 0 ), the influence of porosity and permeability is not considered. The inner surface of the opening is subjected to a support pressure p i , and the inner support pressure p i equals the hydrostatic stress p 0 before the excavation. The elasto-viscoplastic circular tunnel is shown in Figure 1. The analytical model of the surrounding rock consists of three concentric annuli, i.e., the elastic region, the viscoplastic region, and the initial plastic region. To derive the analytical solution, the following assumptions are necessary: the rock mass is isotropic, uniform, and homogeneous; the plane strain condition is considered and the problem illustrated in Figure 1 can be considered to be axisymmetric. For plane strain problems, it can be assumed that the radial stress σ r is the minor principal stress, and the tangential stress σ θ is the major principal stress. Additionally, the strains are assumed to be small. The notation u represents the radial displacement, and the superscripts "c" and "e" represent plastic region and elastic region, respectively. The displacements and stresses are functions of the distance r to the center of the cross section and the time t when gravity is ignored.
The support pressure p i gradually decreases, and the deformation of the rock mass converges, when the roadway is excavated. When p i is lower than the elasto-plastic support pressure p c , the plastic region begins to appear outside the excavation face. After that, the material properties of the plastic region reach the residual value. Therefore, the macro brittle failure of the rock is the result of crack initiation, propagation, and connection in the medium. The elasto-brittleplastic model is shown in Figure 2, where AD represents the ideal plastic model, ACE represents the actual stress-strain curve of the rock, and ABE represents the assumed brittle-plastic model, the assumed brittle plastic model is used to analyze the stress and deformation of surrounding rock. For the problem of circular openings, the tangential stress (σ θ ) and radial stress (σ r ) are the major and minor stresses in circular tunnels, respectively. In this paper, two of the most commonly-used yield criteria are considered, the linear MC criterion and the nonlinear HB criterion where σ 3 is the minor principal stress, σ 1 is the major principal stress,  For axisymmetric problems, the stresses satisfy the differential equation of equilibrium in each region (without considering the body force of the surrounding rock), as shown in Eq. (7) / r r r r The geometric equation of this axisymmetric problem can be denoted as The difficulty of elastic-plastic problems is that the plastic deformation is nonlinear. By the theory of plastic, plastic strain is determined by the plastic potential, and rock mass materials generally obey the MC plastic potential function, which can be shown as Eq. (10) where ψ is the dilatancy angle of the rock mass, and The classical potential theory can be shown in the following equations, the radial and tangential plastic strains can be obtained: Combined the integration results of Eq. (11) with the elastic strain, and considering the influence of roadway excavation, the total strain of the rock mass can be expressed as is the shear modulus, µ is Poisson's ratio, and λ is the plastic strain multiplier. When where i p is the support force, − c R and + c R represent the inner and outer boundaries, respectively, at c r R = . When the radius tends to infinity, the radial stress equals the original in situ stress. Parameters

Surrounding rock model
In the quiescent state, the Bingham viscoplastic model can keep nonzero deviatoric stresses when the stress is less than the yield stress, and no plastic flow appears. Once the stress reaches the critical value, plastic flow occurs [24]. The generalized Bingham body is composed of Hookean, Newtonian, and St. Venant elements, so it is slightly different from the Bingham model. Elastic deformation occurs before applying stress is less than the yield stress, and the material manifests elasto-viscoplastic deformation when the yield stress is reached. Figure 3 shows the generalized Bingham model.
where dot "•" on the top of σ and ε denotes the derivative of time, E is the elastic modulus, s σ is the yield strength of the slider, and η is the viscosity coefficient of the dashpot.
With the aid of the correspondence law, convert E to 2G, convert η to 2η, convert ε to e ij , convert σ to s ij . It is generally considered that the volume deformation of rock mass is elastic, for the purpose of general description, the volume deformation is considered as elastic viscoplastic deformation as the deviator strain tensor. According to the above conversion method, one dimensional constitutive equation can be transformed into three-dimensional constitutive equation, so the three-dimensional equation can be denoted as where ij e is the deviatoric strain tensor; ij s is the deviatoric stress tensor; ij σ is the stress tensor, ij ε is the strain tensor, and both the stress and the strain tensors can be decomposed into the deviatoric stress and hydrostatic pressure, i.e., ; kk σ and kk ε represent the first invariants of the stress and strain, respectively; ij δ is the Kronecker symbol, when i j = , G is the shear modulus; K is the bulk modulus; η is the viscosity coefficient of the deviatoric stress and deviatoric strain; η′ is the viscosity coefficient of the volume strain ( kk ε ) and hydrostatic pressure ( kk σ ).
In a circular opening, Eq. (15) can be unfolded in a polar coordinate system as Eq. (16) where r z θ σ σ σ Θ = + + , which is the volume stress.
By substituting the plane strain condition ( 0 z z ε ε = =  ) into the third equation of Eq. (16), the axial stress can be reformulated as  (16), then using the equilibrium differential equation Eq. (7), can be replaced by Poisson's ratio µ , by simplification, the three-dimensional constitutive equation of this model in a polar coordinate system can be finally obtained as Eq. (19) , and c r σ and c θ σ are the radial stress and tangential stress, respectively, of the rock mass in the plastic region. Parameters r σ and θ σ indicate that rheological deformation occurs when the stress exceeds the plastic yield limit.

Elasto-plastic solutions under the HB yield criterion
The elastic-plastic solutions in the surrounding rock under the MC criterion have been studied, and the results will not be listed here. In the HB rock mass, the plastic region and the elastic region are formed, and the analytical solution of the stress and displacement can be expressed: By combining Eqs. (1), (5), (7) and (13) and C is an integration constant.
(2) Elastic region ( c R r ≤ < ∞ ) The stress, strain, and displacement of the elastic region can be easily obtained as Eqs. (22) and (23) where

Elasto-viscoplastic solutions of the HB surrounding rock
Because only the cross section of the tunnel is studied, sequential excavation is not considered. It is assumed that the long tunnel is excavated at one time here. In the context of the generalized Bingham calculation model, when the tunnel is excavated, elastic deformation occurs instantaneously. Because the sudden stress release on the excavation surface, the stress redistributed in the surrounding rock need to achieve a new equilibrium state, and a viscoplastic region with radius R t generates near the boundary of the tunnel.

Elastic solution of the rock mass at t=0
Because rheological test data indicate that elastic deformation is much larger than plastic deformation [25].
Consequently, the mechanical behavior of rock mass can be assumed to be elastic at 0 t = . Namely, once the tunnel is excavated, only elastic deformation occurs. Hence, when where the superscript "e" and the subscript "0" represent the solutions of the elastic state and the time Combining Eqs. (19), (29) and (30) leads to the differential equation as follows Integrating with respect to r and t, Eq. (31) can be solved as Eq. (32) where ( ) f t is the function of the time t, and ( ) g r is the function of the radius r, yet to be determined.
By means of the first and fourth equations of Eq. (13), and utilizing Eq. (30), the radial and tangential stress can be expressed as It should be noted that the equations of Eq. (33) are valid within the initial plastic region, and it was assumed that the surrounding rock at radius r is always in a plastic state.
The unknown function 1 ( ) g r can be determined by the fourth equation in the initial condition of Eq. (13), namely, Because the condition of displacement continuity at the position of In fact, the deformation of rock mass is the process that the stress is adjusted continuously to reach balance after the tunnel excavation. In the viscoplastic region, the stress and deformation of the rock mass changes with time. When the time t → ∞ , the solutions of this paper and those gotten by the elasto-plastic method are consistent at the corresponding position, and at this time, even if the time increases, the viscoplastic region will not change again.
Eventually, the initial plastic region and the viscoplastic region becomes the plastic region, and the solution of the calculated model corresponds to the foregoing section 4.1.

Calculation methodology
By inputting geometric and mechanical parameters of surrounding rock, ( ) f t and ( ) f t

Example verification
In order to analyze the influence of dilatancy angle on stress and displacement in the plastic region, dilation angles of 0º and 30º are selected. The parameters of the HB criterion are converted into equivalent MC parameters. By doing so, the solutions of the two yield criteria can be compared under similar geological conditions. The specific conversion method [26] is as follows: by equating the areas covered by the two envelopes of the HB and MC criteria, Eq. (43) and Eq. (44) are formed, which may be solved for the equivalent MC rock mass parameters when the HB constants are known [27]. This results in the following equations for the angle of friction φ and cohesive strength c: where, γ is the bulk density of the rock mass and H is the buried depth of the rock mass.
To validate the analytical solutions, the typical input data used by Carranza-Torres [7] were used. After equivalent conversion, the mechanical parameters of the HB and MC criteria are listed in Table 1.

The solutions of functions f(t) and ∫ f(t)dt
The f(t) function in implicit equation (41) is solved by MATLAB software. Figure 5  to zero. ∫ f(t)dt always increase with time. In MC rock mass, when the support force increases from 0MPa to 0.75MPa, the peak value of f(t) decreases from 2.984 to 2.243, and the decrease percentage is 24.83%. However, when the supporting force increases from 0 MPa to 0.75 MPa in HB rock mass, the peak value of f(t) decreased from 3.458 to 2.286, the percentage of decrease was 33.89%. It can be seen that the f(t) function of HB rock mass is more affected by the support force. In any rock mass, the value of ∫ f(t)dt decreases with the increase of support force. In MC and HB rock mass, when the support force is 0 MPa, the peak value of f(t) increases from 2.984 and 3.458 to 4.552 and 4.899, respectively, with the increase percentage of 52.55% and 41.67%. It can be seen that the f(t) function of MC rock mass is more affected by the dilatancy angle.  smaller and smaller. When t=60 days, the circumferential stress in viscoplastic region and initial plastic region is continuous. Fig. 6 The changes in the radial and tangential stresses with time and radius under the HB criterion.

Analytical solutions of stress and displacement
The dilation angle is an important factor affecting the displacement. Figure 7 shows the variation in the radial displacement of the tunnel with time and dimensionless r/a under the MC criterion and HB criterion when the dilation angle is 0º and 30º. Obviously, the radial displacement increases gradually with increasing time. Meanwhile, the displacement changes rapidly at the value of r/a from 1 to 3 and tends to be flat after r/a=3 when ψ=0º, which indicates that large and rapid deformation occurs near the excavation face of the tunnel. For the condition of ψ=30º, the displacement varies rapidly when r/a increases from 1 to 2. Compared with the displacement of ψ=0º, the displacement of the excavation face increases 3.44 times when t=60 days. The radial displacement increases with increasing dilatancy angle. The radial displacement obtained in HB rock mass is larger than that in MC rock mass.

Displacement Verification of the HB yield Criterion
To prove the correctness of the results obtained in this paper, the verification results of the HB rock mass are visualized in Figure 8. After calculation, the state at t = 300 days is selected. Because the viscoplastic region does not change at this time, the final state becomes an elastic-plastic rock mass. Some examples for the HB rock mass are given (Brown [3], Sharan [4], Zhang [8], Lee and Pietruszczak [28], and Park [29]

Influence of support force
The dimensionless final plastic radii (R c /a) with different supporting forces (p i ) are calculated and compared with the analytical results of Zhang [8], as shown in Table 2. As seen from the table, the radii of the plastic region calculated in this paper are consistent with Zhang [8]. With the increase in the supporting force, the dimensionless ρ 0 /a and R c /a both decrease. This law shows that applying the support force can effectively reduce the extension of the plastic region into the surrounding rock of the tunnel. Figure 9 shows the initial and final plastic radii with different supporting forces (p i ) for the HB and MC rock masses. The relationship between the radius of the initial plastic region and the supporting force is linear, and the R-square is as high as 99.9%. Under the two yield criteria, the slope of the straight line is the same. However, the radii of the initial plastic region under the HB yield criterion are larger under the same supporting force. The radii of the final plastic region also decrease with increasing support force. The fitting curve in Figure 9 (b) shows an approximate linear relationship under the MC yield criterion; the R-square is 99.4%. Under the HB yield criterion, the relationship between the support force and radii of the final plastic region is nonlinear. When the support force is less than 0.4MPa, the radii of the final plastic region under the HB yield criterion are larger than that under the MC yield criterion. When the support force is greater than 0.4MPa, the opposite is true.

Influence of Young's modulus
In order to analyze the influence of Young's modulus on the radius of viscoplastic region and the radial stress at the outer radius of initial plastic region, Young's modulus is taken as 2, 5, 10, 20 and 40GPa respectively, the creep time is taken as 4, 10 and 60 days respectively. As shown in Figure 10 (a), in MC rock mass, the radius of viscoplastic region increases with the increase of Young's modulus. Under three different time states, the radius of viscoplastic region increases from 1.573m, 1.677m and 2.095m to 2.178m, 2.246m and 2.262m respectively with the increase of Young's modulus from 2GPa to 40GPa, and the percentages of increase are 38.46%, 33.93% and 7.97% respectively. Similarly, as shown in Figure 10 (b), in HB rock mass, the radius of viscoplastic region increases with the increase of Young's modulus. Under three different time states, the radius of viscoplastic region increases from 1.455m, 1.630m and 2.170m to 2.219m, 2.381m and 2.417m respectively with the increase of Young's modulus from 2GPa to 40GPa, and the percentages of increase are 52.51%, 46.07% and 11.38% respectively. The location of the elastic-plastic interface, that is, the radius of the viscoplastic region, increases with the increase of time and Young's modulus. The larger the young's modulus is, the faster the viscoplastic region develops and the earlier the rock mass reaches equilibrium. The radius of the viscoplastic region in HB rock mass under the same Young's modulus, the viscoplastic radius of HB rock mass develops faster than that of MC rock mass. In MC rock mass, when Young's modulus increases from 2GPa to 40GPa, the radial stress at the outer radius of the initial plastic region decreases with time. Under three different time states, the radial stress decreases from 3.704MPa, 3.131MPa and 1.539MPa to 1.312MPa, 1.412MPa and 1.104MPa respectively, and the percentage of decrease was 64.58%, 54.90% and 28.27%, respectively, as shown in Figure 11 (a). But in the HB rock mass, the radial stress decreases from 4.780MPa, 3.602MPa and 1.356MPa to 1.223MPa, 0.851MPa and 0.780MPa respectively, and the percentage of decrease was 74.41%, 76.37% and 42.48%, respectively, as shown in Figure 11 (b). Through comparative analyses, the radial stress at the outer radius of the initial plastic region in MC and HB rock mass decreases with the increase of time and Young's modulus, but the decrease percentage of HB rock mass is larger than that of MC rock mass. Fig. 11 Effect of Young's modulus on radial stress at the outer radius of initial plastic region for MC (a) and HB (b) rock mass

Conclusions
The paper presented elasto-viscoplastic solutions for the stress and displacement of a circular roadway in an elastobrittle-plastic rock mass with the generalized Bingham model and HB failure criterion. Rock rheology is the result of the continuous accumulation of internal damage, which eventually leads to rock burst. The analytical solution is time dependent, it can obtain the stress and displacement at any time in the process of surrounding rock deformation, this is the biggest difference from the static solution. The HB yield criterion can be used to predict the failure of surrounding rock to a greater extent and provide more reliable guidance for engineering design.
According to the mechanical parameters of HB rock mass, those of MC are deduced. The analytical results obtained under the two yield criteria are compared. When 0 t r R ρ ≤ ≤ , the rock mass exhibits viscoplastic properties, and the radius increases with time. When R t =R c , the rock mass behaves brittle-plastic. The unknown function f(t) and its integral (cf. Figure 5) can be calculated for an arbitrary time form through Eq. (41).
The radial and tangential stresses in the plastic region at the moment which t >0 are always smaller than at t=0 under the HB criterion. Comparing the results of this study with previous research results in the literature, the displacement under the HB rock mass shows very good agreement with the solution of Lee & Pietruszczk, Zhang and Park, which indicates the correctness of the new approach. In this paper, the influence of the dilation angle and the support force is studied. When the dilation angle becomes 30º, the displacement increases noticeably, and the displacement under the HB yield criterion is larger than that under the MC yield criterion at the same time. With the increase in the supporting force, the dimensionless ρ 0 /a and R c /a both decrease. The radius of viscoplastic region