- Email: [email protected]
A numerical large strain solution for circular tunnels excavated in strain-softening rock masses
Computers and Geotechnics 114 (2019) 103142
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
A numerical large strain solution for circular tunnels excavated in strainsoftening rock masses Qiang Zhanga, Hong-Ying Wangb,c, Yu-Jing Jiangb, Meng-Meng Lua, Bin-Song Jianga,
T
⁎
a School of Mechanics & Civil Engineering, State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining & Technology, Xuzhou 221116, China b Geoenvironmental Laboratory, Graduate School of Engineering, Nagasaki University, Nagasaki 8528521, Japan c College of Energy and Transportation Engineering, Jiangsu Vocational Institute of Architectural Technology, Xuzhou 221116, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Circular tunnel Large strain Strain-softening Yield criterion Elastic-plastic analysis
This paper presents a large strain numerical solution of a circular tunnel for the ground reaction curve (GRC) in strain-softening rock masses. Both the stress equilibrium equation and deformation compatibility equation are considered in a current coordinate. A series of examples are employed to validate the proposed solutions, and the results show that: (1) the displacement, softening and residual radii of larger strain solution are smaller than that of small strain solution, and their difference significantly increases with a decrease in the Young’s modulus. (2) the increasing Poisson’s ration and critical shear plastic strain decrease the displacement, softening and residual radii.
1. Introduction The tunnel clearance size and supporting structure parameter are of essential importance in the geotechnical engineering design. A reasonable prediction of the displacement and failure region for the surrounding rock around the tunnel would be the primary influence on the design optimization and economics. In recent years, many studies on elastic-plastic rock masses have been carried out, and many achievements have been obtained [1–10]. Among these studies, scholars mainly pay attention to the Ground Response Curves (GRCs) according to different constitutive models, i.e., the elastic-perfectly (brittle) plastic model [1,2] and strain-softening model [3–8], with numerical and (semi-)analytical methods that follow small strain theory. To improve these models, the deterioration effects of the elastic parameters [9,10], confining pressure [11], seepage [12,13] and out-of-plane stress [14] are further considered. Strictly speaking, those solutions are only appropriate for engineering applications, in which the rocks/structures deform negligibly small, such as a tunnel excavated in the hard rock by assuming that small strain theory is suitable. Unfortunately, the convergence displacement may exceed one-tenth of the initial excavation radius in the surrounding soft strata subjected to high in-situ stress with a significant decrease in the cavity radius. For example, the displacement exceeds 10 ∼ 14% of the initial excavation radius in the Wushaoling tunnel, and that of the Mingyazi tunnel even exceeds 30% of the initial excavation radius, as shown in Fig. 1. The
⁎
displacement predicted using small strain theory deviates significantly from the natural cases because the position of material point is assumed to be fixed, which however has a non-negligible deformation in such large stain cases. Since the supporting pressure is always acting on the deformed excavation radius, considering the stress equilibrium on the current configuration instead of on the initial configuration seems more suitable for large deformation problems. Therefore, the stress and deformation compatibility equations should be always satisfied for the deformed configuration. Considering the mechanical behavior of large strain rock masses, a few attempts have been made to derive the solutions of circular tunnels. Papanastasiou and Durban [15] studied the expansion and contraction behavior of cylindrical cavities in an elastic-perfectly plastic rock mass and derived the differential solutions. Later, Yu and Rowe [16] presented an analytical solution for the problem of a cavity. With the assumption of small strain in the elastic region and considering the out-ofplane stress effect, Vrakas and Anagnostou [17–19] derived the closedform solution of a circular tunnel using the linear Mohr-Coulomb criterion and non-linear envelope criterion. Additionally, a closed-form relationship between the displacement of small and large strain solutions for the cavity contraction/expansion problem was presented considering the general Mohr’s criterion and constant plastic dilatancy. By assuming the elastic rock mass as small stain condition, Park et al. [20] derived the large strain similarity solution for a spherical or circular opening excavated in elastic-perfectly plastic media, and the
Corresponding author. E-mail addresses: [email protected] (Q. Zhang), [email protected] (Y.-J. Jiang), [email protected] (B.-S. Jiang).
https://doi.org/10.1016/j.compgeo.2019.103142 Received 26 July 2018; Received in revised form 2 May 2019; Accepted 17 June 2019 0266-352X/ © 2019 Published by Elsevier Ltd.
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
(a)
(a)
ı
Design contour
Elastic-perfectly plastic u >100cm about 10~14%R
Strain-softening Deformed contour
Elastic-brittle-plastic
İ (b) Elastic-perfectly plastic
-İ3
(b)
Design contour
Strain-softening
u >230cm over 30%R
Elastic-brittle-plastic
İ1 Deformed contour
Fig. 2. Mechanical behavior of rock mass: (a) stress-strain curves, and (b) strain curves.
Fig. 1. Large deformation examples of tunnel contraction: (a) Wushaoling tunnel; (b) Mingyazi tunnel.
Strength drop
Load /kN
comparisons between the small and large strain solutions denoted that the displacement assuming large strain is smaller than that assuming small strain. Considering the large strain behavior of both elastic and plastic rock masses, Zhang et al. [21] derived the differential expressions for a circular opening excavated in elastic-brittle-plastic rock. The examples denoted that the stresses and displacement of the elastic rock mass depend on the elastic parameters, and the difference between the large strain solution and small strain solution may be significant, especially for extremely soft surrounding rock. Guan et al. [22] proposed a finite strain numerical procedure of a circular tunnel in strainsoftening rock mass considering the axial stress effect. Similar to the former study of Vrakas and Anagnostou, the elastic rock mass is assumed following the small strain behavior. Clearly, most of the studies on large strain problems are carried out according to elastic-perfectly (brittle) plastic behavior. For the elasticbrittle-plastic model, a sudden drop in stress immediately happens once the peak stress is reached. Meanwhile, an ideal plastic flow performs when the yield criteria is satisfied for the elastic-perfectly plastic model. The elastic-brittle-plastic model and elastic-perfectly plastic model are particular conditions of mechanical behaviors of rock mass. They are only reasonable for the hard and soft rock mass, respectively. Most of the geotechnical materials generally manifest strain-softening behavior in the post-failure region, whose mechanical behavior is between the above two special cases, as shown in Fig. 2. Fig. 3 shows a typical cyclic loading deformation curve of a rock sample. With the increase in
Envelope curve Test curve
Axial deformation /mm Fig. 3. Typical cyclic loading deformation curve of a rock sample.
deformation, the bearing capacity of the post-failure rock mass gradually decreases. The progressive failure process is very important for stability evaluation of the surrounding rock and optimization of the lining structure in underground engineering. To provide a more reasonable prediction in the GRC and failure region for underground geotechnical engineering, this paper aims at the numerical large strain solutions of circular tunnels excavated in strainsoftening rock masses. The differential expressions are derived for both elastic and plastic rock masses, in which the material position is
2
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
p0
in the Eulerian form as follows:
(
Softening-residual boundary
ıcr r
(
u0
p0
ı0 r0 R0
Rp
⎧ εr = ⎨ εθ = ⎩
u 2 r
(3)
1 [(1 2G 1 [(1 2G
− v )(σr − p0 ) − v (σθ − p0 )] + εrp − v )(σθ − p0 ) − v (σr − p0 )] + εθp
(4)
where G = E/[2(1 + v)] is the shear modulus; E and v are Young’s modulus and Poisson’s ratio, respectively; εrp and εθp are the radial and tangential plastic strains; and εrp ≡ εθp ≡ 0 for the elastic rock mass. With respect to the strain-softening behavior, the softening parameter, which governs the evolution of the material properties, should be determined first. Here, the most widely used plastic shear strain γ p is employed as the softening index, as shown in Eq. (4).
Elastic-softening boundary
p0 Fig. 4. Circular tunnel in an infinite deformable medium.
γ p = εθp − εrp
(5)
The Mohr-Coulomb (MC) plastic potential function is considered (i.e., non-associated plastic flow). According to the plastic potential theory [24], the radial and tangential plastic strains would relate to
deformable. Then, a recursive implementation procedure is proposed, and the proposed solutions are validated by comparisons with numerical solutions. Finally, parametric analyses of Young’s modulus, Poisson’s ratio and the critical plastic shear strain are investigated.
εrp + κ (γ p) εθp = 0
(6)
1 + sin ψ (γ p) , 1 − sin ψ (γ p)
where κ (γ p) = and ψ is the dilation angle. The yield criterion for the plastic rock mass can be expressed as
2. Problem description 2.1. Geometrical relationship
σθ − σr = H (σr , γ p)
Fig. 4 shows a circular tunnel with an initial radius R0 excavated in an infinite isotropic rock mass, and a hydrostatic stress field p0 is imposed throughout the domain prior to the excavation. The assumption of a hydrostatic stress field implies that the body force is neglected, and no rotation of the principal axes occurs for this axisymmetric problem. Then, the problem transforms into a one-dimensional problem in terms of a radial displacement. Additionally, contraction displacement is defined as positive. With the decrease in the supporting pressure σ0 from its initial value p0, the excavation radius gradually converges to r0 with a displacement of u0. For the deformable surrounding rock, the initial position of a material point R is now at the current radius r. Therefore, the radial displacement u of the material point can be defined as
u=R−r
du 2 dr
( ) +⋯ ( ) +⋯
where εr and εθ denote the radial and tangential strains, respectively. In a small strain condition, the higher order terms on the right side of Eq. (3) are neglected, leading to the well-known kinematic equations. The total strain induced by excavation is composed of elastic strain and plastic strain. Based on Hooke’s law, the constitutive equation can be formulated as
Rs
p0
) )
du du 1 ⎧ ⎪ εr = ln 1 + dr = dr − 2 ⎨ ε = ln 1 + u = u − 1 ⎪ θ r r 2 ⎩
(7)
As introduced above, the material parameters vary with the softening parameter. For the linear MC yield criterion, H(σr, γp) can be written as
2c (γ p) cos φ (γ p) 1 + sin φ (γ p) HMC (σr , γ p) = ⎜⎛ − 1⎟⎞ σr + p 1 − sin φ (γ p) ⎠ ⎝ 1 − sin φ (γ )
(8)
where φ is the friction angle and c is the cohesive strength, respectively. For the generalized nonlinear Hoek-Brown (HB) yield criterion, H (σr, γp) becomes
HHB (σr ,
γ p)
= σc
(γ p) ⎛⎜m (γ p) ⎝
(1)
σr + s (γ p) ⎞⎟ σc (γ p) ⎠
a (γ p)
(9)
where σc is the uniaxial strength; m, s and α are three strength parameters, which are related to the geological strength index.
Once the internal supporting pressure σ0 becomes lower than a critical value σrcr , a current plastic region with a radius of Rp develops around the tunnel. In some cases, the plastic region of rock includes another current residual region with a radius Rs besides Rp, in which the material properties reach its residual value.
2.3. Evolution of material parameters
For the axisymmetric plane strain problem, the stress equilibrium equation in the current state can be expressed as follows:
The evolution of the material parameters mainly depends on the lithology and stress state, and it can be obtained through cyclic loadingunloading triaxial compression tests under various confining pressures. For the sake of comparison, a piecewise linear function of γp is assumed in this paper, as shown in Eq. (10).
dσr σ − σθ + r =0 dr r
η (γ p) =
2.2. Governing equations
γp
(2)
where σr and σθ denote the radial stress and tangential stress, respectively. For the present problem, the logarithmic strains can accurately describe the strain-displacement relationship [23]. In the current polar coordinates, the logarithmic radial and tangential strains can be written
⎧ ηp − (ηp − ηr ) γ p ∗ , ⎨ η, r ⎩
0 < γ p < γ p∗ γ p ⩾ γ p∗
(10)
where η denotes any one of the strength parameters, i.e., φ, c, σc, m, s, α and ψ; γp* is the critical plastic shear strain for the residual state of the surrounding rock; and ηp and ηr are the peak and residual values of η, respectively. 3
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
The plastic rock mass should also satisfy the deformation compatibility equation. By using Eqs. (1), (2), (4), (7) and (17), the deformation compatibility equation of a plastic rock mass can be formulated in the original coordinate as follows:
2.4. Boundary conditions Since the supporting pressure is always acting on the excavation boundary, we have σr(r0) = σ0. The stress state at an infinite distance from the tunnel remains undisturbed; thus, it leads to lim σr (r ) = lim σθ (r ) = p0 and lim u (r ) = 0 . The current position of a
r →∞
r →∞
dεθp
r →∞
dR
material point is still unknown before the exact solution is obtained. However, Eq. (1) implies that the material points’ positions are a oneto-one correspondence in the original and deformed coordinates. Then, the boundary condition can also be expressed in the original coordinate as follows:
⎧ R = R0, ⎨ ⎩ R → ∞,
σr = σ0 σr = σθ = p0
R
-
ζH (σr , γ p) H (σr , γ p)/2G + (1 + κ (γ p)) ε p θ e 2GR
1 + sin φ (γ p) ∂HMC (σr , γ p) −1 = 1 − sin φ (γ p) ∂σr
(19)
(11)
∂HHB (σr , γ p) σr = m (γ p) α (γ p) ⎜⎛m (γ p) + s (γ p) ⎟⎞ ∂σr σ ( γ p) c ⎝ ⎠
α (γ p) − 1
(20)
The tangential stress can be easily obtained by substituting the radial stress into the yield criterion. Using Eqs. (1) and (2), the displacement of a plastic rock mass can be obtained as follows: p
u = R {1 − e−[(1 − v)(σθ− p0 ) − v (σr − p0 )]/2G − εθ }
Theoretically, the material parameters at each point should obey the evolution of Eq. (10), and this is the main challenge in the strain-softening analysis. In the finite differential method (FDM), only the integration points exactly satisfy the evolution of the material parameters. The stresses at the other points are estimated by the shape functions and the stresses of the integration points [25]. Therefore, the material parameter for each element is represented by a series of values at the integration points. For the one-dimensional problem of a circular tunnel, the surrounding rock is divided into n annuli, and each annulus is assumed to be isotropic with a constant material parameter, which depends on the plastic shear strain at its outer radius. So far, the stresses and displacements in elastic and plastic surrounding rocks are obtained in a differential equation form. To solve the differential equations of Eqs. (14)–(16) and Eqs. (18)–(21) for the elastic and plastic rock masses, many ordinary differential methods are available, such as the Euler and Runge-Kutta methods. Unlike in the small strain solutions, the stresses depend on the elastic parameters, i.e., Young’s modulus and Poisson’s ratio, suggesting that the explicit stress at the elastic-softening boundary cannot be obtained. In this condition, a numerical method is necessary to solve the differential equations for the elastic and plastic rock masses. To minimize the cumulative error, a fifth-order RungeKutta method is employed, as shown in Eq. (22):
(12)
(13)
By substituting the elastic constitutive equation of Eq. (4) in Eqs. (12) and (13), the stress and deformation compatibility equations of the elastic surrounding rock can be expressed in terms of the radial and tangential stresses as follows:
dσr σ − σr (σθ− σr )/2G e = θ dR R
(14)
dσθ 2G v (σθ − σr ) (σθ− σr )/2G [1 − e (σθ− σr )/2G] + e = dR (1 − v ) R (1 − v ) R
(15)
(21)
4. Numerical implementation procedure
The radial and tangential strains are functions of the displacement with respect to its current position. Therefore, the strains should satisfy the deformation compatibility equation, which can be obtained by substituting the tangential strain into the radial strain in Eq. (3). Then, using the relationship between the original and current coordinates of Eq. (1), the deformation compatibility equation can be expressed in the original coordinate as follows:
35 500 125 2187 11 yi + 1 = yi + ⎛ k1 + k3 + k4 − k5 + k6 ⎞ h 1113 192 6784 84 ⎠ ⎝ 384
Once the stresses are obtained, the displacement of the elastic surrounding rock can be calculated with the tangential geometrical equation given in Eq. (3). Then, the displacement can be expressed in the original coordinate by using Eq. (1), as shown in Eq. (16).
(22)
where h is the step length and
k1 = f (x i , yi ) =
(16)
3.2. Solution in the plastic region By using the same solution approach as that used for the elastic rock mass, the stress equilibrium equation of a plastic rock mass can be reformulated in the original coordinate using Eqs. (1), (3) and (4). Meanwhile, the stress state should satisfy the yield criterion of Eq. (7). By using the yield criterion, the stress equilibrium equation can be formulated as follows:
H (σr , γ p) H (σr , γ p)/2G + (1 + κ (γ p)) ε p dσr θ = 0 e = dR R
)/2G + (1 + κ (γ p)) εθp
(18)
The stress equilibrium equation is expressed in the current (deformed) position, as shown in Eq. (2), but the displacement of the material point is unknown. To solve the present problem, the stress equilibrium equation should be expressed in the original coordinate. Using Eqs. (1) and (3), the stress equilibrium equation can be rewritten in the original coordinate as follows:
u = R {1 − e−[(1 − v)(σθ− p0 ) − v (σr − p0 )]/2G}
p
∂H (σ , γ p)
3.1. Solution in the elastic region
dεθ 1 = (1 − e εθ− εr ) dR R
1 − e H (σr , γ
− 2v + 1; for the MC and HB rock masses, the where ζ = (1 − v ) ∂σr r derivate term can be respectively expressed as follows:
3. Solutions for stress and displacement
dσr σ − σθ εθ− εr e + r =0 dR R
=
d (x i , yi ) dx
(23a)
1 1 k2 = f ⎛x i + h, yi + k1 h⎞ 5 5 ⎝ ⎠
(23b)
3 3 9 k3 = f ⎛x i + h, yi + k1 h + k2 h⎞ 10 40 40 ⎠ ⎝
(23c)
4 44 56 32 k 4 = f ⎛x i + h, yi + k1 h − k2 h + k3 h⎞ 5 45 15 9 ⎝ ⎠
(23d)
8 19372 25360 64448 212 k5 = f ⎛x i + h, yi + k1 h − k2 h + k3 h − k 4 h⎞ 9 6561 2187 6561 729 ⎝ ⎠ (23e)
(17) 4
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
k6 = f ⎛x i + h, ⎝ 9017 355 46732 49 5103 yi + k1 h − k2 h + k3 h + k4 h − k5 h⎞ 3168 33 5247 176 18656 ⎠
annulus at Ri is employed to predict the material properties of the ith annulus with Eq. (10). The tangential plastic strain εθp(−i) is employed as the initial value for the ith annulus for solving the differential equation of Eq. (18). Only in this way can the displacement be continuous; otherwise, the accumulative error will lead to inaccurate results, especially for the elastic-brittle-plastic rock mass. Similar to the calculation method for the elastic-softening radius, the plastic shear strain equals γp* at the softening-residual boundary. For the plastic rock mass, the calculated plastic shear strain is compared with γp* in each step. Once the calculated plastic shear strain is larger than γp*, the inner radius is adjusted to ensure that the plastic shear strain at the inner radius equals γp*. Then, the stress and displacement of each annulus can be recursively calculated until the radius R reaches R0. If the required supporting pressure σ¯r (1) at R = R0 corresponding to the given tangential σ¯θ (n) equals the actual supporting pressure σ0, the calculated results would be corrected. Otherwise, the initial tangential stress σ¯θ (n) is updated until σ¯r (1) = σ0 . Additionally, an adaptive step length method can be used to improve the computation efficiency. For convenience, in this study, a fixed small step length is considered throughout the entire solution to ensure the stability of the solution. The detailed calculation flowchart of the implementation is shown in Fig. 5.
(23f)
To solve the differential equations of the ith annulus, the initial slopes of each variable at the boundary xi should be given first. However, Eq. (11) denotes that the exact boundary condition cannot be accurately given due to the infinite radius. According to Saint-Venant's principle, the deviation would be very small when the boundary condition is replaced by a very large finite radius, i.e., 50R0. Hence, the boundary conditions can be replaced by
⎧ R = R0, ⎨ ⎩ R = 50R 0 ,
σr = σ0 σr = p0
(24)
Then, the present problem transforms into a boundary problem. In the contraction problem of a circular tunnel, the tangential stress is larger than the initial stress in the elastic region because of stress concentration. For the plane strain problem of circular tunnel, the elastic tangential stress is the first principal stress, and the radial stress is the minimum principal stress, that is σθ > σr according to the Lame’s solution. Thus, an initial tangential stress σ¯θ (n) > p0 at R = 50R0 is assumed, and the slope k1 can be obtained for the radial stress and tangential stress. The slope k1 is used to predict the dependent variables at xi + h/5; at this point, the slope k2 can also be determined. Other slopes, i.e., k3−k6, which are in the region (xi, xi+1), can be continuously obtained. Slopes k1−k6 are combined to predict the variables at xi+1 with Eq. (22). Therefore, all of the stresses and displacements for the elastic region can be recursively determined. In the plastic region, the variables are the radial stress and tangential plastic strain, and the radius Rp is still unknown. At the elasticsoftening boundary, the elastic stresses should satisfy the peak yield criterion due to zero plastic strain. For each annulus, the calculated elastic stresses are verified by the peak yield criterion. If the calculated stress (σr(i), σθ(i)) satisfies the peak yield criterion, the elastic-softening boundary is located in the range of (xi, xi+1). The step length is adjusted by dichotomy until the stress (σr(i), σθ(i)) is located on the peak yield surface. Then, the differential equations given in Eqs. (17), (7) and (18) are used to calculate the stresses and tangential plastic strain for the plastic rock mass. With the aid of Eq. (6), the radial plastic strain can be obtained. By substituting the radial and tangential plastic strains into Eq. (5), the plastic shear strain at the inner radius can be obtained. It should be noted that the radial stress and displacement should be continuous, and the stress state of each annulus satisfies the yield criterion. However, the material properties of each annulus are different for the softening rock mass, which may cause a sudden change in the tangential plastic strain at the boundaries of the adjacent annuli. For the (i + 1)th and ith annuli in the softening region, by using the displacement expression of Eq. (21), the tangential plastic strain at R = Ri for the (i + 1)th and ith annuli will be related by
εθp(−i) = εθp(+i) +
1 [(1 2G(i + 1)
−
1 [(1 2G(i)
5. Example verifications 5.1. Verifications for an elastic-perfectly (brittle) plastic rock mass To validate the proposed solutions, the typical input data for the MC and HB rock masses used by Ogawa and Lo [26] are utilized and listed in Table 1. The solutions are validated with the methods proposed by Zhang et al. [21] for both the MC and HB rock masses. For accuracy, a small step length of 0.01R0 and 0.001R0 for the elastic and plastic rock masses are considered, respectively, and an acceptable tolerance of 0.001p0 is employed for the supporting pressure. When a significantly large or small critical plastic shear strain, i.e., γp* = 10 or 1e-5, is employed, the proposed solutions would transform into the solutions of elastic-perfectly plastic and elastic-brittle-plastic rock masses, respectively. The dimensionless plastic radius Rp/R0 and dimensionless displacement u0E/p0R0 are calculated, as listed in Table 2. The calculated results are compared with those using the method of Zhang et al., and the comparison results denote that the proposed solutions are in agreement, with a relative error less than 0.6% for both the dimensionless softening radius and displacement. 5.2. Verifications for a strain-softening rock mass To investigate the strain-softening behavior of a rock mass with large strains, two sets of input data are presented for the MC and HB rock masses, as visualized in Fig. 6 and Fig. 7, respectively. Several values of step length are employed to investigate the convergence of the proposed solutions. For the MC rock mass, the dimensionless displacement u0/R0 is 0.3195 when the step lengths are 0.01R0 and 0.001R0 for the elastic and plastic regions, respectively, and that of 0.3222 corresponds to the step lengths of 1e-4R0 and 1e-5R0. The relative error between the two conditions is less than 8.5‰. Therefore, the results can be considered the actual solutions when the step lengths of 0.01R0 and 0.001R0 are employed for the elastic and plastic rock masses, respectively, which are employed in the following analyses. Fig. 6 shows the GRC and evolutions of Rp and Rs for the MC rock mass. The calculation results are compared with the FDM solutions using the Fast Lagrangian Analysis of Continua (FLAC) software. Since the calculation accuracy of FDM solutions depends on the mesh generation, a minimum element size of 0.1R0 in the concerned region of R < 5R0 is used in the FLAC. Although smaller element sizes result in more accurate solutions, it also causes a localization effect when the element size is small enough, and the softening and residual region would be a series of belt regions [27]. Therefore, an average value of
− v(i + 1) ) H (σr (i), γ(pi + 1) ) + (1 − 2v(i + 1) )(σr (i) − p0 )] − v(i) ) H (σr (i), γ(pi) ) + (1 − 2v(i) )(σr (i) − p0 )] (25)
εθp(−i)
εθp(+i)
and are the tangential plastic strains at Ri for the ith and where (i + 1)th annuli, respectively. Similarly, the plastic shear strains at Ri for the (i + 1)th and ith annuli are also different. Strictly speaking, taking the material parameters of the middle point as the representative values of each annulus is more reasonable. However, this will cause more complex expressions, which should be solved by an iterative algorithm. Zhang et al. [24] studied the influence of four different representative points on the calculation results for the small strain-softening problem, and the four results tend to be uniform as the annuli number increases. To reduce the difficulty of the calculation, the plastic shear strain of the (i + 1)th 5
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
Table 2 Computed results of the dimensionless radius and displacement (data in parentheses are the solutions by using the methods of [21]). Rock mass
Calculation scheme
Rp/R0
MC rock
ebp ep ebp ep
1.7487 1.1642 1.2362 1.6110
HB rock
u0E/(p0R0) (1.7576) (1.1649) (1.2364) (1.6123)
11.9279 (11.9596) 1.5583 (1.5666) 1.1861 (1.1883) 5.9844 (6.0176)
Note: ‘ebp’ denotes the elastic-brittle-plastic behavior, and ‘epp’ denotes the elastic-perfectly plastic behavior.
(a)
0.25
This paper FDM
0.20 R0=1.0, p0=1MPa, E=30MPa, v=0.3
σ0/p0
0.15
p*
γ =0.15, cp=0.2MPa, cr=0.02MPa ϕp=40°, ϕr=20°, ψp=10°, ψr=5°
0.10 0.05 0.00 0.0
0.1
0.2
0.3
0.4
u0/R0
(b)
0.24 0.21 0.18 Rp/R0
σ0/p0
0.15
Rs/R0
0.12 0.09 0.06 0.03 0.00 1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.4
2.6
R/R0 Fig. 5. Detailed flow chart for the sequence of calculations.
Fig. 6. Large deformation solution for the MC rock mass (a) GRC; (b) evolution of Rs and Rp.
Table 1 Geometry and material property parameters. MC rock Radius of opening, R0 (m) Initial stress, p0 (MPa) Supporting pressure, σ0 (MPa) Young’s modulus, E (GPa) Poisson ratio, v(-) Dilation angle, ψp = ψr (Deg) cp (MPa) φp (Deg) cr (MPa) φr (Deg)
the displacement at the excavation surface is taken as the displacement at R0, and the evolutions of Rp and Rs are not listed due to the bifurcation. With the decrease in the supporting pressure, the displacement, softening radius and residual radius gradually increase and reach maximum values when the supporting pressure is completely released. The GRC and evolutions of Rp and Rs for the large strain rock mass are similar to those for the small strain rock mass. The dimensionless displacement u0/R0 is 0.319 and 0.323 for the proposed solution and FDM solution, respectively, with a relative error of 4.4%. The maximum softening and residual radii are 2.25R0 and 1.87R0, respectively. Clearly, the displacement of the proposed total strain solution is slightly smaller than that of the FDM solutions with an incremental velocity approach. This is in agreement with Park’s studies on large strain solutions for an elastic-perfectly plastic rock mass [20]. For the generalized HB rock mass, a composite flow rule with an HB
HB rock 1 1 0 5 0.2 30 0.276 35 0.055 30
Radius of opening, R0 (m) Initial stress, p0 (MPa) Supporting pressure, σ0 (MPa) Young’s modulus, E (GPa) Poisson ratio, v(-) Dilation angle, ψp = ψr (Deg) Uniaxial strength, σc (MPa) mp sp mr sr
1 1 0 5 0.2 30 50 0.2 1e-4 0.05 1e-5
6
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
(a)
0.40 0.35
Proposed solution FDM
0.30
R0=1.0 m, p0=1 MPa, E=40 MPa, v=0.3
σ0/p0
0.25
σp=10 MPa, σr=8 MPa, ap=0.5, ar=0.55
0.20
σ0/p0
(a)
mp=0.5, mr=0.05, sp=1e-4, sr=1e-5
0.15
p*
ψp=5°, ψr=3°, γ =0.15
p0 = 1.0 MPa
0.25
p0 = 2.0 MPa
0.20
p0 = 5.0 MPa
0.05
0.05
0.1
0.2
0.3
0.4
p0 = 1.5 MPa p0 = 3.0 MPa p0 = 10.0 MPa
0.15 0.10
0.00 0.0
0.5
u0/R0
0.2
0.4
0.6
0.8
1.0
u0/R0
(b)
0.40 0.35
0.35 0.30
0.30
0.25
Rp/R0 Rs/R0
0.25
0.20
0.20
σ0/p0
σ0/p0
0.30
0.10
0.00 0.0
(b)
0.35
0.15
0.15 0.10
0.10
0.05
0.05 0.00 1.0
1.5
2.0
2.5
3.0
0.00 1.0
3.5
1.2
1.4
1.6
1.8
2.0
2.2
2.4
2.6
Rp/R0
R/R0
(c)
Fig. 7. Large deformation solution for the HB rock mass (a) GRC; (b) evolution of Rs and Rp.
0.35 0.30
potential function is employed in the FLAC software. To compare these results with the proposed solutions, a new constitutive model based on the MC potential function, in which the radial and tangential plastic strains follow Eq. (6), is developed. When the supporting pressure is completely released, the average dimensionless displacement u0/ R0 = 0.432, and the proposed solution results in u0/R0 = 0.399, Rp/ R0 = 3.122 and Rs/R0 = 2.067, whose displacement is also slightly smaller than that of the FDM solutions with a percentage error of 8.3%, as shown in Fig. 7. In summary, the proposed large strain solutions are in agreement with that of the FDM solutions for both the MC and HB strain-softening rock mass.
0.25
σ0/p0
0.20 0.15 0.10 0.05 0.00 1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.4
Rs/R0
6. Parametric analyses
Fig. 8. Evolutions of u0, Rp and Rs under various p0 and σ0 for MC rock mass: (a) u0; (b) Rp and (c) Rs.
In this section, the influences of supporting pressure and hydrostatic stress on u0, Rp and Rs are studied, and the material property parameters for MC and HB rock mass are the same as that used in Section 5.2 as listed in Figs. 6 and 7, respectively. In this section, those material parameters are considered.
with a decreasing rate and gradually approaches to the limit of R0. However, Rp and Rs do not monotonously increase with the increasing p0. Taking the condition of σ0 = 0 as an example, Rp = 2.22 R0 and Rs = 2.20R0 corresponding to p0 = 2.0 MPa, which are smaller than those of 2.39R0 and 2.22R0 corresponding to p0 = 1.5 MPa. With a further increase in p0, Rp and Rs still decrease. Thus, when the position changes in material point is considered, a very large ground stress may induce a small plastic and residual region. A similar evolutions in u0, Rp and Rs under various σ0 and p0 could be found for HB rock mass as shown in Fig. 9. When p0 > 2.0 MPa, an evident inflection point occurs, and from then on the increasing rate of u0, Rp and Rs become much small. Additionally, with the further increase in p0 (p0 > 2.0 MPa), Rp
6.1. Effects of supporting pressure and hydrostatic stress The evolutions of u0, Rp and Rs for MC and HB rock masses with various of supporting pressure and hydrostatic stress are shown in Figs. 8 and 9, respectively. For MC rock mass, u0, Rp and Rs gradually increase with the decrease in σ0 for a constant p0. Under the condition of p0 < 1.5 MPa, Rp and Rs significantly increase when σ0 < 0.1p0. However, the increasing rate for σ0 < 0.1p0 becomes much small when p0 > 1.5 MPa. On the other hand, with the increase in p0, u0 increases 7
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
(a)
(a)
0.7
1.0
p0 = 1.0 MPa p0 = 1.5 MPa
0.6
0.8
p0 = 2.0 MPa
0.5
p0 = 3.0 MPa p0 = 5.0 MPa
0.6
p0 = 10.0 MPa
σ0/p0
σ0/p0
0.4 0.3 0.2
E=0.1 MPa E=0.3 MPa E=3.0 MPa E=30 MPa E=3.0 GPa
0.4 0.2
0.1 0.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.0
1.2
(b)
1.0
0.6
0.8
1.0
1.0
E=0.1 MPa E=0.4 MPa E=4.0 MPa E=40 MPa E=4.0 GPa
σ 0/ p 0
0.6
0.4
σ0/p0
0.8
0.8
0.5
0.4
0.3
0.2
0.2 0.1
0.0 0.0
0.0 1.0
1.5
2.0
2.5
3.0
large strain rock masses. With the decrease in σ0, u0 gradually increases and reaches the maximum value when σ0 = 0. For the hard rock mass corresponding to a high Young’s modulus, u0 slowly increases after excavation. Once the plastic region develops, u0 rapidly increases. However, for the soft rock mass that corresponds to a low Young’s modulus, u0 dramatically increases at the initial unloading stage, and a significant increase also occurs once the plastic region is formed. For example, when σ0 = 0.9p0, u0 reaches 0.569R0 for the soft rock mass with E = 0.1 MPa, and the elastic displacement at the elastic-plastic boundary even reaches 0.9398R0 and 0.9306R0 for the MC and HB cases, respectively. It should be noted that the ultimate displacement is the excavation radius R0, which is in agreement with the actual condition. The influence of Young’s modulus on the displacement u0 and stress state (σcr r, σcr θ) at the elastic-plastic boundary is shown in Fig. 11 for both large strain and small strain conditions. For the hard rock mass, the displacement u0 is approximately the same for the large strain and small strain rock masses; however, the difference gradually increases with the decrease in Young’s modulus. The displacement from the small strain solutions is acceptable before the residual region develops, with a relative error less of than 10% with the large strain solution. However, for the extremely soft rock mass, the difference in u0 between the large and small strain solutions dramatically increases. u0 of the small strain solution is larger than R0 when Young’s modulus is smaller than approximately 25 MPa and 40 MPa, respectively. When Young’s modulus is 10 MPa, the u0 of the small strain rock mass reaches 5.94 and 30.56 times that of the large strain rock mass for the MC and HB rock masses, respectively. Moreover, the small strain displacements u0 = 371.10R0 and 2594.75R0 for the MC and HB rock masses when E = 0.1 MPa, while the large strain displacements are 0.92R0 and 0.98R0. Obviously,
0.6 0.5 0.4 0.3 0.2 0.1
1.2
1.4
1.6
1.8
0.4
Fig. 10. GRCs with different Young’s moduli: (a) MC rock mass; (b) HB rock mass.
0.7
0.0 1.0
0.2
u0/R0
3.5
Rp/R0
σ0/p0
0.6
0.7 0.6
(c)
0.4
u0/R0
u0/R0
(b)
0.2
2.0
2.2
2.4
2.6
Rs/R0 Fig. 9. Evolutions of u0, Rp and Rs under various p0 and σ0 for HB rock mass: (a) u0; (b) Rp and (c) Rs.
and Rs gradually become smaller and smaller. For example, Rp = 1.95R0 and Rs = 1.94R0 for the condition of p0 = 10.0 MPa, which are 62.50% and 93.72% of those of p0 = 1.0 MPa.
6.2. Effect of Young’s modulus The stress equilibrium equation of the large strain theory illustrates that the elastic parameters affect the stress and displacement distributions. A series of Young’s moduli are employed to investigate its influence on the distribution of stress and displacement. Fig. 10 shows the GRCs with respect to different Young’s moduli for both the MC and HB 8
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
1.0
1.2
u0/R0 0.4
0.8 cr
σr /p0
0.2
5
6
8
9
1.0
Rs/R0 of small strain
1.5
0.0 10
1.0
2.0
0.8
Rp/R0 of small strain
2.0
logE
(b)
Rs/R0 of large strain
3.0 2.5
0.4
7
Rp/R0 of large strain
3.5
R/R0
u0/R0
0.6
4.0 cr
Large strain Small strain
5.0 4.5
1.6
cr
σθ /p0
cr
0.8
0.0
(a)
2.0
σr /p0, σθ /p0
(a)
5
6
7
8
9
10
8
9
10
logE
(b)
1.6
15
cr
σθ /p0
0.2
0.4
0.0
0.0 10
5
6
7
8
9
11 9
Rp/R0 of large strain
7
Rp/R0 of small strain
5
Rs/R0 of small strain
R/R0
cr
σr /p0
cr
0.8
cr
u0/R0
u0/R0
0.4
13 1.2
σr /p0, σθ /p0
Large strain Small strain
0.6
logE
3
Fig. 11. Influence of Young’s modulus on the GRCs and critical stress at the elastic-plastic boundary: (a) MC rock mass; (b) HB rock mass.
1
Rs/R0 of large strain
5
6
7
logE u0 is larger than the initial excavation radius for the small strain solutions. However the results u0 > R0 does not make any sense, since the excavation zone has been fully backfilled when u0 = R0. This misleading result is caused by the assumption of small strain theory, which neglects the position changes during the deformation. Additionally, for the small strain rock mass, the stress state at the elastic-plastic boundary is independent of Young’s modulus and remains a constant value, i.e., 0.20 MPa and 0.34 MPa for the MC and HB rock masses. However, for the large strain behavior, the critical radial stress σrcr and tangential stress σθcr gradually increase with the increase in Young’s modulus and approach the limit values, which equal that of the small strain solution. For the soft rock mass with E = 0.1 MPa, σrcr = 5.37e3p0 andσθcr = 0.88p0 for the MC rock mass, while those of the HB rock mass are 0.12p0 and 0.92p0. The influence of Young’s modulus on Rp and Rs is shown in Fig. 12. For the small strain rock mass, with the decrease in Young’s modulus, Rp initially increases with a slow ratio. Once the residual region forms (logE < 7.7 for MC rock mass and logE < 8.0 for HB rock mass), both Rp and Rs significantly increase, and Rs tends to be equal to Rp. Then, both Rp and Rs maintain a relatively constant value with a very small increasing ratio with the decrease in Young’s modulus. For the employed examples, the relative stable softening radius is 4.44R0 and 12.85R0 for the MC and HB small strain rock masses. For the large strain solutions with a high Young’s modulus, Rp also gradually increases with the decrease in Young’s modulus. Once the residual region develops, Rp and Rs significantly increase, and Rp reaches its maximum value when E = 25.12 MPa and 39.81 MPa for the MC and HB rock masses, respectively. Then, Rp decreases with the further decrease in Young’s modulus and gradually becomes R0. The significant increase in plastic radius is induced by the development of residual region. However, when the Young’s modulus is notably small, the elastic deformation would be very large. A larger elastic deformation would leads to a small plastic region, since the position change is considered in the large strain
Fig. 12. Influence of Young's modulus on the radii of plastic and residual regions: (a) MC rock mass; (b) HB rock mass.
solution. A similar evolution can be found for Rs, and the maximum residual radius corresponds to E = 19.95 MPa for both the MC and HB rock mass. 6.3. Effects of Poisson’s ratio Poisson’s ratio also affects the displacement, as shown in Fig. 13. With the increase in Poisson’s ratio, u0, Rp and Rs gradually decrease. For the MC rock mass, when Poisson’s ratio increases from 0.01 to 0.499, u0 decreases from 0.37R0 to 0.21R0, Rp decreases from 2.28R0 to 2.05R0, and Rs decreases from 1.97R0 to 1.59R0, with a change of 43.24%, 10.09% and 19.29%, respectively. For the HB rock mass, u0, Rp and Rs decrease from 0.46R0, 3.13R0, and 2.16R0 to 0.28R0, 2.98R0 and 1.82R0, corresponding to decreases of 39.13%, 4.79% and 15.74%, respectively. 6.4. Effects of the critical plastic shear strain The elastic-perfectly plastic behavior and elastic-brittle-plastic behavior can be considered two extreme cases of strain-softening behavior. By assigning a critical plastic shear strain of γp* to very large and small values, the proposed strain-softening solutions can reduce to the above two behaviors. Fig. 14 shows the influence of γp* on u0, Rp and Rs. For the MC rock mass, u0 = 0.45R0 and Rp = 2.50R0 when γp* = 0, and those of 5.03 × 10−2R0 and 1.17R0 corresponds to the elastic-perfectly plastic behavior with γp* = 2. The displacement and softening regions decrease by 88.82% and 53.20%. When γp* > 0.24, no residual region occurs. Similarly, for the HB rock mass, u0 and Rp decrease from 0.72R0 and 3.66R0 to 0.07R0 and 1.57R0, respectively, corresponding to a decreases of 90.28% and 51.10%. Additionally, the residual region 9
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
0.36
2.3
0.32
2.1
0.28
1.9
(a)
0.5
0.24
u0/R0
1.7
Rs/R0 0.20 0.0
0.1
0.2
0.3
0.4
u0/R0 Rp/R0
0.4
u0/R0 Rp/R0
3.0
0.6
1.5 0.5
Rs/R0
1.5
0.2
1.0
0.1
0.5
0.0 0.0
0.4
0.8
(b)
3.6
0.45
u0/R0
2.8
Rp/R0 Rs/R0
0.35
u0/R0
u0/R0
R/R0
3.2
0.40
1.2
γ
0.50
2.4
1.6
0.25 0.0
2.0
0.1
0.2
0.3
0.4
0.9
4.0
0.8
u0/R0
3.5
0.7
Rp/R0 Rs/R0
3.0
0.6 0.5
2.5 2.0
0.4
1.5
0.2
1.0
0.1
0.5
0.0 0.0
1.6 0.5
0.0 2.0
p*
0.3
0.30
2.0
0.3
v
(b)
2.5
v
R/R0
2.5
R/R0
0.40
R/R0
u0/R0
(a)
0.4
0.8
1.2
γ
1.6
0.0 2.0
p*
Fig. 13. Influence of Poisson’s ratio on the GRCs and critical stress at the elastic-plastic boundary: (a) MC rock mass; (b) HB rock mass.
Fig. 14. Influence of the critical plastic shear strain on displacement, soften radius and residual radius: (a) MC rock mass, (b) HB rock mass.
disappear when γp* > 0.30. Unfortunately, neither of the two extreme conditions is reasonable because the elastic-perfectly plastic behavior overestimates the bearing capacity of the rock mass and the elasticbrittle-plastic behavior contrarily neglects the bearing capacity of the post-failure rock mass. Only a reasonable consideration of post-failure behavior can provide a theoretical foundation for the optimization of structure design and stability evaluation.
decrease, finally approaching the excavation radius. However, the softening and residual radii of the small strain solutions first significantly increase, and the residual radius tends to be equal to the softening radius; later, they maintain a relatively constant value. With the increase in Poisson’s ratio, the displacement, softening radius and residual radius gradually decrease. With the increase in the critical plastic shear strain, the displacement, softening radius and residual radius also gradually decrease, and the mechanical behavior gradually changes from an elastic-brittle-plastic model to an elastic-perfectly plastic model. Thus, the proposed large strain solutions are more generalized conditions for geotechnical engineering.
7. Conclusions Considering the real deformation behavior, a numerical large strain solution of a circular tunnel excavated in MC and HB strain-softening rock masses are derived with a nonassociated MC potential flow rule. An efficient calculation procedure with a five-order Runge-Kutta method is also presented. By simplifying the boundary conditions at the infinite radius with an assumed stress state, the stresses and displacements can be recursively determined. The proposed solutions are validated by the numerical simulation solutions. In addition, the effects on the elastic parameters and critical plastic shear strain are further investigated. The large strain displacement is smaller than that of the small strain solutions, and the displacement at the excavation radius is no larger than the excavation radius. For the cases without a residual region, the difference between the large strain and small strain solutions is very small with a relative error of displacement less than 10%. In addition, the softening radius of the large stain solution is approximately the same as that of the small strain solution. However, for the cases with a residual region, the displacement of the small strain solutions dramatically increases and may exceed the excavation radius when Young’s modulus is low. With the decrease in Young’s modulus, the softening and residual radii of the large strain solution first increase and then
Acknowledgments The authors are grateful to the financial support from the Fundamental Research Funds for the Central Universities (Grant No. 2018ZZCX04). References [1] Brown ET, Bray JW, Landayi B, Hoek E. Ground response curves for rock tunnels. ASCE J Geotech Eng Div 1983;109:15–39. [2] Zhang Q, Jiang BS, Wu XS, Zhang HQ, Han LJ. Elasto-plastic coupling analysis of circular openings in elasto-brittle-plastic rock mass. Theor Appl Fract Mec 2012;60:60–7. [3] Alejano LR, Alonso E. Considerations of the dilatancy angle in rocks and rock masses. Int J Rock Mech Min Sci 2005;42:481–507. [4] Lee YK, Pietruszczak S. A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass. Tunn Undergr Sp Tech 2008;23:588–99. [5] Park KH, Kim YJ. Analytical solution for a circular opening in an elastic- brittleplastic rock. Int J Rock Mech Min Sci 2006;43:616–22. [6] Park KH, Tontavanich B, Lee JG. A simple procedure for ground response curve of circular tunnel in elastic-strain softening rock masses. Tunn Undergr Sp Tech 2008;23:151–9.
10
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
2014;38:1131–48. [18] Vrakas A. A finite strain solution for the elastoplastic ground response curve in tunnelling rocks with non-linear failure envelopes. Int J Numer Anal Met Geomech 2016;41:1077–90. [19] Vrakas A, Anagnostou G. A simple equation for obtaining finite strain solutions from small strain analyses of tunnels with very large convergences. Geotechnique 2015;65:936–44. [20] Park KH. Large strain similarity solution for a spherical or circular opening excavated in elastic-perfectly plastic media. Int J Numer Anal Met Geomech 2015;39:724–37. [21] Zhang Q, Li C, Jiang BS, Yu LY. Elastoplastic analysis of circular openings in elastobrittle-plastic rock mass based on logarithmic strain. Math Probl Eng 2017;7503912:1–9. [22] Guan K, Zhu WC, Wei J, Liu XG, Niu LL, Wang XR. A finite strain numerical procedure for a circular tunnel in strain-softening rock mass with large deformation. Int J Rock Mech Min Sci 2018;112:266–80. [23] Chadwick P. The quasi-static expansion of a spherical cavity in metals and ideal soils. Q J Mech Appl Math 1959;12:52–71. [24] Zhang Q, Jiang BS, Wang SL, Ge XR, Zhang HQ. Elasto-plastic analysis of a circular opening in strain-softening rock mass. Int J Rock Mech Min Sci 2012;50:38–46. [25] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [26] Ogawa T, Lo KY. Effects of dilatancy and yield criteria on displacements around tunnels. Can Geotech J 1987;24:100–13. [27] Varasa F, Alonso E, Alejano LR, Fdez.-Manína G. Study of bifurcation in the problem of unloading a circular excavation in a strain-softening material. Tunn Undergr Sp Tech 2005;20:311–22.
[7] Sharan SK. Elastic-brittle-plastic analysis of circular openings in Hoek-Brown media. Int J Rock Mech Min Sci 2003;40:817–24. [8] Wang SL, Yin XT, Tang H, Ge XR. A new approach for analyzing circular tunnel in strain-softening rock masses. Int J Rock Mech Min Sci 2010;47:170–8. [9] Zhang Q, Zhang CH, Jiang BS, Li N, Wang YC. Elastoplastic coupling solution of circular openings in strain-softening rock mass considering pressure-dependent effect. Int J Geomech, ASCE 2018;18:04017132. [10] Han JX, Li SC, Li SC, Yang WM. A procedure of strain-softening model for elastoplastic analysis of circular opening considering elasto-plastic coupling. Tunn Undergr Sp Tech 2013;37:128–34. [11] Brown ET, Bray JW, Santarelli FJ. Influence of stress-dependent elastic moduli on stresses and strains around axisymmetric boreholes. Rock Mech Rock Eng 1989;22:189–203. [12] Fahimifar A, Ghadami H, Ahmadvand M. An elasto-plastic model for underwater tunnels considering seepage body forces and strain-softening behavior. Eur J Environ Civ En 2015;19:129–51. [13] Zou JF, Qian ZH. Face-stability analysis of tunnels excavated below groundwater considering coupled flow deformation. Int J Geomech, ASCE 2018;18(8):04018089. [14] Lu AZ, Xu GS, Sun F, Sun WQ. Elasto-plastic analysis of a circular tunnel including the effects of the axial in situ stress. Int J Rock Mech Min Sci 2010;47:50–9. [15] Papanastasiou P, Durban D. Elastoplastic analysis of cylindrical cavity problems in geomaterials. Int J Numer Anal Met Geomech 1997;21:133–49. [16] Yu HS, Rowe RK. Plasticity solutions for soil behaviour around contracting cavities and tunnels. Int J Numer Anal Met Geomech 1999;23:1245–79. [17] Vrakas A, Anagnostou G. A finite strain closed-form solution for the elastoplastic ground response curve in tunnelling. Int J Numer Anal Met Geomech
11
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
A numerical large strain solution for circular tunnels excavated in strainsoftening rock masses Qiang Zhanga, Hong-Ying Wangb,c, Yu-Jing Jiangb, Meng-Meng Lua, Bin-Song Jianga,
T
⁎
a School of Mechanics & Civil Engineering, State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining & Technology, Xuzhou 221116, China b Geoenvironmental Laboratory, Graduate School of Engineering, Nagasaki University, Nagasaki 8528521, Japan c College of Energy and Transportation Engineering, Jiangsu Vocational Institute of Architectural Technology, Xuzhou 221116, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Circular tunnel Large strain Strain-softening Yield criterion Elastic-plastic analysis
This paper presents a large strain numerical solution of a circular tunnel for the ground reaction curve (GRC) in strain-softening rock masses. Both the stress equilibrium equation and deformation compatibility equation are considered in a current coordinate. A series of examples are employed to validate the proposed solutions, and the results show that: (1) the displacement, softening and residual radii of larger strain solution are smaller than that of small strain solution, and their difference significantly increases with a decrease in the Young’s modulus. (2) the increasing Poisson’s ration and critical shear plastic strain decrease the displacement, softening and residual radii.
1. Introduction The tunnel clearance size and supporting structure parameter are of essential importance in the geotechnical engineering design. A reasonable prediction of the displacement and failure region for the surrounding rock around the tunnel would be the primary influence on the design optimization and economics. In recent years, many studies on elastic-plastic rock masses have been carried out, and many achievements have been obtained [1–10]. Among these studies, scholars mainly pay attention to the Ground Response Curves (GRCs) according to different constitutive models, i.e., the elastic-perfectly (brittle) plastic model [1,2] and strain-softening model [3–8], with numerical and (semi-)analytical methods that follow small strain theory. To improve these models, the deterioration effects of the elastic parameters [9,10], confining pressure [11], seepage [12,13] and out-of-plane stress [14] are further considered. Strictly speaking, those solutions are only appropriate for engineering applications, in which the rocks/structures deform negligibly small, such as a tunnel excavated in the hard rock by assuming that small strain theory is suitable. Unfortunately, the convergence displacement may exceed one-tenth of the initial excavation radius in the surrounding soft strata subjected to high in-situ stress with a significant decrease in the cavity radius. For example, the displacement exceeds 10 ∼ 14% of the initial excavation radius in the Wushaoling tunnel, and that of the Mingyazi tunnel even exceeds 30% of the initial excavation radius, as shown in Fig. 1. The
⁎
displacement predicted using small strain theory deviates significantly from the natural cases because the position of material point is assumed to be fixed, which however has a non-negligible deformation in such large stain cases. Since the supporting pressure is always acting on the deformed excavation radius, considering the stress equilibrium on the current configuration instead of on the initial configuration seems more suitable for large deformation problems. Therefore, the stress and deformation compatibility equations should be always satisfied for the deformed configuration. Considering the mechanical behavior of large strain rock masses, a few attempts have been made to derive the solutions of circular tunnels. Papanastasiou and Durban [15] studied the expansion and contraction behavior of cylindrical cavities in an elastic-perfectly plastic rock mass and derived the differential solutions. Later, Yu and Rowe [16] presented an analytical solution for the problem of a cavity. With the assumption of small strain in the elastic region and considering the out-ofplane stress effect, Vrakas and Anagnostou [17–19] derived the closedform solution of a circular tunnel using the linear Mohr-Coulomb criterion and non-linear envelope criterion. Additionally, a closed-form relationship between the displacement of small and large strain solutions for the cavity contraction/expansion problem was presented considering the general Mohr’s criterion and constant plastic dilatancy. By assuming the elastic rock mass as small stain condition, Park et al. [20] derived the large strain similarity solution for a spherical or circular opening excavated in elastic-perfectly plastic media, and the
Corresponding author. E-mail addresses: [email protected] (Q. Zhang), [email protected] (Y.-J. Jiang), [email protected] (B.-S. Jiang).
https://doi.org/10.1016/j.compgeo.2019.103142 Received 26 July 2018; Received in revised form 2 May 2019; Accepted 17 June 2019 0266-352X/ © 2019 Published by Elsevier Ltd.
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
(a)
(a)
ı
Design contour
Elastic-perfectly plastic u >100cm about 10~14%R
Strain-softening Deformed contour
Elastic-brittle-plastic
İ (b) Elastic-perfectly plastic
-İ3
(b)
Design contour
Strain-softening
u >230cm over 30%R
Elastic-brittle-plastic
İ1 Deformed contour
Fig. 2. Mechanical behavior of rock mass: (a) stress-strain curves, and (b) strain curves.
Fig. 1. Large deformation examples of tunnel contraction: (a) Wushaoling tunnel; (b) Mingyazi tunnel.
Strength drop
Load /kN
comparisons between the small and large strain solutions denoted that the displacement assuming large strain is smaller than that assuming small strain. Considering the large strain behavior of both elastic and plastic rock masses, Zhang et al. [21] derived the differential expressions for a circular opening excavated in elastic-brittle-plastic rock. The examples denoted that the stresses and displacement of the elastic rock mass depend on the elastic parameters, and the difference between the large strain solution and small strain solution may be significant, especially for extremely soft surrounding rock. Guan et al. [22] proposed a finite strain numerical procedure of a circular tunnel in strainsoftening rock mass considering the axial stress effect. Similar to the former study of Vrakas and Anagnostou, the elastic rock mass is assumed following the small strain behavior. Clearly, most of the studies on large strain problems are carried out according to elastic-perfectly (brittle) plastic behavior. For the elasticbrittle-plastic model, a sudden drop in stress immediately happens once the peak stress is reached. Meanwhile, an ideal plastic flow performs when the yield criteria is satisfied for the elastic-perfectly plastic model. The elastic-brittle-plastic model and elastic-perfectly plastic model are particular conditions of mechanical behaviors of rock mass. They are only reasonable for the hard and soft rock mass, respectively. Most of the geotechnical materials generally manifest strain-softening behavior in the post-failure region, whose mechanical behavior is between the above two special cases, as shown in Fig. 2. Fig. 3 shows a typical cyclic loading deformation curve of a rock sample. With the increase in
Envelope curve Test curve
Axial deformation /mm Fig. 3. Typical cyclic loading deformation curve of a rock sample.
deformation, the bearing capacity of the post-failure rock mass gradually decreases. The progressive failure process is very important for stability evaluation of the surrounding rock and optimization of the lining structure in underground engineering. To provide a more reasonable prediction in the GRC and failure region for underground geotechnical engineering, this paper aims at the numerical large strain solutions of circular tunnels excavated in strainsoftening rock masses. The differential expressions are derived for both elastic and plastic rock masses, in which the material position is
2
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
p0
in the Eulerian form as follows:
(
Softening-residual boundary
ıcr r
(
u0
p0
ı0 r0 R0
Rp
⎧ εr = ⎨ εθ = ⎩
u 2 r
(3)
1 [(1 2G 1 [(1 2G
− v )(σr − p0 ) − v (σθ − p0 )] + εrp − v )(σθ − p0 ) − v (σr − p0 )] + εθp
(4)
where G = E/[2(1 + v)] is the shear modulus; E and v are Young’s modulus and Poisson’s ratio, respectively; εrp and εθp are the radial and tangential plastic strains; and εrp ≡ εθp ≡ 0 for the elastic rock mass. With respect to the strain-softening behavior, the softening parameter, which governs the evolution of the material properties, should be determined first. Here, the most widely used plastic shear strain γ p is employed as the softening index, as shown in Eq. (4).
Elastic-softening boundary
p0 Fig. 4. Circular tunnel in an infinite deformable medium.
γ p = εθp − εrp
(5)
The Mohr-Coulomb (MC) plastic potential function is considered (i.e., non-associated plastic flow). According to the plastic potential theory [24], the radial and tangential plastic strains would relate to
deformable. Then, a recursive implementation procedure is proposed, and the proposed solutions are validated by comparisons with numerical solutions. Finally, parametric analyses of Young’s modulus, Poisson’s ratio and the critical plastic shear strain are investigated.
εrp + κ (γ p) εθp = 0
(6)
1 + sin ψ (γ p) , 1 − sin ψ (γ p)
where κ (γ p) = and ψ is the dilation angle. The yield criterion for the plastic rock mass can be expressed as
2. Problem description 2.1. Geometrical relationship
σθ − σr = H (σr , γ p)
Fig. 4 shows a circular tunnel with an initial radius R0 excavated in an infinite isotropic rock mass, and a hydrostatic stress field p0 is imposed throughout the domain prior to the excavation. The assumption of a hydrostatic stress field implies that the body force is neglected, and no rotation of the principal axes occurs for this axisymmetric problem. Then, the problem transforms into a one-dimensional problem in terms of a radial displacement. Additionally, contraction displacement is defined as positive. With the decrease in the supporting pressure σ0 from its initial value p0, the excavation radius gradually converges to r0 with a displacement of u0. For the deformable surrounding rock, the initial position of a material point R is now at the current radius r. Therefore, the radial displacement u of the material point can be defined as
u=R−r
du 2 dr
( ) +⋯ ( ) +⋯
where εr and εθ denote the radial and tangential strains, respectively. In a small strain condition, the higher order terms on the right side of Eq. (3) are neglected, leading to the well-known kinematic equations. The total strain induced by excavation is composed of elastic strain and plastic strain. Based on Hooke’s law, the constitutive equation can be formulated as
Rs
p0
) )
du du 1 ⎧ ⎪ εr = ln 1 + dr = dr − 2 ⎨ ε = ln 1 + u = u − 1 ⎪ θ r r 2 ⎩
(7)
As introduced above, the material parameters vary with the softening parameter. For the linear MC yield criterion, H(σr, γp) can be written as
2c (γ p) cos φ (γ p) 1 + sin φ (γ p) HMC (σr , γ p) = ⎜⎛ − 1⎟⎞ σr + p 1 − sin φ (γ p) ⎠ ⎝ 1 − sin φ (γ )
(8)
where φ is the friction angle and c is the cohesive strength, respectively. For the generalized nonlinear Hoek-Brown (HB) yield criterion, H (σr, γp) becomes
HHB (σr ,
γ p)
= σc
(γ p) ⎛⎜m (γ p) ⎝
(1)
σr + s (γ p) ⎞⎟ σc (γ p) ⎠
a (γ p)
(9)
where σc is the uniaxial strength; m, s and α are three strength parameters, which are related to the geological strength index.
Once the internal supporting pressure σ0 becomes lower than a critical value σrcr , a current plastic region with a radius of Rp develops around the tunnel. In some cases, the plastic region of rock includes another current residual region with a radius Rs besides Rp, in which the material properties reach its residual value.
2.3. Evolution of material parameters
For the axisymmetric plane strain problem, the stress equilibrium equation in the current state can be expressed as follows:
The evolution of the material parameters mainly depends on the lithology and stress state, and it can be obtained through cyclic loadingunloading triaxial compression tests under various confining pressures. For the sake of comparison, a piecewise linear function of γp is assumed in this paper, as shown in Eq. (10).
dσr σ − σθ + r =0 dr r
η (γ p) =
2.2. Governing equations
γp
(2)
where σr and σθ denote the radial stress and tangential stress, respectively. For the present problem, the logarithmic strains can accurately describe the strain-displacement relationship [23]. In the current polar coordinates, the logarithmic radial and tangential strains can be written
⎧ ηp − (ηp − ηr ) γ p ∗ , ⎨ η, r ⎩
0 < γ p < γ p∗ γ p ⩾ γ p∗
(10)
where η denotes any one of the strength parameters, i.e., φ, c, σc, m, s, α and ψ; γp* is the critical plastic shear strain for the residual state of the surrounding rock; and ηp and ηr are the peak and residual values of η, respectively. 3
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
The plastic rock mass should also satisfy the deformation compatibility equation. By using Eqs. (1), (2), (4), (7) and (17), the deformation compatibility equation of a plastic rock mass can be formulated in the original coordinate as follows:
2.4. Boundary conditions Since the supporting pressure is always acting on the excavation boundary, we have σr(r0) = σ0. The stress state at an infinite distance from the tunnel remains undisturbed; thus, it leads to lim σr (r ) = lim σθ (r ) = p0 and lim u (r ) = 0 . The current position of a
r →∞
r →∞
dεθp
r →∞
dR
material point is still unknown before the exact solution is obtained. However, Eq. (1) implies that the material points’ positions are a oneto-one correspondence in the original and deformed coordinates. Then, the boundary condition can also be expressed in the original coordinate as follows:
⎧ R = R0, ⎨ ⎩ R → ∞,
σr = σ0 σr = σθ = p0
R
-
ζH (σr , γ p) H (σr , γ p)/2G + (1 + κ (γ p)) ε p θ e 2GR
1 + sin φ (γ p) ∂HMC (σr , γ p) −1 = 1 − sin φ (γ p) ∂σr
(19)
(11)
∂HHB (σr , γ p) σr = m (γ p) α (γ p) ⎜⎛m (γ p) + s (γ p) ⎟⎞ ∂σr σ ( γ p) c ⎝ ⎠
α (γ p) − 1
(20)
The tangential stress can be easily obtained by substituting the radial stress into the yield criterion. Using Eqs. (1) and (2), the displacement of a plastic rock mass can be obtained as follows: p
u = R {1 − e−[(1 − v)(σθ− p0 ) − v (σr − p0 )]/2G − εθ }
Theoretically, the material parameters at each point should obey the evolution of Eq. (10), and this is the main challenge in the strain-softening analysis. In the finite differential method (FDM), only the integration points exactly satisfy the evolution of the material parameters. The stresses at the other points are estimated by the shape functions and the stresses of the integration points [25]. Therefore, the material parameter for each element is represented by a series of values at the integration points. For the one-dimensional problem of a circular tunnel, the surrounding rock is divided into n annuli, and each annulus is assumed to be isotropic with a constant material parameter, which depends on the plastic shear strain at its outer radius. So far, the stresses and displacements in elastic and plastic surrounding rocks are obtained in a differential equation form. To solve the differential equations of Eqs. (14)–(16) and Eqs. (18)–(21) for the elastic and plastic rock masses, many ordinary differential methods are available, such as the Euler and Runge-Kutta methods. Unlike in the small strain solutions, the stresses depend on the elastic parameters, i.e., Young’s modulus and Poisson’s ratio, suggesting that the explicit stress at the elastic-softening boundary cannot be obtained. In this condition, a numerical method is necessary to solve the differential equations for the elastic and plastic rock masses. To minimize the cumulative error, a fifth-order RungeKutta method is employed, as shown in Eq. (22):
(12)
(13)
By substituting the elastic constitutive equation of Eq. (4) in Eqs. (12) and (13), the stress and deformation compatibility equations of the elastic surrounding rock can be expressed in terms of the radial and tangential stresses as follows:
dσr σ − σr (σθ− σr )/2G e = θ dR R
(14)
dσθ 2G v (σθ − σr ) (σθ− σr )/2G [1 − e (σθ− σr )/2G] + e = dR (1 − v ) R (1 − v ) R
(15)
(21)
4. Numerical implementation procedure
The radial and tangential strains are functions of the displacement with respect to its current position. Therefore, the strains should satisfy the deformation compatibility equation, which can be obtained by substituting the tangential strain into the radial strain in Eq. (3). Then, using the relationship between the original and current coordinates of Eq. (1), the deformation compatibility equation can be expressed in the original coordinate as follows:
35 500 125 2187 11 yi + 1 = yi + ⎛ k1 + k3 + k4 − k5 + k6 ⎞ h 1113 192 6784 84 ⎠ ⎝ 384
Once the stresses are obtained, the displacement of the elastic surrounding rock can be calculated with the tangential geometrical equation given in Eq. (3). Then, the displacement can be expressed in the original coordinate by using Eq. (1), as shown in Eq. (16).
(22)
where h is the step length and
k1 = f (x i , yi ) =
(16)
3.2. Solution in the plastic region By using the same solution approach as that used for the elastic rock mass, the stress equilibrium equation of a plastic rock mass can be reformulated in the original coordinate using Eqs. (1), (3) and (4). Meanwhile, the stress state should satisfy the yield criterion of Eq. (7). By using the yield criterion, the stress equilibrium equation can be formulated as follows:
H (σr , γ p) H (σr , γ p)/2G + (1 + κ (γ p)) ε p dσr θ = 0 e = dR R
)/2G + (1 + κ (γ p)) εθp
(18)
The stress equilibrium equation is expressed in the current (deformed) position, as shown in Eq. (2), but the displacement of the material point is unknown. To solve the present problem, the stress equilibrium equation should be expressed in the original coordinate. Using Eqs. (1) and (3), the stress equilibrium equation can be rewritten in the original coordinate as follows:
u = R {1 − e−[(1 − v)(σθ− p0 ) − v (σr − p0 )]/2G}
p
∂H (σ , γ p)
3.1. Solution in the elastic region
dεθ 1 = (1 − e εθ− εr ) dR R
1 − e H (σr , γ
− 2v + 1; for the MC and HB rock masses, the where ζ = (1 − v ) ∂σr r derivate term can be respectively expressed as follows:
3. Solutions for stress and displacement
dσr σ − σθ εθ− εr e + r =0 dR R
=
d (x i , yi ) dx
(23a)
1 1 k2 = f ⎛x i + h, yi + k1 h⎞ 5 5 ⎝ ⎠
(23b)
3 3 9 k3 = f ⎛x i + h, yi + k1 h + k2 h⎞ 10 40 40 ⎠ ⎝
(23c)
4 44 56 32 k 4 = f ⎛x i + h, yi + k1 h − k2 h + k3 h⎞ 5 45 15 9 ⎝ ⎠
(23d)
8 19372 25360 64448 212 k5 = f ⎛x i + h, yi + k1 h − k2 h + k3 h − k 4 h⎞ 9 6561 2187 6561 729 ⎝ ⎠ (23e)
(17) 4
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
k6 = f ⎛x i + h, ⎝ 9017 355 46732 49 5103 yi + k1 h − k2 h + k3 h + k4 h − k5 h⎞ 3168 33 5247 176 18656 ⎠
annulus at Ri is employed to predict the material properties of the ith annulus with Eq. (10). The tangential plastic strain εθp(−i) is employed as the initial value for the ith annulus for solving the differential equation of Eq. (18). Only in this way can the displacement be continuous; otherwise, the accumulative error will lead to inaccurate results, especially for the elastic-brittle-plastic rock mass. Similar to the calculation method for the elastic-softening radius, the plastic shear strain equals γp* at the softening-residual boundary. For the plastic rock mass, the calculated plastic shear strain is compared with γp* in each step. Once the calculated plastic shear strain is larger than γp*, the inner radius is adjusted to ensure that the plastic shear strain at the inner radius equals γp*. Then, the stress and displacement of each annulus can be recursively calculated until the radius R reaches R0. If the required supporting pressure σ¯r (1) at R = R0 corresponding to the given tangential σ¯θ (n) equals the actual supporting pressure σ0, the calculated results would be corrected. Otherwise, the initial tangential stress σ¯θ (n) is updated until σ¯r (1) = σ0 . Additionally, an adaptive step length method can be used to improve the computation efficiency. For convenience, in this study, a fixed small step length is considered throughout the entire solution to ensure the stability of the solution. The detailed calculation flowchart of the implementation is shown in Fig. 5.
(23f)
To solve the differential equations of the ith annulus, the initial slopes of each variable at the boundary xi should be given first. However, Eq. (11) denotes that the exact boundary condition cannot be accurately given due to the infinite radius. According to Saint-Venant's principle, the deviation would be very small when the boundary condition is replaced by a very large finite radius, i.e., 50R0. Hence, the boundary conditions can be replaced by
⎧ R = R0, ⎨ ⎩ R = 50R 0 ,
σr = σ0 σr = p0
(24)
Then, the present problem transforms into a boundary problem. In the contraction problem of a circular tunnel, the tangential stress is larger than the initial stress in the elastic region because of stress concentration. For the plane strain problem of circular tunnel, the elastic tangential stress is the first principal stress, and the radial stress is the minimum principal stress, that is σθ > σr according to the Lame’s solution. Thus, an initial tangential stress σ¯θ (n) > p0 at R = 50R0 is assumed, and the slope k1 can be obtained for the radial stress and tangential stress. The slope k1 is used to predict the dependent variables at xi + h/5; at this point, the slope k2 can also be determined. Other slopes, i.e., k3−k6, which are in the region (xi, xi+1), can be continuously obtained. Slopes k1−k6 are combined to predict the variables at xi+1 with Eq. (22). Therefore, all of the stresses and displacements for the elastic region can be recursively determined. In the plastic region, the variables are the radial stress and tangential plastic strain, and the radius Rp is still unknown. At the elasticsoftening boundary, the elastic stresses should satisfy the peak yield criterion due to zero plastic strain. For each annulus, the calculated elastic stresses are verified by the peak yield criterion. If the calculated stress (σr(i), σθ(i)) satisfies the peak yield criterion, the elastic-softening boundary is located in the range of (xi, xi+1). The step length is adjusted by dichotomy until the stress (σr(i), σθ(i)) is located on the peak yield surface. Then, the differential equations given in Eqs. (17), (7) and (18) are used to calculate the stresses and tangential plastic strain for the plastic rock mass. With the aid of Eq. (6), the radial plastic strain can be obtained. By substituting the radial and tangential plastic strains into Eq. (5), the plastic shear strain at the inner radius can be obtained. It should be noted that the radial stress and displacement should be continuous, and the stress state of each annulus satisfies the yield criterion. However, the material properties of each annulus are different for the softening rock mass, which may cause a sudden change in the tangential plastic strain at the boundaries of the adjacent annuli. For the (i + 1)th and ith annuli in the softening region, by using the displacement expression of Eq. (21), the tangential plastic strain at R = Ri for the (i + 1)th and ith annuli will be related by
εθp(−i) = εθp(+i) +
1 [(1 2G(i + 1)
−
1 [(1 2G(i)
5. Example verifications 5.1. Verifications for an elastic-perfectly (brittle) plastic rock mass To validate the proposed solutions, the typical input data for the MC and HB rock masses used by Ogawa and Lo [26] are utilized and listed in Table 1. The solutions are validated with the methods proposed by Zhang et al. [21] for both the MC and HB rock masses. For accuracy, a small step length of 0.01R0 and 0.001R0 for the elastic and plastic rock masses are considered, respectively, and an acceptable tolerance of 0.001p0 is employed for the supporting pressure. When a significantly large or small critical plastic shear strain, i.e., γp* = 10 or 1e-5, is employed, the proposed solutions would transform into the solutions of elastic-perfectly plastic and elastic-brittle-plastic rock masses, respectively. The dimensionless plastic radius Rp/R0 and dimensionless displacement u0E/p0R0 are calculated, as listed in Table 2. The calculated results are compared with those using the method of Zhang et al., and the comparison results denote that the proposed solutions are in agreement, with a relative error less than 0.6% for both the dimensionless softening radius and displacement. 5.2. Verifications for a strain-softening rock mass To investigate the strain-softening behavior of a rock mass with large strains, two sets of input data are presented for the MC and HB rock masses, as visualized in Fig. 6 and Fig. 7, respectively. Several values of step length are employed to investigate the convergence of the proposed solutions. For the MC rock mass, the dimensionless displacement u0/R0 is 0.3195 when the step lengths are 0.01R0 and 0.001R0 for the elastic and plastic regions, respectively, and that of 0.3222 corresponds to the step lengths of 1e-4R0 and 1e-5R0. The relative error between the two conditions is less than 8.5‰. Therefore, the results can be considered the actual solutions when the step lengths of 0.01R0 and 0.001R0 are employed for the elastic and plastic rock masses, respectively, which are employed in the following analyses. Fig. 6 shows the GRC and evolutions of Rp and Rs for the MC rock mass. The calculation results are compared with the FDM solutions using the Fast Lagrangian Analysis of Continua (FLAC) software. Since the calculation accuracy of FDM solutions depends on the mesh generation, a minimum element size of 0.1R0 in the concerned region of R < 5R0 is used in the FLAC. Although smaller element sizes result in more accurate solutions, it also causes a localization effect when the element size is small enough, and the softening and residual region would be a series of belt regions [27]. Therefore, an average value of
− v(i + 1) ) H (σr (i), γ(pi + 1) ) + (1 − 2v(i + 1) )(σr (i) − p0 )] − v(i) ) H (σr (i), γ(pi) ) + (1 − 2v(i) )(σr (i) − p0 )] (25)
εθp(−i)
εθp(+i)
and are the tangential plastic strains at Ri for the ith and where (i + 1)th annuli, respectively. Similarly, the plastic shear strains at Ri for the (i + 1)th and ith annuli are also different. Strictly speaking, taking the material parameters of the middle point as the representative values of each annulus is more reasonable. However, this will cause more complex expressions, which should be solved by an iterative algorithm. Zhang et al. [24] studied the influence of four different representative points on the calculation results for the small strain-softening problem, and the four results tend to be uniform as the annuli number increases. To reduce the difficulty of the calculation, the plastic shear strain of the (i + 1)th 5
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
Table 2 Computed results of the dimensionless radius and displacement (data in parentheses are the solutions by using the methods of [21]). Rock mass
Calculation scheme
Rp/R0
MC rock
ebp ep ebp ep
1.7487 1.1642 1.2362 1.6110
HB rock
u0E/(p0R0) (1.7576) (1.1649) (1.2364) (1.6123)
11.9279 (11.9596) 1.5583 (1.5666) 1.1861 (1.1883) 5.9844 (6.0176)
Note: ‘ebp’ denotes the elastic-brittle-plastic behavior, and ‘epp’ denotes the elastic-perfectly plastic behavior.
(a)
0.25
This paper FDM
0.20 R0=1.0, p0=1MPa, E=30MPa, v=0.3
σ0/p0
0.15
p*
γ =0.15, cp=0.2MPa, cr=0.02MPa ϕp=40°, ϕr=20°, ψp=10°, ψr=5°
0.10 0.05 0.00 0.0
0.1
0.2
0.3
0.4
u0/R0
(b)
0.24 0.21 0.18 Rp/R0
σ0/p0
0.15
Rs/R0
0.12 0.09 0.06 0.03 0.00 1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.4
2.6
R/R0 Fig. 5. Detailed flow chart for the sequence of calculations.
Fig. 6. Large deformation solution for the MC rock mass (a) GRC; (b) evolution of Rs and Rp.
Table 1 Geometry and material property parameters. MC rock Radius of opening, R0 (m) Initial stress, p0 (MPa) Supporting pressure, σ0 (MPa) Young’s modulus, E (GPa) Poisson ratio, v(-) Dilation angle, ψp = ψr (Deg) cp (MPa) φp (Deg) cr (MPa) φr (Deg)
the displacement at the excavation surface is taken as the displacement at R0, and the evolutions of Rp and Rs are not listed due to the bifurcation. With the decrease in the supporting pressure, the displacement, softening radius and residual radius gradually increase and reach maximum values when the supporting pressure is completely released. The GRC and evolutions of Rp and Rs for the large strain rock mass are similar to those for the small strain rock mass. The dimensionless displacement u0/R0 is 0.319 and 0.323 for the proposed solution and FDM solution, respectively, with a relative error of 4.4%. The maximum softening and residual radii are 2.25R0 and 1.87R0, respectively. Clearly, the displacement of the proposed total strain solution is slightly smaller than that of the FDM solutions with an incremental velocity approach. This is in agreement with Park’s studies on large strain solutions for an elastic-perfectly plastic rock mass [20]. For the generalized HB rock mass, a composite flow rule with an HB
HB rock 1 1 0 5 0.2 30 0.276 35 0.055 30
Radius of opening, R0 (m) Initial stress, p0 (MPa) Supporting pressure, σ0 (MPa) Young’s modulus, E (GPa) Poisson ratio, v(-) Dilation angle, ψp = ψr (Deg) Uniaxial strength, σc (MPa) mp sp mr sr
1 1 0 5 0.2 30 50 0.2 1e-4 0.05 1e-5
6
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
(a)
0.40 0.35
Proposed solution FDM
0.30
R0=1.0 m, p0=1 MPa, E=40 MPa, v=0.3
σ0/p0
0.25
σp=10 MPa, σr=8 MPa, ap=0.5, ar=0.55
0.20
σ0/p0
(a)
mp=0.5, mr=0.05, sp=1e-4, sr=1e-5
0.15
p*
ψp=5°, ψr=3°, γ =0.15
p0 = 1.0 MPa
0.25
p0 = 2.0 MPa
0.20
p0 = 5.0 MPa
0.05
0.05
0.1
0.2
0.3
0.4
p0 = 1.5 MPa p0 = 3.0 MPa p0 = 10.0 MPa
0.15 0.10
0.00 0.0
0.5
u0/R0
0.2
0.4
0.6
0.8
1.0
u0/R0
(b)
0.40 0.35
0.35 0.30
0.30
0.25
Rp/R0 Rs/R0
0.25
0.20
0.20
σ0/p0
σ0/p0
0.30
0.10
0.00 0.0
(b)
0.35
0.15
0.15 0.10
0.10
0.05
0.05 0.00 1.0
1.5
2.0
2.5
3.0
0.00 1.0
3.5
1.2
1.4
1.6
1.8
2.0
2.2
2.4
2.6
Rp/R0
R/R0
(c)
Fig. 7. Large deformation solution for the HB rock mass (a) GRC; (b) evolution of Rs and Rp.
0.35 0.30
potential function is employed in the FLAC software. To compare these results with the proposed solutions, a new constitutive model based on the MC potential function, in which the radial and tangential plastic strains follow Eq. (6), is developed. When the supporting pressure is completely released, the average dimensionless displacement u0/ R0 = 0.432, and the proposed solution results in u0/R0 = 0.399, Rp/ R0 = 3.122 and Rs/R0 = 2.067, whose displacement is also slightly smaller than that of the FDM solutions with a percentage error of 8.3%, as shown in Fig. 7. In summary, the proposed large strain solutions are in agreement with that of the FDM solutions for both the MC and HB strain-softening rock mass.
0.25
σ0/p0
0.20 0.15 0.10 0.05 0.00 1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.4
Rs/R0
6. Parametric analyses
Fig. 8. Evolutions of u0, Rp and Rs under various p0 and σ0 for MC rock mass: (a) u0; (b) Rp and (c) Rs.
In this section, the influences of supporting pressure and hydrostatic stress on u0, Rp and Rs are studied, and the material property parameters for MC and HB rock mass are the same as that used in Section 5.2 as listed in Figs. 6 and 7, respectively. In this section, those material parameters are considered.
with a decreasing rate and gradually approaches to the limit of R0. However, Rp and Rs do not monotonously increase with the increasing p0. Taking the condition of σ0 = 0 as an example, Rp = 2.22 R0 and Rs = 2.20R0 corresponding to p0 = 2.0 MPa, which are smaller than those of 2.39R0 and 2.22R0 corresponding to p0 = 1.5 MPa. With a further increase in p0, Rp and Rs still decrease. Thus, when the position changes in material point is considered, a very large ground stress may induce a small plastic and residual region. A similar evolutions in u0, Rp and Rs under various σ0 and p0 could be found for HB rock mass as shown in Fig. 9. When p0 > 2.0 MPa, an evident inflection point occurs, and from then on the increasing rate of u0, Rp and Rs become much small. Additionally, with the further increase in p0 (p0 > 2.0 MPa), Rp
6.1. Effects of supporting pressure and hydrostatic stress The evolutions of u0, Rp and Rs for MC and HB rock masses with various of supporting pressure and hydrostatic stress are shown in Figs. 8 and 9, respectively. For MC rock mass, u0, Rp and Rs gradually increase with the decrease in σ0 for a constant p0. Under the condition of p0 < 1.5 MPa, Rp and Rs significantly increase when σ0 < 0.1p0. However, the increasing rate for σ0 < 0.1p0 becomes much small when p0 > 1.5 MPa. On the other hand, with the increase in p0, u0 increases 7
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
(a)
(a)
0.7
1.0
p0 = 1.0 MPa p0 = 1.5 MPa
0.6
0.8
p0 = 2.0 MPa
0.5
p0 = 3.0 MPa p0 = 5.0 MPa
0.6
p0 = 10.0 MPa
σ0/p0
σ0/p0
0.4 0.3 0.2
E=0.1 MPa E=0.3 MPa E=3.0 MPa E=30 MPa E=3.0 GPa
0.4 0.2
0.1 0.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.0
1.2
(b)
1.0
0.6
0.8
1.0
1.0
E=0.1 MPa E=0.4 MPa E=4.0 MPa E=40 MPa E=4.0 GPa
σ 0/ p 0
0.6
0.4
σ0/p0
0.8
0.8
0.5
0.4
0.3
0.2
0.2 0.1
0.0 0.0
0.0 1.0
1.5
2.0
2.5
3.0
large strain rock masses. With the decrease in σ0, u0 gradually increases and reaches the maximum value when σ0 = 0. For the hard rock mass corresponding to a high Young’s modulus, u0 slowly increases after excavation. Once the plastic region develops, u0 rapidly increases. However, for the soft rock mass that corresponds to a low Young’s modulus, u0 dramatically increases at the initial unloading stage, and a significant increase also occurs once the plastic region is formed. For example, when σ0 = 0.9p0, u0 reaches 0.569R0 for the soft rock mass with E = 0.1 MPa, and the elastic displacement at the elastic-plastic boundary even reaches 0.9398R0 and 0.9306R0 for the MC and HB cases, respectively. It should be noted that the ultimate displacement is the excavation radius R0, which is in agreement with the actual condition. The influence of Young’s modulus on the displacement u0 and stress state (σcr r, σcr θ) at the elastic-plastic boundary is shown in Fig. 11 for both large strain and small strain conditions. For the hard rock mass, the displacement u0 is approximately the same for the large strain and small strain rock masses; however, the difference gradually increases with the decrease in Young’s modulus. The displacement from the small strain solutions is acceptable before the residual region develops, with a relative error less of than 10% with the large strain solution. However, for the extremely soft rock mass, the difference in u0 between the large and small strain solutions dramatically increases. u0 of the small strain solution is larger than R0 when Young’s modulus is smaller than approximately 25 MPa and 40 MPa, respectively. When Young’s modulus is 10 MPa, the u0 of the small strain rock mass reaches 5.94 and 30.56 times that of the large strain rock mass for the MC and HB rock masses, respectively. Moreover, the small strain displacements u0 = 371.10R0 and 2594.75R0 for the MC and HB rock masses when E = 0.1 MPa, while the large strain displacements are 0.92R0 and 0.98R0. Obviously,
0.6 0.5 0.4 0.3 0.2 0.1
1.2
1.4
1.6
1.8
0.4
Fig. 10. GRCs with different Young’s moduli: (a) MC rock mass; (b) HB rock mass.
0.7
0.0 1.0
0.2
u0/R0
3.5
Rp/R0
σ0/p0
0.6
0.7 0.6
(c)
0.4
u0/R0
u0/R0
(b)
0.2
2.0
2.2
2.4
2.6
Rs/R0 Fig. 9. Evolutions of u0, Rp and Rs under various p0 and σ0 for HB rock mass: (a) u0; (b) Rp and (c) Rs.
and Rs gradually become smaller and smaller. For example, Rp = 1.95R0 and Rs = 1.94R0 for the condition of p0 = 10.0 MPa, which are 62.50% and 93.72% of those of p0 = 1.0 MPa.
6.2. Effect of Young’s modulus The stress equilibrium equation of the large strain theory illustrates that the elastic parameters affect the stress and displacement distributions. A series of Young’s moduli are employed to investigate its influence on the distribution of stress and displacement. Fig. 10 shows the GRCs with respect to different Young’s moduli for both the MC and HB 8
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
1.0
1.2
u0/R0 0.4
0.8 cr
σr /p0
0.2
5
6
8
9
1.0
Rs/R0 of small strain
1.5
0.0 10
1.0
2.0
0.8
Rp/R0 of small strain
2.0
logE
(b)
Rs/R0 of large strain
3.0 2.5
0.4
7
Rp/R0 of large strain
3.5
R/R0
u0/R0
0.6
4.0 cr
Large strain Small strain
5.0 4.5
1.6
cr
σθ /p0
cr
0.8
0.0
(a)
2.0
σr /p0, σθ /p0
(a)
5
6
7
8
9
10
8
9
10
logE
(b)
1.6
15
cr
σθ /p0
0.2
0.4
0.0
0.0 10
5
6
7
8
9
11 9
Rp/R0 of large strain
7
Rp/R0 of small strain
5
Rs/R0 of small strain
R/R0
cr
σr /p0
cr
0.8
cr
u0/R0
u0/R0
0.4
13 1.2
σr /p0, σθ /p0
Large strain Small strain
0.6
logE
3
Fig. 11. Influence of Young’s modulus on the GRCs and critical stress at the elastic-plastic boundary: (a) MC rock mass; (b) HB rock mass.
1
Rs/R0 of large strain
5
6
7
logE u0 is larger than the initial excavation radius for the small strain solutions. However the results u0 > R0 does not make any sense, since the excavation zone has been fully backfilled when u0 = R0. This misleading result is caused by the assumption of small strain theory, which neglects the position changes during the deformation. Additionally, for the small strain rock mass, the stress state at the elastic-plastic boundary is independent of Young’s modulus and remains a constant value, i.e., 0.20 MPa and 0.34 MPa for the MC and HB rock masses. However, for the large strain behavior, the critical radial stress σrcr and tangential stress σθcr gradually increase with the increase in Young’s modulus and approach the limit values, which equal that of the small strain solution. For the soft rock mass with E = 0.1 MPa, σrcr = 5.37e3p0 andσθcr = 0.88p0 for the MC rock mass, while those of the HB rock mass are 0.12p0 and 0.92p0. The influence of Young’s modulus on Rp and Rs is shown in Fig. 12. For the small strain rock mass, with the decrease in Young’s modulus, Rp initially increases with a slow ratio. Once the residual region forms (logE < 7.7 for MC rock mass and logE < 8.0 for HB rock mass), both Rp and Rs significantly increase, and Rs tends to be equal to Rp. Then, both Rp and Rs maintain a relatively constant value with a very small increasing ratio with the decrease in Young’s modulus. For the employed examples, the relative stable softening radius is 4.44R0 and 12.85R0 for the MC and HB small strain rock masses. For the large strain solutions with a high Young’s modulus, Rp also gradually increases with the decrease in Young’s modulus. Once the residual region develops, Rp and Rs significantly increase, and Rp reaches its maximum value when E = 25.12 MPa and 39.81 MPa for the MC and HB rock masses, respectively. Then, Rp decreases with the further decrease in Young’s modulus and gradually becomes R0. The significant increase in plastic radius is induced by the development of residual region. However, when the Young’s modulus is notably small, the elastic deformation would be very large. A larger elastic deformation would leads to a small plastic region, since the position change is considered in the large strain
Fig. 12. Influence of Young's modulus on the radii of plastic and residual regions: (a) MC rock mass; (b) HB rock mass.
solution. A similar evolution can be found for Rs, and the maximum residual radius corresponds to E = 19.95 MPa for both the MC and HB rock mass. 6.3. Effects of Poisson’s ratio Poisson’s ratio also affects the displacement, as shown in Fig. 13. With the increase in Poisson’s ratio, u0, Rp and Rs gradually decrease. For the MC rock mass, when Poisson’s ratio increases from 0.01 to 0.499, u0 decreases from 0.37R0 to 0.21R0, Rp decreases from 2.28R0 to 2.05R0, and Rs decreases from 1.97R0 to 1.59R0, with a change of 43.24%, 10.09% and 19.29%, respectively. For the HB rock mass, u0, Rp and Rs decrease from 0.46R0, 3.13R0, and 2.16R0 to 0.28R0, 2.98R0 and 1.82R0, corresponding to decreases of 39.13%, 4.79% and 15.74%, respectively. 6.4. Effects of the critical plastic shear strain The elastic-perfectly plastic behavior and elastic-brittle-plastic behavior can be considered two extreme cases of strain-softening behavior. By assigning a critical plastic shear strain of γp* to very large and small values, the proposed strain-softening solutions can reduce to the above two behaviors. Fig. 14 shows the influence of γp* on u0, Rp and Rs. For the MC rock mass, u0 = 0.45R0 and Rp = 2.50R0 when γp* = 0, and those of 5.03 × 10−2R0 and 1.17R0 corresponds to the elastic-perfectly plastic behavior with γp* = 2. The displacement and softening regions decrease by 88.82% and 53.20%. When γp* > 0.24, no residual region occurs. Similarly, for the HB rock mass, u0 and Rp decrease from 0.72R0 and 3.66R0 to 0.07R0 and 1.57R0, respectively, corresponding to a decreases of 90.28% and 51.10%. Additionally, the residual region 9
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
0.36
2.3
0.32
2.1
0.28
1.9
(a)
0.5
0.24
u0/R0
1.7
Rs/R0 0.20 0.0
0.1
0.2
0.3
0.4
u0/R0 Rp/R0
0.4
u0/R0 Rp/R0
3.0
0.6
1.5 0.5
Rs/R0
1.5
0.2
1.0
0.1
0.5
0.0 0.0
0.4
0.8
(b)
3.6
0.45
u0/R0
2.8
Rp/R0 Rs/R0
0.35
u0/R0
u0/R0
R/R0
3.2
0.40
1.2
γ
0.50
2.4
1.6
0.25 0.0
2.0
0.1
0.2
0.3
0.4
0.9
4.0
0.8
u0/R0
3.5
0.7
Rp/R0 Rs/R0
3.0
0.6 0.5
2.5 2.0
0.4
1.5
0.2
1.0
0.1
0.5
0.0 0.0
1.6 0.5
0.0 2.0
p*
0.3
0.30
2.0
0.3
v
(b)
2.5
v
R/R0
2.5
R/R0
0.40
R/R0
u0/R0
(a)
0.4
0.8
1.2
γ
1.6
0.0 2.0
p*
Fig. 13. Influence of Poisson’s ratio on the GRCs and critical stress at the elastic-plastic boundary: (a) MC rock mass; (b) HB rock mass.
Fig. 14. Influence of the critical plastic shear strain on displacement, soften radius and residual radius: (a) MC rock mass, (b) HB rock mass.
disappear when γp* > 0.30. Unfortunately, neither of the two extreme conditions is reasonable because the elastic-perfectly plastic behavior overestimates the bearing capacity of the rock mass and the elasticbrittle-plastic behavior contrarily neglects the bearing capacity of the post-failure rock mass. Only a reasonable consideration of post-failure behavior can provide a theoretical foundation for the optimization of structure design and stability evaluation.
decrease, finally approaching the excavation radius. However, the softening and residual radii of the small strain solutions first significantly increase, and the residual radius tends to be equal to the softening radius; later, they maintain a relatively constant value. With the increase in Poisson’s ratio, the displacement, softening radius and residual radius gradually decrease. With the increase in the critical plastic shear strain, the displacement, softening radius and residual radius also gradually decrease, and the mechanical behavior gradually changes from an elastic-brittle-plastic model to an elastic-perfectly plastic model. Thus, the proposed large strain solutions are more generalized conditions for geotechnical engineering.
7. Conclusions Considering the real deformation behavior, a numerical large strain solution of a circular tunnel excavated in MC and HB strain-softening rock masses are derived with a nonassociated MC potential flow rule. An efficient calculation procedure with a five-order Runge-Kutta method is also presented. By simplifying the boundary conditions at the infinite radius with an assumed stress state, the stresses and displacements can be recursively determined. The proposed solutions are validated by the numerical simulation solutions. In addition, the effects on the elastic parameters and critical plastic shear strain are further investigated. The large strain displacement is smaller than that of the small strain solutions, and the displacement at the excavation radius is no larger than the excavation radius. For the cases without a residual region, the difference between the large strain and small strain solutions is very small with a relative error of displacement less than 10%. In addition, the softening radius of the large stain solution is approximately the same as that of the small strain solution. However, for the cases with a residual region, the displacement of the small strain solutions dramatically increases and may exceed the excavation radius when Young’s modulus is low. With the decrease in Young’s modulus, the softening and residual radii of the large strain solution first increase and then
Acknowledgments The authors are grateful to the financial support from the Fundamental Research Funds for the Central Universities (Grant No. 2018ZZCX04). References [1] Brown ET, Bray JW, Landayi B, Hoek E. Ground response curves for rock tunnels. ASCE J Geotech Eng Div 1983;109:15–39. [2] Zhang Q, Jiang BS, Wu XS, Zhang HQ, Han LJ. Elasto-plastic coupling analysis of circular openings in elasto-brittle-plastic rock mass. Theor Appl Fract Mec 2012;60:60–7. [3] Alejano LR, Alonso E. Considerations of the dilatancy angle in rocks and rock masses. Int J Rock Mech Min Sci 2005;42:481–507. [4] Lee YK, Pietruszczak S. A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass. Tunn Undergr Sp Tech 2008;23:588–99. [5] Park KH, Kim YJ. Analytical solution for a circular opening in an elastic- brittleplastic rock. Int J Rock Mech Min Sci 2006;43:616–22. [6] Park KH, Tontavanich B, Lee JG. A simple procedure for ground response curve of circular tunnel in elastic-strain softening rock masses. Tunn Undergr Sp Tech 2008;23:151–9.
10
Computers and Geotechnics 114 (2019) 103142
Q. Zhang, et al.
2014;38:1131–48. [18] Vrakas A. A finite strain solution for the elastoplastic ground response curve in tunnelling rocks with non-linear failure envelopes. Int J Numer Anal Met Geomech 2016;41:1077–90. [19] Vrakas A, Anagnostou G. A simple equation for obtaining finite strain solutions from small strain analyses of tunnels with very large convergences. Geotechnique 2015;65:936–44. [20] Park KH. Large strain similarity solution for a spherical or circular opening excavated in elastic-perfectly plastic media. Int J Numer Anal Met Geomech 2015;39:724–37. [21] Zhang Q, Li C, Jiang BS, Yu LY. Elastoplastic analysis of circular openings in elastobrittle-plastic rock mass based on logarithmic strain. Math Probl Eng 2017;7503912:1–9. [22] Guan K, Zhu WC, Wei J, Liu XG, Niu LL, Wang XR. A finite strain numerical procedure for a circular tunnel in strain-softening rock mass with large deformation. Int J Rock Mech Min Sci 2018;112:266–80. [23] Chadwick P. The quasi-static expansion of a spherical cavity in metals and ideal soils. Q J Mech Appl Math 1959;12:52–71. [24] Zhang Q, Jiang BS, Wang SL, Ge XR, Zhang HQ. Elasto-plastic analysis of a circular opening in strain-softening rock mass. Int J Rock Mech Min Sci 2012;50:38–46. [25] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [26] Ogawa T, Lo KY. Effects of dilatancy and yield criteria on displacements around tunnels. Can Geotech J 1987;24:100–13. [27] Varasa F, Alonso E, Alejano LR, Fdez.-Manína G. Study of bifurcation in the problem of unloading a circular excavation in a strain-softening material. Tunn Undergr Sp Tech 2005;20:311–22.
[7] Sharan SK. Elastic-brittle-plastic analysis of circular openings in Hoek-Brown media. Int J Rock Mech Min Sci 2003;40:817–24. [8] Wang SL, Yin XT, Tang H, Ge XR. A new approach for analyzing circular tunnel in strain-softening rock masses. Int J Rock Mech Min Sci 2010;47:170–8. [9] Zhang Q, Zhang CH, Jiang BS, Li N, Wang YC. Elastoplastic coupling solution of circular openings in strain-softening rock mass considering pressure-dependent effect. Int J Geomech, ASCE 2018;18:04017132. [10] Han JX, Li SC, Li SC, Yang WM. A procedure of strain-softening model for elastoplastic analysis of circular opening considering elasto-plastic coupling. Tunn Undergr Sp Tech 2013;37:128–34. [11] Brown ET, Bray JW, Santarelli FJ. Influence of stress-dependent elastic moduli on stresses and strains around axisymmetric boreholes. Rock Mech Rock Eng 1989;22:189–203. [12] Fahimifar A, Ghadami H, Ahmadvand M. An elasto-plastic model for underwater tunnels considering seepage body forces and strain-softening behavior. Eur J Environ Civ En 2015;19:129–51. [13] Zou JF, Qian ZH. Face-stability analysis of tunnels excavated below groundwater considering coupled flow deformation. Int J Geomech, ASCE 2018;18(8):04018089. [14] Lu AZ, Xu GS, Sun F, Sun WQ. Elasto-plastic analysis of a circular tunnel including the effects of the axial in situ stress. Int J Rock Mech Min Sci 2010;47:50–9. [15] Papanastasiou P, Durban D. Elastoplastic analysis of cylindrical cavity problems in geomaterials. Int J Numer Anal Met Geomech 1997;21:133–49. [16] Yu HS, Rowe RK. Plasticity solutions for soil behaviour around contracting cavities and tunnels. Int J Numer Anal Met Geomech 1999;23:1245–79. [17] Vrakas A, Anagnostou G. A finite strain closed-form solution for the elastoplastic ground response curve in tunnelling. Int J Numer Anal Met Geomech
11