Strain-degradation analysis of tunnel structures based on twin-shear unified strength criterion

This paper focuses on the stability analysis of tunnel structures with a circular hollow section based on the with a strain-degradation model using unified strength theory. A simplified numerical approach is proposed for analyzing the elasto-plastic behavior of surrounding rocks. A constitutive equation is proposed based on the unified strength theory and the strain degradation behavior was taken into account as well with deviatoric plastic strain chosen as the degradation parameter. Meanwhile, the plastic potential function can be obtained with the associated flow rule, by which the surrounding rocks can be classified as plastic zone and elastic zone. The solution in plastic zones adopted the differential method subdividing surrounding rock into infinitesimal concentric annuli, in which the radial stresses decrease monotonically with a decrease in radial coordinate. The relationship of the increment of stresses and the increment of strains in adjacent annuli can be derived from equilibrium, constitutive and geometrical equations. The numerical solution of each annulus is calculated from the outmost annulus at the elastic-plastic interface. In addition, the impact of the intermediate principal stress factor and the critical degradation parameter on the solution, and the factors of influence on plastic radius are discussed.


INTRODUCTION
Analysis of stresses and displacements around circular hollow sections in isotropic rock masses has been one of the critical issues in civil engineering.The analytical solutions can be applied in many fields, such as excavation calculation of shield tunnel, the support computation of large drilling engineering and some underground structures.These technical problems necessitate the development of elasto-plastic methods to analyze the stability of underground structures (Brown et al. 1983, Alonso et al. 2003, and David et al. 2014).
A literature review has shown that theoretical analysis methods of ground reaction curves for underground openings excavated in elastic-plastic materials are diversified and have been fully developed (Alonso et al. 2003, Brown et al. 1983, Carranza-Torres and Fairhurst 1999, 2000, Chen et al. 1999, Duncan Fama 1993, Pan and Edwin 1996, Park and Kim 2006, Sharan 2003, 2005, 2008).Those numerical methods provide an accurate prediction of the displacements around a circular opening in the elastic-perfectly plastic rock mass and propose analytical solutions for practical applications.However, much more attention should be paid to the strain degradation, which describes a gradual degradation behavior of materials in the post-failure region.Most of the existing applications are based on the elasticbrittle-plastic model and elastic-perfectly plastic model that are only responsible for the hard rock and soft rock, respectively.The significance of the strain degradation in rock masses has been found and the study of circular tunnels was further carried out by considering the effect of strain degradation (Lee and Pietruszczak 2008, Bai et al. 2012, Guan and Bai 2018, Wang and Qu 2018, Zhang et al. 2018, 2019).In addition, a new formulation for the analysis of degradation was proposed by David et. al (2018), which was based on the extended lumped damage mechanics, introducing ideas of fracture and damage mechanics into the concept of plastic hinge.
The Mohr-Coulomb (M-C) and Hoek-Brown (H-B) failure criteria assumed that rock strength is not significantly influenced by the intermediate principal stress.However, some other researches showed that the rock strength is generally correlated to the intermediate principal stress, which proves the limitation of the M-C and H-B failure mechanisms (Zhao et al. 2018, Vachaparampil and Ghassemi 2017, Lade and Rodriguez 2014, Rodriguez and Lade 2013, Yang et al. 2016, Prashant and Penumadu 2004, and Xu et al. 2017).Yu (2018) proposed a unified strength theory (UST) to account for the influence of the intermediate principal stress and shear stress by introducing a parameter b related to the intermediate principal stress σ2.
There were several studies focused on the approaches related to the theory of UST.Wang et al. (2005) proposed that the effects of various strengths on the dynamic load-carrying capacity are greater than those in the static plastic limit state based on the responses of a circular plate using UST.A formula for computing the ultimate bearing capacity of a shallow foundation based on UST was given by Fan et al. (2005).The results indicated that the strength of foundations is greater than the results calculated of M-C when considering the intermediate principal stress.In addition, the UST has been employed to predict the strength of concrete-filled steel tube (CFST) structures (Guan et al. 2019).The formula of bearing capacity and shear capacity of CFST connection was improved based on the UST (Wang et al. 2015, andShi et al. 2016).Most of the failure criteria are included in the UST while b takes different values (Wang and Qu 2018) and the UST has great advantages to be adopted for analyzing the stress state of rocks composed of various properties.The studies on the circular openings based on the UST have been implemented, and the stresses and displacements of the plastic zone were given (Ma et al. 2009, andZhang et al. 2010).Nonetheless, they neglected the influences of the degradation of strength parameters leading a limitation to the materials with strain softening.
This paper develops a new numerical solution for predicting the strength, deformation and the strain degradation behaviors of circular hollow rock masses with the isotropic property.In the elastic zone, strength parameters can be assumed as constant values, and then linearly correlate to deviatoric plastic strain after rock yielding.The stress and strain fields are given by Lamé solution in the elastic zone which will serve as a priori in the next step of the plastic zone.The plastic zone can be subdivided into a finite number of concentric annuli whose thicknesses are not constant but determined internally to comply with the equilibrium equation.The radial stress corresponding to minor principal stress symmetrically decreases from the yield surface to the excavation boundary of each annulus.Then the hoop stress is calculated in terms of the unified failure criterion and increment of plastic strains can be obtained successively through the compatibility equations in each plastic ring starting from the outmost one, respectively.Furthermore, the proposed UST approach is calibrated by comparing the results of numerical analysis (Pan et al., 2010) with the traditional solution using the generalized failure criteria.Finally, the parametric analyses are performed to quantify the influences of the intermediate principal stress on the critical deviatoric plastic strain, strength parameters and supporting pressure in the plastic zone.

Constitutive Relationship
Figure 1 shows a circular hollow tunnel structure excavated to anisotropic rock mass with radius L. This area sustains the hydrostatic stress field σ0 prior to the excavation.When the internal supporting force pi is smaller the limit state pic, a plastic zone will be triggered around the hollow area expressed by the plastic radius Rp.Furthermore, the plastic zone will consist of a degradation area and residual area with a radius of Rs when the strain degradation behavior is accounted for, as shown in Figure 1.As above mentioned, the rock masses are assumed to be governed by the UST yield criterion given by Yu ( 2018) The relationship among the shear strength τ0, the uniaxial tensile strength σt and the uniaxial compressive strength σc can be determined as follows: For solid materials, if the equation (1) and equation (2) are represented by the principal stress, the rock cohesion C and rock friction angle φ can be written as where σ1, σ2 and σ3 indicate the major, intermediate and minor principal stresses, respectively.b is a UST parameter especially introduced in the UST to take account of the influence of the intermediate principal stress on various failure criteria (0≤b≤1).

Associated flow rule
The plastic deformation law obtained by taking the yield criterion as the potential function is called the associated flow rule, in which the development of plastic deformation is proportional to the normal direction of the yield surface.This hardening process is called isotropic hardening being widely used for elastic-perfectly plastic materials.In this paper, the associated flow rule was adopted.However, due to the generality of the unified strength theory, the non-associated flow rule can also be established by changing the parameter b to a different value with the yield criterion.Therefore, the associated flow rule for the united theory can be simplified to determine whether parameter b is consistent between the yield criterion and the potential function, and further research is needed in this area.In the plane strain state, the intermediate principal stress σ2 is assumed as half of the maximum principal stress σ1 and the minimum principal stress σ3 (Xu andYu 2006, andYu 1998).Therefore, It is apparent to see equation ( 7) should be adopted.If the equation ( 7) is converted to a polar coordinate system, so the plastic potential function becomes where X and Y are the strength parameters defined in terms of the friction angle φ (η), cohesion C (η) and the UST parameter b as follows:

Evolution law of strength parameters
Strain degradation behavior describes a gradual degradation mode of materials and can be used for strain degradation rock mass.The strain degradation model can be understood as the transition from the elastic-brittle-plastic model to the elastic-perfectly plastic model shown in Figure 2, which shows a gradual loss of the strength capacity distinguished with the rapid decline of the elastic-brittle-plastic model and considerable plastic deformation of the elastic-perfectly plastic model.The degradation parameter, which is adopted to evaluate the strain degradation behavior, is defined as η = γ p = εθ p -εr p , where εθ p and εr p are the plastic parts of the hoop and radial strains, respectively.The critical Latin American Journal of Solids and Structures, 2020, 17(2), e251 5/21 deviatoric plastic strain η* dominates the size of the strain degradation region and is related to the rock quality.Furthermore, it should be noted that the critical deviatoric plastic strain η* is a direct reflection of loss modulus, increasing with the decrease of loss modulus.In particular, it is supposed to be the elastic-brittle-plastic rock while η* =0 and elastic-perfectly plastic with η* =∞.  and Structures, 2020, 17(2), e251  6/21 As shown in Figure 3, the strain degradation region appears once the deviatoric plastic strain η is greater than zero.Subsequently, the strength parameters gradually transform with the accumulation of plastic strain.Therefore, the evolution law of strength parameters can be established based on the degradation parameter η shown in equation ( 12) and equation ( 13).where Cp and Cr represent the peak and residual strengths respectively, similar to the associated φp and φr.η* is the critical deviatoric plastic strain from where the residual behavior is initiated.
The linear relationship is utilized to simplify the problem.From the physical point of view, the occurrence of micro-cracks during the rock excavation increases barriers between rock masses which leads to the increase of the friction angle.On the other hand, the formation and propagation of micro-cracks affect the integrity of rock mass, hence the decrease of cohesion.There obtained based on the cohesion weaken friction strengthen (CWFS) model.

Stress and strain fields of the elastic zone
In the elastic zone, the initial strain can be obtained from Hooke's law as follows: Latin American Journal of Solids and Structures, 2020, 17(2), e251 7/21 Based on the assumption of strain remains in-plane 0 0 z   , the hydrostatic stress field can be treated as in elasticity σr 0 = σθ 0 = σ0.Therefore, Equation ( 14) can be rewritten as: The stress, strain and displacement relationship can be obtained through the equilibrium equation, constitutive equation, and geometric equation as follows: When introducing the boundary condition ρ = Rp, σr = σR, ρ = ∞, σr = σ0 into Equations ( 16), ( 17), ( 18), ( 19) and ( 20), the stress and strain fields can be obtained by: Latin American Journal of Solids and Structures, 2020, 17(2), e251 8/21 The stability of an underground space structure is measured by the relative increment of stress and strain.The strain condition in the elastic zone can be expressed by: In the elasto-plastic interface, introducing the boundary conditions and Equation ( 23) results in the equation below, which will serve as the initial value for the subsequent calculation in the plastic zone.
The analytical expression σR is available by solving the following equation It should be noted that the critical pressure pic is equal to R  and the plastic zone is formed only when the initial pressure pi is lower than a critical value pic.

Stress and strain fields of the plastic zone
It is assumed that the plastic zone is composed of n concentric annuli as shown in Figure 4, where the i th annulus is bounded by two circles of the normalized radii ρ(i-1) and ρ(i).The radial stress σr along each ring decreases monotonically from σR when ρ(0) = Rp to pi when ρ(n) = L.The incremental radial stress can be calculated as follows: Latin American Journal of Solids and Structures, 2020, 17(2), e251 9/21 The corresponding hoop stress derived by Equation ( 9) becomes: Then the hoop stress increment can be calculated: Figure 4: Plastic zone with a finite number of annuli.

Stress-strain relationship
The strain condition of the plastic zone is given by the flow rule: The relationship between the radial strain p r d and hoop strain p d   is obtained by eliminating d of Equation (39):

Geometric formula
The compatibility equation for the plane-strain problems is obtained by (Brown et al. 1983): Latin American Journal of Solids and Structures, 2020, 17(2), e251 11/21 where strain is composed of elastic and plastic strain:      , so the Equation ( 41) can be rewritten as: The relationship between e r  and e   is solved by Equations ( 14) and ( 28): ( ) According to the difference method, formula (44), ( 45) and ( 46) are available: Hence the plastic strain increment p r   and p    are calculated as: Then, the total strain and stress at the i th annulus are obtained as: Latin American Journal of Solids and Structures, 2020, 17(2), e251 12/21 Figure 5 plots the detailed calculation process of stress, strain and displacement around the circular tunnel.First, the parameters (L, E, ν, b, Cp, Cr, φp, φr, σ0, pi, η*, n, σR) which are determined by the properties of the material or can be directly calculated are input.Then the solution in the elastic zone can be obtained as well as in the elasto-plastic interfaces.Whether to enter the residual zone needs to be estimated, which is dominated by η*, before resolving the stress and strain in the plastic zone, and the procedure ends if i > n in the excavation face.

Convergence analysis
Numerical simulations of physical processes generally involve solving some differential equations on a computational domain too complicated to solve analytically.Therefore, a numerical solution is an effective and viable option.The Finite Difference Method (FDM) employed in this paper is a way to solve differential equations numerically, which is not the only implement method, alternatives include the Finite Volume Method (FVM) and Finite Element Method (FEM), and also various mesh-free approaches.However, FDM has been widely used not only because of the controllable precision but also to both derive and implement numerically simply.The convergence is closely related to the divided number of the plastic region.Figure 6 compares the convergences of the normalized displacement for M-C rock mass while the finite element numbers are 30, 300 and 3000, respectively.The soft rock (Park and Kim 2006) was used and the detailed data shown in Table 1.   and Structures, 2020, 17(2), e251 13/21 Higher precision can be reached with the smaller division element as well as better continuity based on Figure 6.However, this improvement in precision shows nonlinear growth while the number elements for the plastic zone changes accordingly.The rapid reduction of the normalized displacement exists when the divided number increases from 30 to 300, compared with imperceptible changes from 300 to 3000.Therefore, a suggestive value of the divided number should be clarified with the consideration of precision and computational efficiency.
Figure 7 shows the sensitivity analysis between convergence and computational efficiency.It can be seen that the radius of the plastic zone tends to be stable while the divided element number exceeds 100 and infinite approximates the analytical solution.Figure 7 (b) presents a tendency of computational efficiency with increasing divided number n, in which the computing time is approximately exponential with the number of elements.The number of divided elements is suggested between 200~400 that can meet the high precision and high-efficiency computations.In the following, the element number of 300 will be adopted for parametric analysis.

Verifications for M-C and H-B rock masses
The M-C and H-B failure criteria have been widely applied in rock mechanics from which many new models were verified based on experimental versus numerical comparisons.Park and Kim (2006) adopted the linear M-C model to investigate the stress-strain responses around the circular hollow section, and a generalized nonlinear model based on the H-B failure criterion was proposed (Sharan 2008).However, their models are not capable to incorporate the influence of the intermediate principal stress and cyclic effect (Chen et al. 2017).The proposed model with a unified failure criterion is equivalent to the one with the M-C criterion when b = 0, which neglects the intermediate principal stress.This model can be therefore verified by comparing its analysis results with that of the nonlinear model using the H-B failure criterion.
Figure 8 compares the normalized displacement between the analytical results and the exact solution by Park and Kim 2006).The basic mechanical properties of this experiment are presented in Table 1.In order to calculate the displacement surrounding the hollow portion of the elastic-perfectly plastic rock masses, the critical deviatoric plastic strain η* can be given according to the limit state shown in Figure 2 (c).Figure 8 compares the results of the proposed method with test results by Park and Kim (2006), reasonable agreements can be confirmed that the proposed analytical approach enables to accurately predict the displacement distribution of M-C rock masses.
Latin American Journal of Solids and Structures, 2020, 17(2 Distributions of radial displacements u in the plastic zone around the opening are shown in Figure 9, in which the results obtained by using the new model are compared with the Sharan's analytical solutions (Tu et al. 2018, Lee and Pietruszczak 2008, and Zhang et al. 2012).Table 2 summarizes the material properties of the brittle-plastic rock masses where mb, s and a are the parameters in the generalized H-B failure criterion; mbr, sr and ar are the residual values for mb, s and a, respectively.A geological strength index (GSI) is introduced in the H-B failure criterion for rock masses.The equivalent shear strength parameters, i.e., Cp, Cr, φp and φr are the constants determined by linear regression analysis using the RocLab Program (Zhang et al. 2012).η* =0 is adopted with the elastic-brittle-plastic rock masses.The error as shown in Figure 9 is more notable than that in Figure 8, which may be caused from the error between the UST material properties (Cp, Cr, φp and φr) with H-B material properties (mb, s, a, mbr, sr and ar).The detailed values of the error are presented in Table 3, from which the proposed approach based on the UST demonstrates reasonable accuracy compared to the one with the H-B failure criterion.Thus, the great consistency is shown in Figure 8 and Figure 9 indicating the unified solution includes the results of H-B and M-C when b=0.   and Structures, 2020, 17(2), e251 15/21

Model verification
An elastic-perfectly plastic constitutive model based on the UST was proposed by Pan et al. (2010), who compiled this model in the finite element program ABAQUS.The numerical results of a bi-dimensional circular tunnel structure are compared with the analytical results by using the proposed method.The material and geometric properties of the tunnel structure model are summarized in Table 4.The strain degradation was considered by defining the value of η*= ∞ for the numerical simulation.Figure 10 shows the stress distributions obtained by the new analytical model and the numerical results (pi = 0 MPa).The influence of the parameter b (0, 0.5 and 1) on the r/L responses were evaluated.The r/L ratio corresponding to the peak stress slightly decreased when the value of b increases from 1 to 0, and high accuracy was obtained by the proposed analytical methods.
Table 4 Material properties of a two-dimensional deep circular tunnel case (Pan et al., 2010).

Influence of parameters on plastic radius and stress
Parametric analyses were performed to evaluate the influences of the intermediate principal stress on the critical plastic deviation strain, strength parameters (C and φ), and supporting forces in the plastic zone, the validated analytical model in the previous section is utilized to control the variables.

Intermediate principal stress
Soft rock mass (Park and Kim 2006) is adopted and the effect of strain degradation is taken into consideration with the variable η*= 12e-2 shown in Table 5.On the other hand, the intermediate principal stress is controlled by the parameter b (0~1).In this part, the stresses and radial displacements of the circular hollow section using various values of b (0, 0.5 and 1) are analyzed.The actual distributions of the stresses and radial displacements must lie somewhere in between the two extremes represented by b=0 and b=1, respectively.
Figure 11 shows the relation between the distributions of the radial stress σr and hoop stress σθ versus the r/L ratio when applied to demonstrate the sensitivity of the parameter b with the normalized radial displacements (uE/rσ0).
From Figures 11 and 12, it can be seen that the stresses and radial displacements are all influenced by the parameter b, while the hoop stresses are more significant.The peak values in Figure 11 are the hoop stresses at the elastic-plastic interface moving left with the increase in the value of b, which means the reduction in the radius of the plastic area.For instance, the normalized radial displacements reduce greatly with the increase in the value of b.The uE/rσ0 value is decreased by 28.9%, from 1.73 to 1.23, when b=1 compared with that of b=0.It can be concluded that the latent potentialities of rock mass are underestimated while the rock mass is governed by Mohr-Coulomb or Hoek-Brown failure criteria (b=0).Park and Kim (2006)

Critical plastic deviation strain (strain degradation)
The study on critical plastic deviation strain is implemented using the data of Sharan (2008) shown in Table 2. Particularly, four critical plastic deviation strains corresponding to elastic-brittle-plastic (η*=0), elastic-perfectly plastic (η*=∞) and strain degradation behavior (η*=8e-2, η*=12e-2) of a circular opening are compared in Figure 13.It can be found that the r/L ratio corresponding to the peak stress significantly increases for the elastic-brittle-plastic model relative to that for the elastic-perfectly plastic model.From Figure 13, it can be seen that the variation of critical plastic deviation strain is obvious between the plastic zone and elastic zone, where the peak hoop stress is the demarcation point.The radial displacement u of the brittle-plastic model reaches the largest 0.76m, but the radial stress decreases contrarily.Therefore, these comparative results indicate that smaller plastic zone is likely formed in rock masses when incorporating the effect of strain degradation which might correspondingly result in larger instability potential.

Strength parameters C and φ
The Example of Pan (Table 6) is adopted in this subsection to investigate the influence of strength parameters C and φ.Cohesion ranges from 1 MPa to 8 MPa with a constant frictional angle φ (30°) is employed to assess the effect of parameter C. On the contrary, the influence of parameter φ is evaluated with this variable changing from

Supporting compression force
Influences of the supporting force pi on displacement distribution evaluated as illustrated in Figure 15 and analytical models using two values of b (0, 1) are compared.Material and geometric properties of the tunnel structure model are summarized in Table 4 in which the constant supporting force pi is substituted by a variable value ranging from 0 MPa to 10 MPa.The double drop lines illustrate that the normalized radial displacement reduces with the supporting pressure pi increasing.The effects of b and pi are compared and the results indicate that the supporting pressure plays a more important role in maintaining the stability of the opening.It can be seen that taking the intermediate principal stress into account decreases the normalized radial displacement around 0.78, which is a contrast with that of 3.02 for pi.However, the influence of b cannot be completely ignored.When the normalized radial displacement is constant (1.7), the supporting pressure will shrink by 33.3%, from 6.3 MPa to 4.2 MPa, considering the intermediate principal stress effect.Thus, the supporting pressure is overestimated while the rock mass is governed by Mohr-Coulomb or Hoek-Brown failure criteria (b=0).

Figure 1 :
Figure 1: Plastic zone form around the circular hollow rock mass.

Figure 2 :
Figure 2: Schematic graph of excavation problem and stress-strain relationship: (a) for elastic-perfectly-plastic rock mass, (b) for strain degradation rock mass, (c) for elastic-brittle-plastic rock mass.

Figure 3 :
Figure 3: Evolution law of strength parameters φ and C.

Figure 5 :
Figure 5: Flowchart for calculating stress and displacement.

Figure 6 :
Figure 6: Normalized displacement distribution with M-C rock mass

Figure 7 :
Figure 7: Sensitivity analysis between convergence and computational efficiency: (a) the influence of divided number on plastic zone radius; (b) the influence of divided number on computational time.

Figure 8 :
Figure 8: Comparison of displacement distribution with M-C rock mass

Figure 9 :
Figure 9: Comparison of displacement distribution with H-B rock mass.

Figure 10 :
Figure 10: Comparisons on the stress distributions of a cylindrical tunnel structure.

Figure 13 :
Figure 13: Elastic-brittle-plastic, elastic-perfectly plastic and strain degradation behavior of a circular opening in UST rock mass.

Table 6 Figure 14 :
Figure 14: Correlation of strength parameters C and φ on displacement distribution.

Figure 15 :
Figure 15: Effect of pi on displacement distribution with different b.

Table 3
Comparisons on standard errors between the proposed model and the model with H-B criterion

Table 5
Data of soft rock provided in