Time-Dependent Behavior of a Circular Symmetrical Tunnel Supported with Rockbolts

: Under the effect of initial stress and excavation disturbance, there exists interaction between rock mass and rockbolt in deeply buried tunnels. In order to fully explore the mechanism of rock mass supported with rockbolts, this article studied the time-dependent behavior of the rock mass supported with discretely mechanically or frictionally coupled (DMFC) rockbolts. The interaction model elastic solutions under distributed force model were analyzed, then the viscoelastic analytical solutions were conducted to describe the rheological properties of the coupling model, and the solutions were acquired by setting the constitutive models of the rockbolt and rock mass in terms of a one-dimensional Kelvin model and a three-dimensional Burgers model based on material properties and dimension. Several examples were performed and the inﬂuence of initial stress σ 0 , the viscosity parameters η 1 and η 2 of the three-dimensional Burgers model as well as the pre-tension T 0 on reinforcement effect were analyzed. According to the proposed model, the smaller η 2 is or the larger the pre-tension T 0 is, the more effective the support effect. However, when the pre-tension is too large, the support effect is no longer signiﬁcantly enhanced. In addition, the early reinforcement effect is controlled by the ﬁrst creep stage in the Burgers model while the ultimate support effect is mainly inﬂuenced by the viscosity coefﬁcient of the second creep stage in the Burgers model. This research can provide an important theoretical reference to guide the parameter design of rockbolt reinforcement engineering in a circular symmetrical tunnel.


Introduction
Rockbolts have been extensively employed for rock reinforcement in civil and mining engineering for many years [1]. However, the interaction mechanism between rockbolt and rockmass is still not clear at present [2]. Some engineering materials have certain rheological characteristics [3]. Due to the time-dependent behavior of the material, the rock mass displacement and stress may increase after excavation, and may last for several months or years [4]. Understanding the rheological mechanism of rock mass supported with rockbolts is essential for engineering design optimization.
Various approaches have been established to describe the time-dependent behavior of rock mass or rockbolt based on analytical solutions [5][6][7], empirical approaches [8][9][10], and numerical Symmetry 2018, 10, 381 2 of 16 methods [11][12][13]. Li et al. [14] studied the effect of the stress, creep coefficient, and geometry parameters on the stress relaxation of rockbolts. However, the time-dependence of rock mass was ignored. Wang et al. [15] presented the viscoelastic solutions of the coupling model with a distributed force model through using several common rheological models, since the selected rheological model does not consider material properties and dimensional effects. Therefore, it cannot fully reflect the true rheological behavior of the coupling model. Sharifzadeh et al. [16] simulated the long-term deformation of a tunnel excavated in a weak rock mass, and noted that the instability of the tunnel may be attributed to the stress induced by the creep behavior of the surrounding rock. Above-research have primarily been focused on the evolution of rock mass or rockbolt, few studies have been conducted on the coupling rheological behavior of the interaction model [15].
The research work in this paper is to study the coupling rheological behavior of the rock mass supported with discretely mechanically or frictionally coupled (DMFC) rockbolts based on reasonable rheological model and dimension. The elastic and viscoelasticity solutions were first discussed, and the parametric study was subsequently conducted. Finally, the long-term tunnel stability was evaluated.

Principle of the Coupling Model
The rockmass would deform under the initial stress and support force of the rockbolt. Here, we assumed that the rockbolt-rock mass will not be decoupled, so the deformation of the rockbolt is coordinated with the rock mass, and the coupling process is depicted in Figure 1. methods [11][12][13]. Li et al. [14] studied the effect of the stress, creep coefficient, and geometry parameters on the stress relaxation of rockbolts. However, the time-dependence of rock mass was ignored. Wang et al. [15] presented the viscoelastic solutions of the coupling model with a distributed force model through using several common rheological models, since the selected rheological model does not consider material properties and dimensional effects. Therefore, it cannot fully reflect the true rheological behavior of the coupling model. Sharifzadeh et al. [16] simulated the long-term deformation of a tunnel excavated in a weak rock mass, and noted that the instability of the tunnel may be attributed to the stress induced by the creep behavior of the surrounding rock. Above-research have primarily been focused on the evolution of rock mass or rockbolt, few studies have been conducted on the coupling rheological behavior of the interaction model [15]. The research work in this paper is to study the coupling rheological behavior of the rock mass supported with discretely mechanically or frictionally coupled (DMFC) rockbolts based on reasonable rheological model and dimension. The elastic and viscoelasticity solutions were first discussed, and the parametric study was subsequently conducted. Finally, the long-term tunnel stability was evaluated.

Principle of the Coupling Model
The rockmass would deform under the initial stress and support force of the rockbolt. Here, we assumed that the rockbolt-rock mass will not be decoupled, so the deformation of the rockbolt is coordinated with the rock mass, and the coupling process is depicted in Figure 1. The reinforcement mechanism for DMFC rockbolts can be divided into three models: the distribution force model, the point load model, and the equivalent material model. Distribution force models regard the force of the rockbolt to rock mass as the tangential uniform distribution of the surface force. The point load model simplifies a circular tunnel as a spring device, which reduces the force to a pair of opposing forces of equal magnitude. The equivalent material model divides the rock mass into a reinforced and an original zone. The rock mass and rockbolt in the reinforced area are regarded as homogeneous media. However, the equivalent material model neglects the mechanism of the action between the rockbolt and the rock mass in the reinforced zone; thus it cannot reflect the true rheological state of the rockbolt and rock mass. This paper employs the distribution force calculation considering the analytic complexity of the point load model and the hypothesis error of the equivalent material model. The distributions of the three models are shown in Figure 2. The reinforcement mechanism for DMFC rockbolts can be divided into three models: the distribution force model, the point load model, and the equivalent material model. Distribution force models regard the force of the rockbolt to rock mass as the tangential uniform distribution of the surface force. The point load model simplifies a circular tunnel as a spring device, which reduces the force to a pair of opposing forces of equal magnitude. The equivalent material model divides the rock mass into a reinforced and an original zone. The rock mass and rockbolt in the reinforced area are regarded as homogeneous media. However, the equivalent material model neglects the mechanism of the action between the rockbolt and the rock mass in the reinforced zone; thus it cannot reflect the true rheological state of the rockbolt and rock mass. This paper employs the distribution force calculation considering the analytic complexity of the point load model and the hypothesis error of the equivalent material model. The distributions of the three models are shown in Figure 2.

Elastic Solutions
The relative distributions of the rockbolts and rock mass is shown in Figure 3. The following assumptions are taken: (a) the tunnel is circular and deep; (b) the support force in tunnel opening P0 and the reinforcement force Pρ are uniformly distributed on the excavation surface and the interface between the reinforced zone and the original zone; (c) the problem is axisymmetric, and the lateral pressure coefficient Ka = 1; and (d) the deformation is small.
For the axisymmetric problem under polar coordinates, the general stress formula and displacement formulae can be expressed as:

Elastic Solutions
The relative distributions of the rockbolts and rock mass is shown in Figure 3. The following assumptions are taken: (a) the tunnel is circular and deep; (b) the support force in tunnel opening P 0 and the reinforcement force P ρ are uniformly distributed on the excavation surface and the interface between the reinforced zone and the original zone; (c) the problem is axisymmetric, and the lateral pressure coefficient K a = 1; and (d) the deformation is small.
For the axisymmetric problem under polar coordinates, the general stress formula and displacement formulae can be expressed as: where σ ρ , σ θ and u ρ are respectively the radial stress, tangential stress and displacement; E is Young's modulus; µ is Poisson's ratio; ρ is the radial coordinate; A, B, and C are coefficients. where σρ, σθ and uρ are respectively the radial stress, tangential stress and displacement; E is Young's modulus; μ is Poisson's ratio; ρ is the radial coordinate; A, B, and C are coefficients. When the tunnel is not excavated, the displacement elastic solution of the rock mass under the initial stress is as follows: In the reinforced zone, the mechanical problem can be regarded in terms of the elastic mechanics of the rear wall cylinder under internal and external pressure: the internal force arises from the support of the rockbolt. A mechanical analysis diagram is shown in Figure 4b. The general forms expressions for the stress and displacement in the reinforcement zone (r < ρ < R) can be written as: where σρ1, σθ1 and uρ1 are respectively the radial stress, tangential stress, and displacement in the reinforced zone and A1 and C1 are the pending coefficients.
In the original zone, the mechanical problem can be regarded in terms of the elastic mechanics of an infinite region affected by the initial stress σ0 and pressure on the inside of a circular hole. A mechanical analysis diagram is shown in Figure 4a. The general forms expressions of the stress and displacement in the original zone (ρ > R) are as follows: When the tunnel is not excavated, the displacement elastic solution of the rock mass under the initial stress is as follows: In the reinforced zone, the mechanical problem can be regarded in terms of the elastic mechanics of the rear wall cylinder under internal and external pressure: the internal force arises from the support of the rockbolt. A mechanical analysis diagram is shown in Figure 4b. The general forms expressions for the stress and displacement in the reinforcement zone (r < ρ < R) can be written as: where σ ρ1 , σ θ 1 and u ρ1 are respectively the radial stress, tangential stress, and displacement in the reinforced zone and A 1 and C 1 are the pending coefficients.
In the original zone, the mechanical problem can be regarded in terms of the elastic mechanics of an infinite region affected by the initial stress σ 0 and pressure on the inside of a circular hole. A mechanical analysis diagram is shown in Figure 4a. The general forms expressions of the stress and displacement in the original zone (ρ > R) are as follows: where σ ρ2 , σ θ 2 and u ρ2 are respectively the radial stress, tangential stress, and displacement in the original zone, and A and C are coefficients.
Symmetry 2018, 10, x FOR PEER REVIEW 5 of 17 where σρ2, σθ2 and uρ2 are respectively the radial stress, tangential stress, and displacement in the original zone, and ' A and ' C are coefficients. The stress boundary conditions of the pre-tension rockbolts at the tunnel face can be expressed as follows: where P0 is the support force in tunnel opening; Pρ is the reinforcement force; L is the length of the free part of the rockbolt; L  is the deformation of the rockbolt in the axial direction, x  is the pre-tension length of the rockbolt; Sz is the rockbolt spacing in the longitudinal direction; Sθ is the rockbolt spacing in the tangential direction; The other boundary conditions can be expressed as follows: The stress boundary conditions of the pre-tension rockbolts at the tunnel face can be expressed as follows: where P 0 is the support force in tunnel opening; P ρ is the reinforcement force; L is the length of the free part of the rockbolt; L is the deformation of the rockbolt in the axial direction, x is the pre-tension length of the rockbolt; S z is the rockbolt spacing in the longitudinal direction; S θ is the rockbolt spacing in the tangential direction; u 0 ρ1(ρ=R) and u 0 ρ1(ρ=r) are the displacements of the original state without support; A b is the area of the rockbolt cross section; and E b is the elastic modulus of the rockbolt.
The other boundary conditions can be expressed as follows: The elastic solutions of the rockbolt and rock mass can be presented by combining the displacement and stress conditions. Solutions have been discussed [17][18][19][20].

Viscoelastic Solutions
Wang et al. [15] selected several common rheological models to explain the rheological characteristic of rockbolts and rock masses; the Maxwell and the Kelvin model are respectively used to describe the rheological properties of the rockbolts and rock masses in this paper. While the Maxwell model cannot characterize the elastic after-effect and deformation limit characteristics of rock mass, the Kelvin model cannot characterize the transient and viscous flow of a rock mass. However, in practical engineering, the rock mass will have transient deformation and viscous flow due to the disturbance of excavation [3]. Therefore, there may be some drawbacks in fully expressing the rheology of a rock mass based on the selected models. Choosing a rheological model that is more consistent with the material properties to reflect the time-dependency of the rockbolt-rock mass coupling model fully is necessary.

Rheological Model of Rockbolt and Definition of Operator Function
With respect to the material properties, the rockbolt is primarily composed of a hot rolled ribbed bar, the amount of creep for steel is small, and the creep of the steel is dominated by stable creep. With respect to the scale effect, the axial stiffness of the rockbolt is much greater than the tangential stiffness, and the rockbolt can be regarded as an ideal one-dimensional viscoelastic material. Based on the analysis of the rheological properties of the steel bar, this paper uses a one-dimensional Kelvin model to describe the rheological properties of the rockbolt.
The variation between the viscoelasticity and the one-dimensional elastic parameters of the rockbolt is shown in Equation (8), where Q b (D) and P b (D) are the operator functions of the one-dimensional Kelvin rheological model. The space parameter transformation based on the Laplace transform is The Laplace form of the one-dimensional Kelvin model operator is where η b is the rockbolt's viscosity coefficient, s is the variable obtained from the Laplace transformation.

Rheological Model of Rock Mass and Definition of Operator Function
With respect to the material properties, the characteristics of the rock mass are complicated and changeable in underground engineering. The Burgers model consists of the Maxwell model and the Kelvin model; the model can describe rheological properties such as elastic deformation, creep, stress relaxation, elastic effect and viscous flow. With respect to the scale effect, the size of one direction is much longer than the other two directions for some engineering fields, the problem should be considered as three-dimensional problem [21]. In tunnel engineering, the excavation direction is substantially longer than the other two directions; thus, the rock mass should be regarded as a three-dimensional viscoelastic material. The one-dimensional constitutive equation needs to be extended to a three-dimensional form, and the three-dimensional viscoelastic constitutive relation can be derived from the corresponding relation of the one-dimensional viscoelastic constitutive relation. The distribution of rheological model diagram is shown in Figure 5. For the rock mass, the three-dimensional viscoelastic constitutive equation can be obtained by developing the one-dimensional constitutive equation based on the generalized Hooke's law The relationships among the shear modulus, bulk modulus, Young's modulus, and Poisson's ratio μ are Therefore, the constitutive models of the elastic and viscoelastic materials can be written as where Sij is the partial stress tensor; eij is the partial strain tensor; σij is the stress tensor; ij ε is the strain tensor; and '( ) For the rock mass, the three-dimensional viscoelastic constitutive equation can be obtained by developing the one-dimensional constitutive equation based on the generalized Hooke's law The relationships among the shear modulus, bulk modulus, Young's modulus, and Poisson's ratio µ are Therefore, the constitutive models of the elastic and viscoelastic materials can be written as where S ij is the partial stress tensor; e ij is the partial strain tensor; σ ij is the stress tensor; ε ij is the strain where s is the Laplace variable, and the operator functions of the Laplace space form in three-dimensional Burgers model can be expressed as where η 1 is the viscous coefficient of the first creep stage; η 2 is the viscous coefficient of the second creep stage; G 0 is the elastic shear modulus of the rock mass; G 1 is the corresponding viscoelastic shear modulus of the rock mass; and K is the elastic bulk modulus of the rock mass. The Laplace space solutions can be obtained by replacing the elastic parameters of the elastic solutions with the viscoelastic parameters and substituting the Laplace operator function of the selected rheological model, which enables the viscoelastic problem to be solved. The viscoelastic solution can be obtained by the Laplace inverse.
The parameters of the elastic solutions are replaced by viscoelastic parameters and substituted into the operator functions of the corresponding rheological model. The main process is expressed as follows

Analytical Solutions in the Laplace Space
The Laplace space analytical solutions can be described by substituting in the operator function of the corresponding rheological model, and the stress and displacement fields of the rock mass are expressed in Appendix A (A1)-(A6).
Then, the tunnel surface support stress can be calculated as follows Therefore, the axial force of the rockbolt can be written as where the tensioned rockbolt deformation is . Hence, the axial force of rockbolt can be expressed in Appendix A (A7).

Axial Force Changes over Time
The corresponding parameter values are provided in Table 1. The analytical curves can be obtained by substituting the parameters into the corresponding expression and performing the inverse Laplace transformation.  As shown in Figure 6a, the influences of η 1 and σ 0 on the axial force evolve over time. The results indicate that η 1 affects not only the axial force initial value but also the rockbolt rheological state. When η 1 = 3.0 × 10 10 and 5.0 × 10 10 Pa·s, the axial force decreases to a stable value over time, the rockbolt rheological state is stress relaxation. When η 1 = 7.0 × 10 10 Pa·s, the axial force increases to a stable value over time, the rheological state of the rockbolt is stable creep. The axial force will converge to a stable value regardless of whether the rheological state of a rockbolt is creep or stress relaxation; thus η 1 does not affect the stable value of the axial force and only influences the early support effect. When η 1 and η 2 remain unchanged, a larger σ 0 will generate a larger initial force and higher stability. Figure 6b indicates that the influences of η 2 and σ 0 on the axial force evolve over time.
The results indicate that η 2 affects not only the initial value and stable value of axial force but also the rheological state of the rockbolt; thus, η 2 affects not only the early reinforced effect but also the final support effect. When η 2 = 1.0 × 10 11 Pa·s, 1.5 × 10 11 Pa·s, and 2.0 × 10 11 Pa·s, the axial force decreases to a stable value over time; thus, the rheological state of the rockbolt is stress relaxation. A smaller value of η 2 produces a larger initial value and convergent value of the axial force. When η 2 = 3.0 × 10 11 Pa·s, the axial force increases to a stable value over time, and the convergent value is less than η 2 = 1.0 × 10 11 Pa·s, 1.5 × 10 11 Pa·s, and 2.0 × 10 11 Pa·s; thus, the rheological state of the rockbolt is creep. As shown in Figure 6c, the influence of η b on axial force evolves over time. When η b = 1.0 × 10 20 Pa·s and 3.0 × 10 20 Pa·s, the axial force decreases to a stable value over time, the rockbolt rheological state is stress relaxation. When η b = 5.0 × 10 20 Pa·s, the axial force is gradually reduced to zero and then increases to a stable value over time; the stable value is the same as when η b = 1.0 × 10 20 Pa·s and 3.0 × 10 20 Pa·s. When η b = 6.0 × 10 20 Pa·s, the axial force of the rockbolt increases to a stable value over time; thus, the rockbolt rheological state is creep, and the initial value is the same as when η 1 = 3.0 × 10 10 Pa·s. Figure 6d reveals that the influence of T 0 on the axial force evolves over time; larger values of T 0 produce a greater absolute value of initial axial force. Thus, the larger is T 0 , the better the support effect. When T 0 = 4.0 × 10 4 N, 5.0 × 10 4 N, and 5.5 × 10 4 N, the axial force absolute value gradually decreases over time and converges to a fixed value. A larger value of T 0 yields a higher fix value. But, when the T 0 is too large, the fix value will not increase significantly. At this time, the rockbolt support to the rock mass is in the stage of relaxation.
Pa·s. Figure 6d reveals that the influence of T0 on the axial force evolves over time; larger values of T0 produce a greater absolute value of initial axial force. Thus, the larger is T0, the better the support effect. When T0 = 4.0 × 10 4 N, 5.0 × 10 4 N, and 5.5 × 10 4 N, the axial force absolute value gradually decreases over time and converges to a fixed value. A larger value of T0 yields a higher fix value. But, when the T0 is too large, the fix value will not increase significantly. At this time, the rockbolt support to the rock mass is in the stage of relaxation.

Radial Stress Changes over Time
As shown in Figure 7, the case considered in this test primarily investigated how the viscosity coefficients η1, η2 and σ0 influenced the radial stress. Figure 7a shows the jump value fitting curve between the reinforced and the original zone when the time is zero. When η2 is fixed, the fitting curve is When η1 is fixed, the fitting curve is As shown in Figure 7a, η1 and η2 are negatively correlated with the radial stress jump value, and the slope of y1 is greater than the slope of y2, which indicates that the effect of η1 on the radial stress jump value is larger than the effect of η2. The radial stress jumps at the junction of the reinforced zone and the original zone due to the support of the rockbolt. Figure 7b shows the radial stress fitting curve in the excavation surface when the time is zero; when η2 is fixed, the fitting curve is 3 When η1 is fixed, the fitting curve is Note that η1 and η2 are negatively correlated with the radial

Radial Stress Changes over Time
As shown in Figure 7, the case considered in this test primarily investigated how the viscosity coefficients η 1 , η 2 and σ 0 influenced the radial stress. Figure 7a shows the jump value fitting curve between the reinforced and the original zone when the time is zero. When η 2 is fixed, the fitting curve is y 1 = −4334.5 × x + 49929. R 2 = 0.9198. When η 1 is fixed, the fitting curve is y 2 = −1160.3 × x + 27924. R 2 = 0.9513. As shown in Figure 7a, η 1 and η 2 are negatively correlated with the radial stress jump value, and the slope of y 1 is greater than the slope of y 2 , which indicates that the effect of η 1 on the radial stress jump value is larger than the effect of η 2 . The radial stress jumps at the junction of the reinforced zone and the original zone due to the support of the rockbolt. Figure 7b shows the radial stress fitting curve in the excavation surface when the time is zero; when η 2 is fixed, the fitting curve is y 3 = −19850 × x + 214942. R 2 = 0.936. When η 1 is fixed, the fitting curve is y 4 = −3148.1 × x + 109336. R 2 = 0.9678. Note that η 1 and η 2 are negatively correlated with the radial stress in excavation surface. The slope of y 3 is greater than y 4 , which indicates that the effect of η 1 on the excavation surface radial stress is also larger than the effect of η 2 . Therefore, the magnitude of the radial stress in the excavation surface of the tunnel is equivalent to the supporting force of the rock mass in the coupling model, and the greater the support force is, the better the support effect. The smaller η 1 or η 2 is, the better the support effect in the excavation moment. Figure 7c shows the evolution of the radial stress of the reinforced zone monitoring point ρ = 4 m over time (influence of η 1 and σ 0 ). The absolute value of radial stress gradually decreases over time because the stress of the rock mass is gradually released over time. When the parameters η 1 and η 2 are constant, a larger σ 0 yields a greater absolute value of the initial and convergence values. When the parameter η 2 is fixed, a larger η 1 produces a smaller absolute value of the initial radial stress. The radial stress converges to a fixed value over time, and the fixed value is independent of η 1 . A smaller η 1 produces a larger change rate. Figure 7d shows the evolution of the radial stress over time for the reinforced zone monitoring point ρ = 4 m (influence of η 2 and σ 0 ). The absolute value of radial stress gradually decreases over time; when the parameters η 1 and η 2 are fixed, a greater value of σ 0 yields a higher absolute value of initial radial stress and convergence value. When η 1 is constant, a smaller η 2 produces a larger absolute value of the initial and convergence values. Thus, the smaller is the η 2 , the better is the support effect. η 2 does not have a significant effect on the rate of change, which remains fairly constant. the excavation surface radial stress is also larger than the effect of η2. Therefore, the magnitude of the radial stress in the excavation surface of the tunnel is equivalent to the supporting force of the rock mass in the coupling model, and the greater the support force is, the better the support effect. The smaller η1 or η2 is, the better the support effect in the excavation moment. Figure 7c shows the evolution of the radial stress of the reinforced zone monitoring point ρ = 4 m over time (influence of η1 and σ0). The absolute value of radial stress gradually decreases over time because the stress of the rock mass is gradually released over time. When the parameters η1 and η2 are constant, a larger σ0 yields a greater absolute value of the initial and convergence values. When the parameter η2 is fixed, a larger η1 produces a smaller absolute value of the initial radial stress. The radial stress converges to a fixed value over time, and the fixed value is independent of η1. A smaller η1 produces a larger change rate. Figure 7d shows the evolution of the radial stress over time for the reinforced zone monitoring point ρ = 4 m (influence of η2 and σ0). The absolute value of radial stress gradually decreases over time; when the parameters η1 and η2 are fixed, a greater value of σ0 yields a higher absolute value of initial radial stress and convergence value. When η1 is constant, a smaller η2 produces a larger absolute value of the initial and convergence values. Thus, the smaller is the η2, the better is the support effect. η2 does not have a significant effect on the rate of change, which remains fairly constant.  η1, σ0, and η2, σ0).

Tangential Stress Changes over Time
The tangential stress at the monitoring point ρ = 6 m of the reinforced zone is shown in Figure  8a,b. Figure 8a shows the influence of η2 on the tangential stress in the reinforced zone over time. The

Tangential Stress Changes over Time
The tangential stress at the monitoring point ρ = 6 m of the reinforced zone is shown in Figure 8a,b. Figure 8a shows the influence of η 2 on the tangential stress in the reinforced zone over time. The tangential stress gradually decreases over time and then converges to a stable value. When η 1 is fixed, a larger η 2 generates a larger initial and stable values of the tangential stress. Figure 8b shows the influence of η 1 on the reinforced zone tangential stress over time. The tangential stress also decreases over time and finally converges to a stable value. The constant value is independent of η 1 , and a smaller η 1 yields a smaller initial value but a larger change rate. Figure 8c shows the influence of η 2 on the tangential stress of the original zone over time: η 1 influences not only the initial value but also the stable value; a smaller η 2 yields greater initial and stable values of the tangential stress. Figure 8d shows the influence of η 1 on the tangential stress in the original zone over time. The tangential stress also decreases over time, and η 1 influences only the initial value. The tangential stress finally converges to a certain value. The constant value is independent of η 1 , and a smaller η 1 produces not only a greater initial value but also a larger change rate. tangential stress gradually decreases over time and then converges to a stable value. When η1 is fixed, a larger η2 generates a larger initial and stable values of the tangential stress. Figure 8b shows the influence of η1 on the reinforced zone tangential stress over time. The tangential stress also decreases over time and finally converges to a stable value. The constant value is independent of η1, and a smaller η1 yields a smaller initial value but a larger change rate. Figure 8c shows the influence of η2 on the tangential stress of the original zone over time: η1 influences not only the initial value but also the stable value; a smaller η2 yields greater initial and stable values of the tangential stress. Figure 8d shows the influence of η1 on the tangential stress in the original zone over time. The tangential stress also decreases over time, and η1 influences only the initial value. The tangential stress finally converges to a certain value. The constant value is independent of η1, and a smaller η1 produces not only a greater initial value but also a larger change rate.  Figure 9a shows the displacement of the rock mass along the radial distance when the time is zero. When the radial distance increases, the displacement exponentially decreases and converges to zero along the depth of the rock mass. In the moment of excavation, when η2 is fixed, a smaller η1 causes a larger displacement at any position. If η1 is fixed, a smaller η2 yields a larger displacement. Figure 9b shows the displacement when t = 4 d along the radial distance. The amount of displacement significantly changes from the excavation moment. Figure 9c shows the displacement evolution rule of ρ = 7 m over time, the displacement gradually increases and exponentially converges to a fixed value. A smaller η1 produces larger absolute convergence values. A smaller η2 yields larger absolute convergence values. Therefore, important support measures should be taken when η1 or η2 is small.  Figure 9a shows the displacement of the rock mass along the radial distance when the time is zero. When the radial distance increases, the displacement exponentially decreases and converges to zero along the depth of the rock mass. In the moment of excavation, when η 2 is fixed, a smaller η 1 causes a larger displacement at any position. If η 1 is fixed, a smaller η 2 yields a larger displacement. Figure 9b shows the displacement when t = 4 d along the radial distance. The amount of displacement significantly changes from the excavation moment. Figure 9c shows the displacement evolution rule of ρ = 7 m over time, the displacement gradually increases and exponentially converges to a fixed value. A smaller η 1 produces larger absolute convergence values. A smaller η 2 yields larger absolute convergence values. Therefore, important support measures should be taken when η 1 or η 2 is small.

Conclusions
During tunnel excavation, the stresses and deformation of rock mass and the axial force of rockbolt will change over time. This paper established a rheological model of the rock mass supported with DMFC rockbolts. The main findings can be summarized as follows: (1) The axial force of DMFC rockbolts are positively correlated with the support force at the excavation face in a tunnel, and the greater the convergence value of the axial force is, the better the support effect. In addition, the greater the pre-tension of rockbolt, the better the reinforcement effect, however, when the pre-tension is too large, the rock bolt support effect will not increase significantly. (2) The η1 of the three-dimensional Burgers model influences the early support effect, η2 of three-dimensional Burgers model affects both the early and the ultimate reinforcement effect. In addition, there is a significant negative correlation between rock mass displacement and η1 or η2. Therefore, important support measures should be taken when η1 or η2 is small. (3) In this paper, the interaction model elastic solutions were solved based on the distributed force model. However, the axial force of the rockbolt is similar to the concentrated force to rock mass (Bobet 2006), Hence, a more suitable theoretical model should be explored to solve the coupling model in future. (4) Continuously Mechanically Coupled (CMC) rockbolts applications are more extensive than DMFC rockbolts, the reasonable and simplified method for the theoretical model of CMC rockbolts can be further studied based on this model, which lays a foundation of the preliminary research for solving the theoretical model of CMC rockbolts.

Conclusions
During tunnel excavation, the stresses and deformation of rock mass and the axial force of rockbolt will change over time. This paper established a rheological model of the rock mass supported with DMFC rockbolts. The main findings can be summarized as follows: (1) The axial force of DMFC rockbolts are positively correlated with the support force at the excavation face in a tunnel, and the greater the convergence value of the axial force is, the better the support effect. In addition, the greater the pre-tension of rockbolt, the better the reinforcement effect, however, when the pre-tension is too large, the rock bolt support effect will not increase significantly. (2) The η 1 of the three-dimensional Burgers model influences the early support effect, η 2 of three-dimensional Burgers model affects both the early and the ultimate reinforcement effect. In addition, there is a significant negative correlation between rock mass displacement and η 1 or η 2 . Therefore, important support measures should be taken when η 1 or η 2 is small. (3) In this paper, the interaction model elastic solutions were solved based on the distributed force model. However, the axial force of the rockbolt is similar to the concentrated force to rock mass (Bobet 2006), Hence, a more suitable theoretical model should be explored to solve the coupling model in future. (4) Continuously Mechanically Coupled (CMC) rockbolts applications are more extensive than DMFC rockbolts, the reasonable and simplified method for the theoretical model of CMC rockbolts can be further studied based on this model, which lays a foundation of the preliminary research for solving the theoretical model of CMC rockbolts.
Author Contributions: All the authors contributed to publishing this paper. W.H. and G.W. contributed to the formulation of the overarching research goals and aims; C.L. contributed to the theoretical formula derivation of the model; H.L. and K.W. provided language and pictures supports.
Funding: This research was funded by National Natural Science Foundation of China (No. 51479108).

Conflicts of Interest:
The authors declare no conflicts of interest.