- Email: [email protected]
Similarity solution for a spherical or circular opening in elastic-strain softening rock mass
International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
Contents lists available at ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Similarity solution for a spherical or circular opening in elastic-strain softening rock mass K.H. Park n College of Engineering & Technology, University of Dar es Salaam, P.O. Box 35131, Dar es Salaam, Tanzania
art ic l e i nf o
a b s t r a c t
Article history: Received 24 December 2013 Received in revised form 5 July 2014 Accepted 7 July 2014
This paper deals with the similarity solution for a spherical or circular opening excavated in elastic-strain softening rock mass compatible with a linear Mohr–Coulomb (M–C) or a nonlinear Hoek–Brown (H–B) yield criterion. A similarity solution for stresses and displacement is presented by replacing the partial differential equations from stress equilibrium, constitutive law, and consistency equations with firstorder ordinary differential equations. The Runge–Kutta (R–K) method is used to solve those first-order ordinary differential equations. Some measures are discussed to solve numerical instability problems in the use of R–K driver with adaptive steps. For comparison, the simple numerical stepwise procedure for a spherical opening is also presented by modifying the previous procedure for a circular opening. Three data sets are used to show the accuracy and practical application of the proposed methods. The results show the importance in choosing the solver for the system of ordinary differential equations and the initial values in the similarity solution. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Spherical opening Circular opening Similarity solution Elastic-strain softening rock Hoek–Brown yield criterion
1. Introduction Prediction of the stresses and displacements around a circular or spherical opening in the rock mass at great depth is an important problem in tunneling and underground space development such as the design of tunnels, boreholes and mine shafts. A number of closed-form analytical solutions and numerical procedures have been developed by considering the elastic– perfectly plastic, elastic–brittle plastic, and elastic-strain softening models of material behavior with the linear Mohr–Coulomb (M–C) and nonlinear Hoek–Brown (H–B) criteria [1–15]. For an elastic-strain softening model, Brown et al. [1] presented a numerical stepwise procedure for the stresses and displacements in the H–B media, while the procedure was further modified by Park et al. [10]. Lee and Pietruszczak [11] obtained the displacement and stress solutions of circular and spherical openings by gradually relieving the supporting pressure using a finite difference numerical procedure. Wang et al. [12,14] presented a numerical procedure by simplifying the strain-softening process as a series of stress drop and plastic flow, and then solving a series of elastic–brittle plastic solutions. On the other hand, the similarity solution method has been used for the analysis of a circular opening. In the method, partial
n
Tel.: þ 255 22 241 0128; fax: þ 255 22 241 0741. E-mail address: [email protected]
http://dx.doi.org/10.1016/j.ijrmms.2014.07.003 1365-1609/& 2014 Elsevier Ltd. All rights reserved.
differential equations are replaced by a small number of first-order ordinary differential equations, which can be solved by using relatively simple computational methods. Carranza-Torres and Fairhurst [4] presented the self-similarity solution for a circular tunnel excavated in elastic–perfectly plastic H–B rock mass. Carranza-Torres [7] extended the solution to the elastic–brittle plastic generalized H–B rock mass. The Runge–Kutta (R–K) method was used to integrate the second-order differential equation for radial displacement in the plastic region. For an elastic-strain softening model, Alonso et al. [6] proposed the self-similarity solution by solving the system of ordinary differential equations of equilibrium, persistence, radial displacement velocity and flow rule. The Runge–Kutta–Fehlberg (R–K–F) method was used to integrate the system of equations. However, the instable result in a certain case using the H–B yield criterion was reported [10]. The main objective of this paper is to develop similarity solution and numerical stepwise procedure for a spherical or circular opening excavated in elastic-strain softening rock mass compatible with M–C or H–B yield criterion. The paper consists of five parts. The first part deals with the definition of the problem and the necessary governing equations for spherical and circular openings. In the second part, the similarity solution is developed simultaneously for circular and spherical openings using H–B and M–C yield criteria. The system of first-order ordinary differential equations is obtained from the stress equilibrium, constitutive, and consistency equations. In the third part, the instability problems in numerical implementation of the similarity solution are discussed.
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
Strain Softening
3
Some measures are undertaken to overcome those instability problems. The fourth part describes the numerical stepwise procedure for a spherical opening, which is modified from Brown et al. [1] and Park et al. [10] methods for a circular opening. In the final part, the accuracy and practical application of the similarity solution and numerical procedure are illustrated by solving some examples. The results from similarity solution and numerical stepwise procedure are compared for three data sets using M–C and H–B yield criteria. The effect of η* (the value of the softening parameter controlling the transition between the softening and residual stages) in the similarity solution is investigated.
1
152
2. Definition of the problem Fig. 1 shows a circular or spherical opening being excavated in a continuous, homogeneous, isotropic, initially elastic rock mass subjected to a hydrostatic stress po. The opening surface is subjected to an internal pressure pi. As pi is gradually reduced, the radial displacement occurs and a plastic region develops around the opening when pi is less than the initial yield stress. The material behavior of elastic-strain softening model used in this study is shown in Fig. 2. There are three different zones around the opening: the elastic zone, the softening zone, and the residual zone. After initial yielding, the strength of rock drops gradually with increasing strain and follows the post-yield softening behavior. Two most commonly used yield criteria for the rock are considered in this study: the nonlinear H–B yield criterion, a σr f ðσ θ ; σ r ; ηÞ ¼ σ θ σ r σ ci mðηÞ þ sðηÞ ð1Þ
σ ci
and the linear M–C yield criterion, pffiffiffiffiffiffiffiffiffiffi f ðσ θ ; σ r ; ηÞ ¼ σ θ αðηÞσ r 2 cðηÞ αðηÞ
ð2Þ
with the plastic potential function such as, gðσ θ ; σ r ; ηÞ ¼ σ θ βðηÞσ r
ð3Þ
where σr is the radial stress, σθ is the circumferential stress, σci is the uniaxial compressive strength of the rock, m, s and a are the semi-empirical parameters that characterize the rock mass,
1
elastic
softening
residual
Fig. 2. Elastic-strain softening model [1,10].
α ¼ ð1 þ sin ϕÞ=ð1 sin ϕÞ, c is the cohesion of the rock, ϕ is the friction angle of the rock, β ¼ ð1 þ sin ψ Þ=ð1 sin ψ Þ, ψ is the dilation angle, and η is the softening parameter. For the deterioration process of material properties, it is assumed that the material properties decrease linearly with strain from their peak values at η ¼0 to the residual ones at η ¼ η*, such as [6,10] ( wp ðwp wr Þηηn for 0 o η o ηn w¼ ð4Þ wr for η Z ηn where w represents any one of the material parameters, such as c, ϕ, ψ , m and s, and η*¼the value of the softening parameter controlling the transition between the softening and residual stages. The subscripts p and r indicate the peak and residual values of the material parameters, respectively. The softening parameter η can be defined by two different ways, by a function of internal variables or in an incremental way. In this study the softening parameter is defined as the difference between the major and minor principal plastic strains, which reflects the plastic shear strain
η ¼ εpθ εpr ¼ εθ εr ðεeθ εer Þ
ð5Þ
in which superscripts e and p indicate the elastic and plastic parts, respectively. The solutions for circular and spherical openings can be developed simultaneously, with the value of the dimensional parameter k being 1 for a circular opening and 2 for a spherical opening. Three equations are needed to describe the contraction of the opening. 2.1. Stress equilibrium Assuming a state of plane strain or spherically symmetry around the opening, the equilibrium equation in polar coordinate system is given by ∂σ r k þ ðσ r σ θ Þ ¼ 0 ∂r r
ð6Þ
2.2. Constitutive law The standard linear relation for an isotropic elastic material can be written in the rate form Fig. 1. A circular or spherical opening in an infinite medium.
2Gε_ e ¼ Μσ_
ð7Þ
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
e _e _ _ e T is the vector of elastic strain rates, where ε ¼T εr ; εθ σ_ ¼ σ_ r ; σ_ θ is the vector of stress components, G is the shear modulus, superscript e indicates the elastic part, and M is a dimensionless compliance matrix such as, $ % M 11 M 12 1 ð2 kÞυ kυ 1 ð8Þ M¼ ¼ M 21 M 22 1 þðk 1Þν υ 1υ and υ is the Poisson ratio. The plastic strain rates are expressed in terms of the derivatives
of a plastic potential function g σ θ ; σ r ; η and given by ∂g ∂σ r
ε_ pr ¼ λ
kε_ pθ ¼ λ
!
p r pffiffiffiffiffiffiffiffiffiffi ðcp cr Þ cðηÞ 2 cos ϕðηÞ ðϕ ϕ Þ σ r þ pffiffiffiffiffiffiffiffiffiffi þ 2 αðηÞ
2 n η ηn αðηÞ 1 sin ϕðηÞ
∂f ¼ ∂η
ð23Þ d¼
∂g=∂σ θ 1 ¼ ∂g=∂σ r βðηÞ
ð24Þ
3. Similarity solution
∂g ∂σ θ
ð10Þ
kε_ p ∂g=∂σ r ¼ pθ ∂g=∂σ θ ε_ r
ð11Þ
where d is a simple function of dilation angle. Substituting Eq. (7) into Eq. (11) results in Ar σ_ r þ Aθ σ_ θ þ Au ε_ r ¼ Ao ε_ θ
ð12Þ
where Ar ¼
ð22Þ
ð9Þ
where λ is the plastic multiplier and superscript p indicates the plastic part. By eliminating the plastic multiplier in Eqs. (9) and (10), the non-associated flow rule can be expressed as d¼
∂f ¼1 ∂σ θ
153
1 k 1 k k M 11 þ M 21 ; Aθ ¼ M 12 þ M 22 ; Au ¼ 1; Ao ¼ 2G d 2G d d ð13Þ
2.3. Consistency condition A further rate equation, valid in the plastic region, is obtained by differentiation of the yield function associated with a given solid material with respect to time, such as ∂f ∂f ∂f σ_ r þ σ_ θ þ η_ ¼ 0 ð14Þ ∂σ r ∂σ θ ∂η By considering the rate of softening parameter as ∂η ∂η ∂η ∂η η_ ¼ σ_ r þ σ_ θ þ ε_ r þ ε_ ∂σ r ∂σ θ ∂ εr ∂ εθ θ
ð15Þ
the consistency equation can be written as Br σ_ r þ Bθ σ_ θ þ Bu ε_ r ¼ Bo ε_ θ
ð16Þ
where ∂f ∂f ∂η ∂f ∂f ∂η ∂f ∂η ∂f ∂η þ ; B ¼ þ ; Bu ¼ ; Bo ¼ Br ¼ ∂σ r ∂η ∂σ r θ ∂σ θ ∂ η ∂σ θ ∂η ∂εr ∂ η ∂ εθ ð17Þ The derivatives for the softening region are as follows. For the H–B yield criterion, a 1 ∂f σr ¼ 1 a mðηÞ mðηÞ þsðηÞ ð18Þ ∂σ r σ ci ∂f ¼1 ∂σ θ
ð19Þ
a 1 σr σ r ðmp mr Þ ðsp sr Þ mðηÞ þ sðηÞ þ n n
∂f ¼ a σ ci ∂η σ ci
σ ci
η
η
ð20Þ
and for the M–C yield criterion, ∂f ¼ αðηÞ ∂σ r
ð21Þ
In a similarity solution, the partial differential Eqs. (6), (12) and (16) that occur in the basic equilibrium, constitutive, and consistency equations are replaced by ordinary derivatives. Consider the single dimensionless radial coordinate ρ,
ρ¼
r r ¼ re τ
ð25Þ
and the non-dimensionalized forms of displacement, u, and stresses u~ ¼
2G u σr σ ; σ~ r ¼ ; σ~ θ ¼ θ p1 τ p1 p1
ð26Þ
where re is the radius at the elastic–plastic interface, τ is the fictitious time variable, and p1 is a suitably chosen reference, such as p1 ¼ po p1y
ð27Þ
and p1y is the initial yielding stress. From Eq. (25) and the chain rule for differentiation, it follows that the space and time derivatives occurring in the governing partial differential equations can be expressed in terms of the single ordinary derivative d/dρ, such as ∂ 1 d ρ d ¼ ; ðÞ ¼ ð28Þ ∂r r e dρ r e dρ Using Eq. (28), the stress equilibrium equation, Eq. (6) is replaced as dσ~ r k þ ðσ~ r σ~ θ Þ ¼ 0 dρ ρ Considering that p du~ u_ ¼ 1 u~ ρ dρ 2G
ε_ r ¼
2 ∂u_ p ρ d u~ ¼ 1 ∂r 2G r e dρ2
u_ r
ε_ θ ¼ ¼
p1 u~ du~ 2Gr e ρ dρ
ð29Þ
ð30Þ
ð31Þ
ð32Þ
the constitutive equation, Eq. (12), and the consistency equation, Eq. (16), are replaced as 2 dσ~ r ~ dσ~ θ A~ u d u~ A~ o 1 du~ u~ A~ r þ Aθ þ ¼ dρ dρ 2G dρ2 2G ρ dρ ρ2
ð33Þ
2 dσ~ r ~ dσ~ θ B~ u d u~ B~ o 1 du~ u~ þ Bθ þ ¼ 2 B~ r 2 dρ dρ 2G dρ 2G ρ dρ ρ
ð34Þ
As Eqs. (33) and (34) are the second-order equations, they can be expressed as a system of four first-order ordinary differential
154
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
where
equations as 2
1 6 A~ 6 r 6 6~ 4 Br 0
0 A~
0
B~ θ
A~ u 2G B~ u 2G
0
0
θ
9 8 8
9 3> dσ~ r > > k σ~ θ σ~ r > dρ > ρ > > 0 > > > > > > > >
> > > > > > > dσ~ θ > > 7 ~ A~ o 1 du~ u = = < < 0 7 dρ 2G ρ dρ ρ2 7 7> dυ~ > ¼ > B~ ~ ~ > 0 5> > > > dρ > > o 1 du u > > 2 > > > > > > > > > du~ > > 2G ρ dρ ρ > 1 : ; ; : dρ υ~
σo ¼
4. Numerical implementation of similarity solution
ð36aÞ
" # dσ~ θ 1 ðA~ o B~ u A~ u B~ o Þ υ~ u~ dσ~ r ~ ~ ~ ~ ¼ ðA r B u A u B r Þ 2G ρ ρ2 dρ dρ A~ θ B~ u A~ u B~ θ ð36bÞ " # dυ~ 2G ðA~ θ B~ o A~ o B~ θ Þ υ~ u~ dσ~ r ¼ 2 þ ðA~ r B~ θ A~ θ B~ r Þ 2G ρ ρ dρ A~ θ B~ u A~ u B~ θ dρ ð36cÞ du~ ¼ υ~ dρ
ð36dÞ
In the residual region, using ∂f ¼ 0-B~ u ¼ B~ o ¼ 0 ∂η
ð37Þ
Eqs. (36b) and (36c) are further simplified as dσ~ θ B~ r dσ~ r ¼ dρ B~ θ dρ " ! # dυ~ 2G A~ o υ~ u~ dσ~ r B~ r ¼ 2 A~ r A~ θ dρ A~ u 2G ρ ρ dρ B~ θ
ð38aÞ
the necessary initial values at
ð41Þ ð42Þ
ρ ¼1 can be obtained as
σ~ r ð1Þ ¼ σ o 1
ð43Þ
1 k
ð44Þ
σ~ θ ð1Þ ¼ σ o þ υ~ ð1Þ ¼ 1 ~ uð1Þ ¼
1 k
along with the fifth-order formula:
2825 18; 575 13; 525 277 1 k1 þ k3 þ k4 þ k5 þ k6 h 27; 648 48; 384 55; 296 14; 336 4
ð49Þ where
Using the solutions for stress and displacement in the outer elastic region given by
r 1 þ k σ r ¼ po þ ðp po Þ e ð40Þ r
ðp pÞ r e 1 þ k r u¼ o 2kG r
(1) In Eq. (18) of the H–B yield criterion, ∂f =∂σ r -1 as σr-0 in the case of sr ¼ 0. This problem happens in the residual region when σr r0, and can induce instable results as pi/po-0. This case was reported in the previous study for a circular opening when applying the data Set 1 [10]. (2) The numerical instability problem due to small step size and stiffness can arise at the elastic–softening and softening– residual interfaces. The R–K–F method uses the following fourth-order estimate [16]: 37 250 125 512 yi þ 1 ¼ yi þ k1 þ k3 þ k4 þ k6 h ð48Þ 378 621 594 1771
ð38bÞ
3.1. Initial conditions at the elastic–plastic boundary
ðp po Þ r e 1 þ k k r
In order to obtain the similarity solution from Eqs. (36a)–(36d), (38a) and (38b), many standard ordinary differential equation system solver library routines are available. Alonso et al. [6] used the R–K–F method and the program in the MATLAB environment under the function name ‘ode45’. Although Alonso et al. [6] showed good results for the M–C yield criterion, the following instability problems should be addressed for the H–B yield criterion in numerical implementation:
yi þ 1 ¼ yi þ
Eqs. (36a)–(36d), (38a) and (38b) can be solved as an initial value problem starting at the elastic–plastic boundary (ρ ¼ 1). Many standard ordinary differential equation system solver library routines are available for this purpose. The softening parameter η is also replaced as p u~ ð39Þ η~ ¼ 1 υ~ þ ðM11 M21 Þσ~ r þ ðM12 M22 Þσ~ θ 2G ρ
σ θ ¼ po
ð47Þ
ð35Þ
or simply dσ~ r k ¼ ðσ~ σ~ r Þ dρ ρ θ
po po p1y
ð45Þ ð46Þ
k1 ¼ f ðxi ; yi Þ ¼
dy ðx ; y Þ dx i i
ð50aÞ
1 1 k2 ¼ f xi þ h; yi þ k1 h 5 5
ð50bÞ
3 3 9 k3 ¼ f xi þ h; yi þ k1 h þ k2 h 10 40 40
ð50cÞ
3 3 9 6 k4 ¼ f xi þ h; yi þ k1 h k2 h þ k3 h 5 10 10 5
ð50dÞ
11 5 70 35 k5 ¼ f xi þ h; yi k1 h þ k2 h k3 h þ k4 h 54 2 27 27
ð50eÞ
7 1631 175 575 k1 h þ k2 h þ k3 h k6 ¼ f xi þ h; yi þ 8 55; 296 512 13; 824 44; 275 253 k4 hþ k5 h þ 110; 592 4096
ð50fÞ
and h is the step size. First, the slopes for all variables at the initial value (xi) are obtained. These slopes (a set of k1's) are then used to make predictions of dependent variable at the point of xi þh/5. These point values are in turn used to compute a set of slopes at that point (the k2's). A set of k3–k6 is continuously obtained at different points between xi and xi þ 1. Finally, the k's are combined to make the final prediction in Eqs. (48) and (49). If the elastic–softening or softening–residual interface is located between xi and xi þ 1, the large local truncation error can be estimated. In order to adjust the step size, a certain criterion can be used, for example the R–K driver with adaptive steps (subroutine odeint, available in [17]). However, if the system of
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
equation involves rapidly changing components together with slowly changing ones, the stiffness problem can arise in the solution. It might suppose that the adaptive step-size routine might offer a solution for this problem. However, this is not the case, because the stability requirement will still necessitate using very small steps throughout the entire solution [16]. The error in the softening region can be accumulated in the calculation of the residual region. (3) In the case of very small value of η* (that is, the elastic–brittle plastic behavior), the numerical instability can arise because initial values of σ~ θ ð1Þ and υ~ ð1Þ, such as Eqs. (44) and (45), are not valid in the elastic–brittle plastic model. In order to overcome those problems, in this study the classical fourth-order R–K method is used without any adaptive step size yi þ 1 ¼ yi þ 16 ðk1 þ 2k2 þ2k3 þ k4 Þh k1 ¼ f ðxi ; yi Þ
ð52aÞ
1 1 k2 ¼ f xi þ h; yi þ k1 h 2 2
ð52bÞ
1 1 k3 ¼ f xi þ h; yi þ k2 h 2 2
ð52cÞ
k4 ¼ f ðxi þ h; yi þ k3 hÞ
ð52dÞ
The subroutine for the classical fourth-order R–K method is further modified to consider the case of σr r 0 in the residual region. First the prediction of yi þ 1 is obtained at xi using Euler's method, ð53Þ
If σr(i þ 1) r0, then the values of ρi þ 1 and ui þ 1 are obtained by the interpolation. In order to consider the initial values for very small value of η*, the first prediction of η(2) is obtained from the prediction of y(2) using Euler's method. If η(2) Z η*, then the following initial values for the elastic–brittle plastic model, instead of Eqs. (44) and (45), are used: a σ m p σ~ ð1Þ σ~ θ ð1Þ ¼ σ~ r ð1Þ þ ci r 1 r þ sr ð54Þ p1 σ ci ~ þ 2GA~ r ðσ~ r ð1Þ p~ o Þ þ 2GA~ θ ðσ~ θ ð1Þ p~ o Þ υ~ ð1Þ ¼ Ao uð1Þ
ð55Þ
The use of Eq. (55) was verified in [18].
5. Numerical stepwise procedure The simple numerical stepwise procedure, proposed in the previous studies [1,10], is modified to obtain the solution for a spherical opening. In the same way, the dimensional parameter k is used as 1 for a circular opening and 2 for a spherical opening. By dividing the plastic region with a number of thin annular rings (Fig. 3), one can obtain the radial stress at rj for the H–B yield criterion [1,10] qffiffiffiffiffiffiffiffiffiffiffiffiffi σ r ðjÞ ¼ b b2 a ð56Þ where a ¼ σ 2rðj 1Þ 4ℓk
2 1 2m c
b ¼ σ r ðj 1Þ þ ℓk mσ c 2
ℓ¼
r j 1 r j r j 1 þr j
σ σ rðj 1Þ þ sσ 2c
Fig. 3. Typical annulus in plastic region.
ð51Þ
yi þ 1 ¼ yi þ f ðxi ; yi Þh
155
ð57Þ ð58Þ
2 ð59Þ
Table 1 Data set (H–B criterion).
Young's modulus, E (MPa) Poisson's ratio, ν Initial stress, po (MPa) Radius of tunnel, ri (m) σci (MPa) mp sp mr sr ψ peak (deg) ψ r (deg) η
Set 1 [1]
Set 2 [6]
1380 0.25 3.31 5.35 27.6 0.5 0.001 0.1 0 19.47 5.22 0.004742
1380 0.25 3.31 3.0 3.31 0.5 0.001 0.1 0.00001 15.0 5.0 0.0125
using the radii given by 2εθðj 1Þ εrðj 1Þ εrðjÞ rj ¼ rj 1 2εθðjÞ εrðj 1Þ εrðjÞ
ð60Þ
Let the circumferential strain increment be dεθ ¼ εθðjÞ εθðj 1Þ
ð61Þ
Then, the radial strain increment at rj can be obtained as dεrðjÞ ¼ dεerðj 1Þ βðdεθðjÞ dεeθðj 1Þ Þ
ð62Þ
and dεθðj 1Þ are the radial and the circumferential where dε elastic strain increments at rj 1, respectively, and β ¼ ð1 þ
sin ψ Þ= 1 sin ψ . By selecting an arbitrary small value for dεθ, the values of εθ(j), εr(j), rj and uj can be calculated using Eqs. (61), (62) and (60) in that order. The process is then repeated several times to determine completely the strains and displacements in the plastic region.The radial stress at the elastic–plastic interface, in which r¼r1 ¼ re, is given by e rðj 1Þ
e
σ re ¼ po Mσ c
ð63Þ
in which
8sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9 2 < 2 = m p 1 k 1 þ k p o þ sp mp M¼ m2p þ 4 ; 2 1þk : k σc
ð64Þ
Using Eq. (63) as the starting point, successive values of σr(j) can be calculated from Eq. (56) for the radii determined from Eq. (60). The softening parameter at rj can be approximately estimated as
η ¼ εθðj 1Þ εrðj 1Þ þ If η exceeds condition.
1 ðM 11 M 21 Þσ rðj 1Þ ðM 12 M 22 Þσ θðj 1Þ 2G ð65Þ
η*, the strength of the rock mass reaches residual
156
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
The proposed procedure can be easily implemented into the computer program to obtain the ground response curve, by simply modifying the program in reference [10].
Table 2 Data Set 3 (M–C criterion). Radius of tunnel, ri Initial stress, po Young's modulus, E Poisson's ratio, ν η*
3.0 m 20 MPa 10 GPa 0.25 0.008
cp ϕp cr ϕr ψp ¼ ψr
1.0 MPa 301 0.7 MPa 221 3.751
6. Application In order to show the applicability and accuracy of the proposed solutions, the results of the radial displacements, the variations of re/ri and rs/ri, (which can show the spreading of the plastic radius with decreasing pressures) and the radial and circumferential stresses in the plastic region are compared using three data sets. Tables 1 and 2 show the data set for H–B and M–C criteria, respectively. The step size, h¼0.01, is used for similarity solution. Figs. 4–6 show the results of the radial displacements at the opening surface and the variations of re/ri and rs/ri, and the radial and circumferential stresses within the plastic region for each data set. In the figures, SS indicates the results from similarity solution,
Fig. 4. Results for data Set 1. (a) Ground response curve. (b) Variation of re/ri. (c) Variation of rs/ri. (d) Variation of radial stress for pi ¼ 0. (e) Variation of circumferential stress for pi ¼0.
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
while FD indicates those from numerical stepwise procedure. From these figures, the following observations can be made: (1) Good agreement between results from similarity solution and numerical stepwise procedure can be seen for data Set 2 and Set 3 in both circular and spherical openings. (2) For data Set 1, good agreement of the ground response curve between two solutions is obtained. It is noticed that the result of ground response curve from similarity solution in this study is better than that of the previous study (Fig. 4 in [10]). However, the differences in the variations of rs/ri and circumferential stress for pi ¼ 0 from two solutions are shown in Fig. 4 (c) and (e) for both circular and spherical openings. As can be seen in Fig. 4(d) and (e), the variation of radial stresses in the softening region is soft, while the circumferential stresses in
157
the softening region is suddenly changed. The results of similarity solution are obtained using the initial values of Eqs. (54) and (55), while the numerical stepwise procedure shows the elastic-strain softening behavior. (3) Considering the data sets used in this study, the displacements induced by opening from the circular solution are about 4.3– 4.7 times larger than those from the spherical solution.
6.1. Effect of η* As mentioned in Section 4 and shown in Fig. 4(e), the numerical instability can arise in the similarity solution because of very small value of η* or sharp change of circumferential stress in the softening region. The effect of η* is investigated using the
Fig. 5. Results for data Set 2. (a) Ground response curve. (b) Variation of re/ri. (c) Variation of rs/ri. (d) Variation of radial stress for pi ¼ 0. (e) Variation of circumferential stress for pi ¼ 0.
158
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
Fig. 6. Results for data Set 3. (a) Ground response curve. (b) Variation of re/ri. (c) Variation of rs/ri. (d) Variation of radial stress for pi ¼ 0. (e) Variation of circumferential stress for pi ¼0.
data Set 1. Fig. 7 shows the results of ground response curve and circumferential stress for pi ¼0 with different values of η*. As η* decreases, the displacement at the opening increases and the slope of circumferential stress in the softening region is sharply changed. When η* ¼0.005, the results are obtained using the initial values of Eqs. (54) and (55).
7. Conclusions The similarity solution and numerical stepwise procedure for stresses and displacement have been presented for circular and spherical openings excavated in elastic-strain softening rock
masses compatible with M–C and H–B yield criteria. The accuracy and practical application of the proposed methods are illustrated by solving some examples. Although good agreement from similarity solution and numerical stepwise procedure can be obtained for the M–C yield criterion, the numerical instability can arise for the H–B yield criterion in the similarity solution. In order to overcome those instability problems, some measures are undertaken as follows: (1) The classical fourthorder R–K method without any adaptive step size is used to solve the system of ordinary differential equations. (2) Euler's method is used to consider the case of σr r0 in the residual region. (3) The initial values of σ~ θ ð1Þ and υ~ ð1Þ, such as Eqs. (54) and (55), are used to consider the case of very small value of η*.
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
159
Fig. 7. Effect of η* (data Set 1). (a) Ground response curve. (b) Variation of circumferential stress for pi ¼ 0.
Although the methods proposed in this paper are limited in scope, it appears to be useful for the preliminary design of a spherical or circular rock opening to obtain the ground response curve. The methods can also be used for the validation of numerical methods, such as the finite element method. References [1] Brown ET, Bray JW, Ladanyi B, Hoek E. Ground response curves for rock tunnels. J Geotech Eng ASCE 1983;109:15–39. [2] Detournay E. Elastoplastic model of a deep tunnel for a rock with variable dilatancy. Rock Mech Rock Eng 1986;19:99–108. [3] Wang Y. Ground response of circular tunnel in poorly consolidated rock. J Geotech Eng ASCE 1996;122:703–1008. [4] Carranza-Torres C, Fairhurst C. The elasto-plastic response of underground excavations in rock masses that satisfy the Hoek–Brown failure criterion. Int J Rock Mech Min Sci 1999;36:777–809. [5] Sharan SK. Elastic–brittle–plastic analysis of circular openings in Hoek–Brown media. Int J Rock Mech Min Sci 2003;40:817–24. [6] Alonso E, Alejano LR, Varas F, Fdez-Manin G, Carranza-Torres C. Ground response curves for rock masses exhibiting strain-softening behavior. Int J Numer Anal Methods Geomech 2003;27:1153–85. [7] Carranza-Torres C. Elasto-plastic solution of tunnel problems using the generalized form of the Hoek–Brown failure criterion. Int J Rock Mech Min Sci 2004;41:480–1.
[8] Sharan SK. Exact and approximate solutions for displacements around circular openings in elastic–brittle–plastic Hoek–Brown rock. Int J Rock Mech Min Sci 2005;42:542–9. [9] Park KH, Kim YJ. Analytical solution for a circular opening in an elastic–brittle– plastic rock. Int J Rock Mech Min Sci 2006;43:616–22. [10] Park KH, Tontavanich B, Lee JG. A simple procedure for ground response curve of circular tunnel in elastic-strain softening rock masses. Tunn Underg Space Technol 2008;23:151–9. [11] Lee YK, Pietruszczak S. A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass. Tunn Underg Space Technol 2008;23:588–99. [12] 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. [13] Wang SL, Yin SD. A closed-form solution for a spherical cavity in the elastic– brittle–plastic medium. Tunn Underg Space Technol 2011;26:236–41. [14] Wang SL, Zheng H, Li CG, Ge XR. A finite element implementation of strainsoftening rock mass. Int J Rock Mech Min Sci 2011;48:67–76. [15] Wang SL, Yin SD, Wu ZJ. Strain-softening analysis of a spherical cavity. Int J Numer Anal Methods Geomech 2012;36:182–202. [16] Chapra SC, Canale RP. Numerical methods for engineers. 4th ed.. New York: McGraw-Hill; 2002. [17] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran 77. Cambridge: Cambridge University Press; 1997. [18] Park KH. Similarity solution of a spherical or circular opening in elastic–brittle plastic rock mass. 2014. [in preparation].
Contents lists available at ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Similarity solution for a spherical or circular opening in elastic-strain softening rock mass K.H. Park n College of Engineering & Technology, University of Dar es Salaam, P.O. Box 35131, Dar es Salaam, Tanzania
art ic l e i nf o
a b s t r a c t
Article history: Received 24 December 2013 Received in revised form 5 July 2014 Accepted 7 July 2014
This paper deals with the similarity solution for a spherical or circular opening excavated in elastic-strain softening rock mass compatible with a linear Mohr–Coulomb (M–C) or a nonlinear Hoek–Brown (H–B) yield criterion. A similarity solution for stresses and displacement is presented by replacing the partial differential equations from stress equilibrium, constitutive law, and consistency equations with firstorder ordinary differential equations. The Runge–Kutta (R–K) method is used to solve those first-order ordinary differential equations. Some measures are discussed to solve numerical instability problems in the use of R–K driver with adaptive steps. For comparison, the simple numerical stepwise procedure for a spherical opening is also presented by modifying the previous procedure for a circular opening. Three data sets are used to show the accuracy and practical application of the proposed methods. The results show the importance in choosing the solver for the system of ordinary differential equations and the initial values in the similarity solution. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Spherical opening Circular opening Similarity solution Elastic-strain softening rock Hoek–Brown yield criterion
1. Introduction Prediction of the stresses and displacements around a circular or spherical opening in the rock mass at great depth is an important problem in tunneling and underground space development such as the design of tunnels, boreholes and mine shafts. A number of closed-form analytical solutions and numerical procedures have been developed by considering the elastic– perfectly plastic, elastic–brittle plastic, and elastic-strain softening models of material behavior with the linear Mohr–Coulomb (M–C) and nonlinear Hoek–Brown (H–B) criteria [1–15]. For an elastic-strain softening model, Brown et al. [1] presented a numerical stepwise procedure for the stresses and displacements in the H–B media, while the procedure was further modified by Park et al. [10]. Lee and Pietruszczak [11] obtained the displacement and stress solutions of circular and spherical openings by gradually relieving the supporting pressure using a finite difference numerical procedure. Wang et al. [12,14] presented a numerical procedure by simplifying the strain-softening process as a series of stress drop and plastic flow, and then solving a series of elastic–brittle plastic solutions. On the other hand, the similarity solution method has been used for the analysis of a circular opening. In the method, partial
n
Tel.: þ 255 22 241 0128; fax: þ 255 22 241 0741. E-mail address: [email protected]
http://dx.doi.org/10.1016/j.ijrmms.2014.07.003 1365-1609/& 2014 Elsevier Ltd. All rights reserved.
differential equations are replaced by a small number of first-order ordinary differential equations, which can be solved by using relatively simple computational methods. Carranza-Torres and Fairhurst [4] presented the self-similarity solution for a circular tunnel excavated in elastic–perfectly plastic H–B rock mass. Carranza-Torres [7] extended the solution to the elastic–brittle plastic generalized H–B rock mass. The Runge–Kutta (R–K) method was used to integrate the second-order differential equation for radial displacement in the plastic region. For an elastic-strain softening model, Alonso et al. [6] proposed the self-similarity solution by solving the system of ordinary differential equations of equilibrium, persistence, radial displacement velocity and flow rule. The Runge–Kutta–Fehlberg (R–K–F) method was used to integrate the system of equations. However, the instable result in a certain case using the H–B yield criterion was reported [10]. The main objective of this paper is to develop similarity solution and numerical stepwise procedure for a spherical or circular opening excavated in elastic-strain softening rock mass compatible with M–C or H–B yield criterion. The paper consists of five parts. The first part deals with the definition of the problem and the necessary governing equations for spherical and circular openings. In the second part, the similarity solution is developed simultaneously for circular and spherical openings using H–B and M–C yield criteria. The system of first-order ordinary differential equations is obtained from the stress equilibrium, constitutive, and consistency equations. In the third part, the instability problems in numerical implementation of the similarity solution are discussed.
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
Strain Softening
3
Some measures are undertaken to overcome those instability problems. The fourth part describes the numerical stepwise procedure for a spherical opening, which is modified from Brown et al. [1] and Park et al. [10] methods for a circular opening. In the final part, the accuracy and practical application of the similarity solution and numerical procedure are illustrated by solving some examples. The results from similarity solution and numerical stepwise procedure are compared for three data sets using M–C and H–B yield criteria. The effect of η* (the value of the softening parameter controlling the transition between the softening and residual stages) in the similarity solution is investigated.
1
152
2. Definition of the problem Fig. 1 shows a circular or spherical opening being excavated in a continuous, homogeneous, isotropic, initially elastic rock mass subjected to a hydrostatic stress po. The opening surface is subjected to an internal pressure pi. As pi is gradually reduced, the radial displacement occurs and a plastic region develops around the opening when pi is less than the initial yield stress. The material behavior of elastic-strain softening model used in this study is shown in Fig. 2. There are three different zones around the opening: the elastic zone, the softening zone, and the residual zone. After initial yielding, the strength of rock drops gradually with increasing strain and follows the post-yield softening behavior. Two most commonly used yield criteria for the rock are considered in this study: the nonlinear H–B yield criterion, a σr f ðσ θ ; σ r ; ηÞ ¼ σ θ σ r σ ci mðηÞ þ sðηÞ ð1Þ
σ ci
and the linear M–C yield criterion, pffiffiffiffiffiffiffiffiffiffi f ðσ θ ; σ r ; ηÞ ¼ σ θ αðηÞσ r 2 cðηÞ αðηÞ
ð2Þ
with the plastic potential function such as, gðσ θ ; σ r ; ηÞ ¼ σ θ βðηÞσ r
ð3Þ
where σr is the radial stress, σθ is the circumferential stress, σci is the uniaxial compressive strength of the rock, m, s and a are the semi-empirical parameters that characterize the rock mass,
1
elastic
softening
residual
Fig. 2. Elastic-strain softening model [1,10].
α ¼ ð1 þ sin ϕÞ=ð1 sin ϕÞ, c is the cohesion of the rock, ϕ is the friction angle of the rock, β ¼ ð1 þ sin ψ Þ=ð1 sin ψ Þ, ψ is the dilation angle, and η is the softening parameter. For the deterioration process of material properties, it is assumed that the material properties decrease linearly with strain from their peak values at η ¼0 to the residual ones at η ¼ η*, such as [6,10] ( wp ðwp wr Þηηn for 0 o η o ηn w¼ ð4Þ wr for η Z ηn where w represents any one of the material parameters, such as c, ϕ, ψ , m and s, and η*¼the value of the softening parameter controlling the transition between the softening and residual stages. The subscripts p and r indicate the peak and residual values of the material parameters, respectively. The softening parameter η can be defined by two different ways, by a function of internal variables or in an incremental way. In this study the softening parameter is defined as the difference between the major and minor principal plastic strains, which reflects the plastic shear strain
η ¼ εpθ εpr ¼ εθ εr ðεeθ εer Þ
ð5Þ
in which superscripts e and p indicate the elastic and plastic parts, respectively. The solutions for circular and spherical openings can be developed simultaneously, with the value of the dimensional parameter k being 1 for a circular opening and 2 for a spherical opening. Three equations are needed to describe the contraction of the opening. 2.1. Stress equilibrium Assuming a state of plane strain or spherically symmetry around the opening, the equilibrium equation in polar coordinate system is given by ∂σ r k þ ðσ r σ θ Þ ¼ 0 ∂r r
ð6Þ
2.2. Constitutive law The standard linear relation for an isotropic elastic material can be written in the rate form Fig. 1. A circular or spherical opening in an infinite medium.
2Gε_ e ¼ Μσ_
ð7Þ
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
e _e _ _ e T is the vector of elastic strain rates, where ε ¼T εr ; εθ σ_ ¼ σ_ r ; σ_ θ is the vector of stress components, G is the shear modulus, superscript e indicates the elastic part, and M is a dimensionless compliance matrix such as, $ % M 11 M 12 1 ð2 kÞυ kυ 1 ð8Þ M¼ ¼ M 21 M 22 1 þðk 1Þν υ 1υ and υ is the Poisson ratio. The plastic strain rates are expressed in terms of the derivatives
of a plastic potential function g σ θ ; σ r ; η and given by ∂g ∂σ r
ε_ pr ¼ λ
kε_ pθ ¼ λ
!
p r pffiffiffiffiffiffiffiffiffiffi ðcp cr Þ cðηÞ 2 cos ϕðηÞ ðϕ ϕ Þ σ r þ pffiffiffiffiffiffiffiffiffiffi þ 2 αðηÞ
2 n η ηn αðηÞ 1 sin ϕðηÞ
∂f ¼ ∂η
ð23Þ d¼
∂g=∂σ θ 1 ¼ ∂g=∂σ r βðηÞ
ð24Þ
3. Similarity solution
∂g ∂σ θ
ð10Þ
kε_ p ∂g=∂σ r ¼ pθ ∂g=∂σ θ ε_ r
ð11Þ
where d is a simple function of dilation angle. Substituting Eq. (7) into Eq. (11) results in Ar σ_ r þ Aθ σ_ θ þ Au ε_ r ¼ Ao ε_ θ
ð12Þ
where Ar ¼
ð22Þ
ð9Þ
where λ is the plastic multiplier and superscript p indicates the plastic part. By eliminating the plastic multiplier in Eqs. (9) and (10), the non-associated flow rule can be expressed as d¼
∂f ¼1 ∂σ θ
153
1 k 1 k k M 11 þ M 21 ; Aθ ¼ M 12 þ M 22 ; Au ¼ 1; Ao ¼ 2G d 2G d d ð13Þ
2.3. Consistency condition A further rate equation, valid in the plastic region, is obtained by differentiation of the yield function associated with a given solid material with respect to time, such as ∂f ∂f ∂f σ_ r þ σ_ θ þ η_ ¼ 0 ð14Þ ∂σ r ∂σ θ ∂η By considering the rate of softening parameter as ∂η ∂η ∂η ∂η η_ ¼ σ_ r þ σ_ θ þ ε_ r þ ε_ ∂σ r ∂σ θ ∂ εr ∂ εθ θ
ð15Þ
the consistency equation can be written as Br σ_ r þ Bθ σ_ θ þ Bu ε_ r ¼ Bo ε_ θ
ð16Þ
where ∂f ∂f ∂η ∂f ∂f ∂η ∂f ∂η ∂f ∂η þ ; B ¼ þ ; Bu ¼ ; Bo ¼ Br ¼ ∂σ r ∂η ∂σ r θ ∂σ θ ∂ η ∂σ θ ∂η ∂εr ∂ η ∂ εθ ð17Þ The derivatives for the softening region are as follows. For the H–B yield criterion, a 1 ∂f σr ¼ 1 a mðηÞ mðηÞ þsðηÞ ð18Þ ∂σ r σ ci ∂f ¼1 ∂σ θ
ð19Þ
a 1 σr σ r ðmp mr Þ ðsp sr Þ mðηÞ þ sðηÞ þ n n
∂f ¼ a σ ci ∂η σ ci
σ ci
η
η
ð20Þ
and for the M–C yield criterion, ∂f ¼ αðηÞ ∂σ r
ð21Þ
In a similarity solution, the partial differential Eqs. (6), (12) and (16) that occur in the basic equilibrium, constitutive, and consistency equations are replaced by ordinary derivatives. Consider the single dimensionless radial coordinate ρ,
ρ¼
r r ¼ re τ
ð25Þ
and the non-dimensionalized forms of displacement, u, and stresses u~ ¼
2G u σr σ ; σ~ r ¼ ; σ~ θ ¼ θ p1 τ p1 p1
ð26Þ
where re is the radius at the elastic–plastic interface, τ is the fictitious time variable, and p1 is a suitably chosen reference, such as p1 ¼ po p1y
ð27Þ
and p1y is the initial yielding stress. From Eq. (25) and the chain rule for differentiation, it follows that the space and time derivatives occurring in the governing partial differential equations can be expressed in terms of the single ordinary derivative d/dρ, such as ∂ 1 d ρ d ¼ ; ðÞ ¼ ð28Þ ∂r r e dρ r e dρ Using Eq. (28), the stress equilibrium equation, Eq. (6) is replaced as dσ~ r k þ ðσ~ r σ~ θ Þ ¼ 0 dρ ρ Considering that p du~ u_ ¼ 1 u~ ρ dρ 2G
ε_ r ¼
2 ∂u_ p ρ d u~ ¼ 1 ∂r 2G r e dρ2
u_ r
ε_ θ ¼ ¼
p1 u~ du~ 2Gr e ρ dρ
ð29Þ
ð30Þ
ð31Þ
ð32Þ
the constitutive equation, Eq. (12), and the consistency equation, Eq. (16), are replaced as 2 dσ~ r ~ dσ~ θ A~ u d u~ A~ o 1 du~ u~ A~ r þ Aθ þ ¼ dρ dρ 2G dρ2 2G ρ dρ ρ2
ð33Þ
2 dσ~ r ~ dσ~ θ B~ u d u~ B~ o 1 du~ u~ þ Bθ þ ¼ 2 B~ r 2 dρ dρ 2G dρ 2G ρ dρ ρ
ð34Þ
As Eqs. (33) and (34) are the second-order equations, they can be expressed as a system of four first-order ordinary differential
154
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
where
equations as 2
1 6 A~ 6 r 6 6~ 4 Br 0
0 A~
0
B~ θ
A~ u 2G B~ u 2G
0
0
θ
9 8 8
9 3> dσ~ r > > k σ~ θ σ~ r > dρ > ρ > > 0 > > > > > > > >
> > > > > > > dσ~ θ > > 7 ~ A~ o 1 du~ u = = < < 0 7 dρ 2G ρ dρ ρ2 7 7> dυ~ > ¼ > B~ ~ ~ > 0 5> > > > dρ > > o 1 du u > > 2 > > > > > > > > > du~ > > 2G ρ dρ ρ > 1 : ; ; : dρ υ~
σo ¼
4. Numerical implementation of similarity solution
ð36aÞ
" # dσ~ θ 1 ðA~ o B~ u A~ u B~ o Þ υ~ u~ dσ~ r ~ ~ ~ ~ ¼ ðA r B u A u B r Þ 2G ρ ρ2 dρ dρ A~ θ B~ u A~ u B~ θ ð36bÞ " # dυ~ 2G ðA~ θ B~ o A~ o B~ θ Þ υ~ u~ dσ~ r ¼ 2 þ ðA~ r B~ θ A~ θ B~ r Þ 2G ρ ρ dρ A~ θ B~ u A~ u B~ θ dρ ð36cÞ du~ ¼ υ~ dρ
ð36dÞ
In the residual region, using ∂f ¼ 0-B~ u ¼ B~ o ¼ 0 ∂η
ð37Þ
Eqs. (36b) and (36c) are further simplified as dσ~ θ B~ r dσ~ r ¼ dρ B~ θ dρ " ! # dυ~ 2G A~ o υ~ u~ dσ~ r B~ r ¼ 2 A~ r A~ θ dρ A~ u 2G ρ ρ dρ B~ θ
ð38aÞ
the necessary initial values at
ð41Þ ð42Þ
ρ ¼1 can be obtained as
σ~ r ð1Þ ¼ σ o 1
ð43Þ
1 k
ð44Þ
σ~ θ ð1Þ ¼ σ o þ υ~ ð1Þ ¼ 1 ~ uð1Þ ¼
1 k
along with the fifth-order formula:
2825 18; 575 13; 525 277 1 k1 þ k3 þ k4 þ k5 þ k6 h 27; 648 48; 384 55; 296 14; 336 4
ð49Þ where
Using the solutions for stress and displacement in the outer elastic region given by
r 1 þ k σ r ¼ po þ ðp po Þ e ð40Þ r
ðp pÞ r e 1 þ k r u¼ o 2kG r
(1) In Eq. (18) of the H–B yield criterion, ∂f =∂σ r -1 as σr-0 in the case of sr ¼ 0. This problem happens in the residual region when σr r0, and can induce instable results as pi/po-0. This case was reported in the previous study for a circular opening when applying the data Set 1 [10]. (2) The numerical instability problem due to small step size and stiffness can arise at the elastic–softening and softening– residual interfaces. The R–K–F method uses the following fourth-order estimate [16]: 37 250 125 512 yi þ 1 ¼ yi þ k1 þ k3 þ k4 þ k6 h ð48Þ 378 621 594 1771
ð38bÞ
3.1. Initial conditions at the elastic–plastic boundary
ðp po Þ r e 1 þ k k r
In order to obtain the similarity solution from Eqs. (36a)–(36d), (38a) and (38b), many standard ordinary differential equation system solver library routines are available. Alonso et al. [6] used the R–K–F method and the program in the MATLAB environment under the function name ‘ode45’. Although Alonso et al. [6] showed good results for the M–C yield criterion, the following instability problems should be addressed for the H–B yield criterion in numerical implementation:
yi þ 1 ¼ yi þ
Eqs. (36a)–(36d), (38a) and (38b) can be solved as an initial value problem starting at the elastic–plastic boundary (ρ ¼ 1). Many standard ordinary differential equation system solver library routines are available for this purpose. The softening parameter η is also replaced as p u~ ð39Þ η~ ¼ 1 υ~ þ ðM11 M21 Þσ~ r þ ðM12 M22 Þσ~ θ 2G ρ
σ θ ¼ po
ð47Þ
ð35Þ
or simply dσ~ r k ¼ ðσ~ σ~ r Þ dρ ρ θ
po po p1y
ð45Þ ð46Þ
k1 ¼ f ðxi ; yi Þ ¼
dy ðx ; y Þ dx i i
ð50aÞ
1 1 k2 ¼ f xi þ h; yi þ k1 h 5 5
ð50bÞ
3 3 9 k3 ¼ f xi þ h; yi þ k1 h þ k2 h 10 40 40
ð50cÞ
3 3 9 6 k4 ¼ f xi þ h; yi þ k1 h k2 h þ k3 h 5 10 10 5
ð50dÞ
11 5 70 35 k5 ¼ f xi þ h; yi k1 h þ k2 h k3 h þ k4 h 54 2 27 27
ð50eÞ
7 1631 175 575 k1 h þ k2 h þ k3 h k6 ¼ f xi þ h; yi þ 8 55; 296 512 13; 824 44; 275 253 k4 hþ k5 h þ 110; 592 4096
ð50fÞ
and h is the step size. First, the slopes for all variables at the initial value (xi) are obtained. These slopes (a set of k1's) are then used to make predictions of dependent variable at the point of xi þh/5. These point values are in turn used to compute a set of slopes at that point (the k2's). A set of k3–k6 is continuously obtained at different points between xi and xi þ 1. Finally, the k's are combined to make the final prediction in Eqs. (48) and (49). If the elastic–softening or softening–residual interface is located between xi and xi þ 1, the large local truncation error can be estimated. In order to adjust the step size, a certain criterion can be used, for example the R–K driver with adaptive steps (subroutine odeint, available in [17]). However, if the system of
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
equation involves rapidly changing components together with slowly changing ones, the stiffness problem can arise in the solution. It might suppose that the adaptive step-size routine might offer a solution for this problem. However, this is not the case, because the stability requirement will still necessitate using very small steps throughout the entire solution [16]. The error in the softening region can be accumulated in the calculation of the residual region. (3) In the case of very small value of η* (that is, the elastic–brittle plastic behavior), the numerical instability can arise because initial values of σ~ θ ð1Þ and υ~ ð1Þ, such as Eqs. (44) and (45), are not valid in the elastic–brittle plastic model. In order to overcome those problems, in this study the classical fourth-order R–K method is used without any adaptive step size yi þ 1 ¼ yi þ 16 ðk1 þ 2k2 þ2k3 þ k4 Þh k1 ¼ f ðxi ; yi Þ
ð52aÞ
1 1 k2 ¼ f xi þ h; yi þ k1 h 2 2
ð52bÞ
1 1 k3 ¼ f xi þ h; yi þ k2 h 2 2
ð52cÞ
k4 ¼ f ðxi þ h; yi þ k3 hÞ
ð52dÞ
The subroutine for the classical fourth-order R–K method is further modified to consider the case of σr r 0 in the residual region. First the prediction of yi þ 1 is obtained at xi using Euler's method, ð53Þ
If σr(i þ 1) r0, then the values of ρi þ 1 and ui þ 1 are obtained by the interpolation. In order to consider the initial values for very small value of η*, the first prediction of η(2) is obtained from the prediction of y(2) using Euler's method. If η(2) Z η*, then the following initial values for the elastic–brittle plastic model, instead of Eqs. (44) and (45), are used: a σ m p σ~ ð1Þ σ~ θ ð1Þ ¼ σ~ r ð1Þ þ ci r 1 r þ sr ð54Þ p1 σ ci ~ þ 2GA~ r ðσ~ r ð1Þ p~ o Þ þ 2GA~ θ ðσ~ θ ð1Þ p~ o Þ υ~ ð1Þ ¼ Ao uð1Þ
ð55Þ
The use of Eq. (55) was verified in [18].
5. Numerical stepwise procedure The simple numerical stepwise procedure, proposed in the previous studies [1,10], is modified to obtain the solution for a spherical opening. In the same way, the dimensional parameter k is used as 1 for a circular opening and 2 for a spherical opening. By dividing the plastic region with a number of thin annular rings (Fig. 3), one can obtain the radial stress at rj for the H–B yield criterion [1,10] qffiffiffiffiffiffiffiffiffiffiffiffiffi σ r ðjÞ ¼ b b2 a ð56Þ where a ¼ σ 2rðj 1Þ 4ℓk
2 1 2m c
b ¼ σ r ðj 1Þ þ ℓk mσ c 2
ℓ¼
r j 1 r j r j 1 þr j
σ σ rðj 1Þ þ sσ 2c
Fig. 3. Typical annulus in plastic region.
ð51Þ
yi þ 1 ¼ yi þ f ðxi ; yi Þh
155
ð57Þ ð58Þ
2 ð59Þ
Table 1 Data set (H–B criterion).
Young's modulus, E (MPa) Poisson's ratio, ν Initial stress, po (MPa) Radius of tunnel, ri (m) σci (MPa) mp sp mr sr ψ peak (deg) ψ r (deg) η
Set 1 [1]
Set 2 [6]
1380 0.25 3.31 5.35 27.6 0.5 0.001 0.1 0 19.47 5.22 0.004742
1380 0.25 3.31 3.0 3.31 0.5 0.001 0.1 0.00001 15.0 5.0 0.0125
using the radii given by 2εθðj 1Þ εrðj 1Þ εrðjÞ rj ¼ rj 1 2εθðjÞ εrðj 1Þ εrðjÞ
ð60Þ
Let the circumferential strain increment be dεθ ¼ εθðjÞ εθðj 1Þ
ð61Þ
Then, the radial strain increment at rj can be obtained as dεrðjÞ ¼ dεerðj 1Þ βðdεθðjÞ dεeθðj 1Þ Þ
ð62Þ
and dεθðj 1Þ are the radial and the circumferential where dε elastic strain increments at rj 1, respectively, and β ¼ ð1 þ
sin ψ Þ= 1 sin ψ . By selecting an arbitrary small value for dεθ, the values of εθ(j), εr(j), rj and uj can be calculated using Eqs. (61), (62) and (60) in that order. The process is then repeated several times to determine completely the strains and displacements in the plastic region.The radial stress at the elastic–plastic interface, in which r¼r1 ¼ re, is given by e rðj 1Þ
e
σ re ¼ po Mσ c
ð63Þ
in which
8sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9 2 < 2 = m p 1 k 1 þ k p o þ sp mp M¼ m2p þ 4 ; 2 1þk : k σc
ð64Þ
Using Eq. (63) as the starting point, successive values of σr(j) can be calculated from Eq. (56) for the radii determined from Eq. (60). The softening parameter at rj can be approximately estimated as
η ¼ εθðj 1Þ εrðj 1Þ þ If η exceeds condition.
1 ðM 11 M 21 Þσ rðj 1Þ ðM 12 M 22 Þσ θðj 1Þ 2G ð65Þ
η*, the strength of the rock mass reaches residual
156
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
The proposed procedure can be easily implemented into the computer program to obtain the ground response curve, by simply modifying the program in reference [10].
Table 2 Data Set 3 (M–C criterion). Radius of tunnel, ri Initial stress, po Young's modulus, E Poisson's ratio, ν η*
3.0 m 20 MPa 10 GPa 0.25 0.008
cp ϕp cr ϕr ψp ¼ ψr
1.0 MPa 301 0.7 MPa 221 3.751
6. Application In order to show the applicability and accuracy of the proposed solutions, the results of the radial displacements, the variations of re/ri and rs/ri, (which can show the spreading of the plastic radius with decreasing pressures) and the radial and circumferential stresses in the plastic region are compared using three data sets. Tables 1 and 2 show the data set for H–B and M–C criteria, respectively. The step size, h¼0.01, is used for similarity solution. Figs. 4–6 show the results of the radial displacements at the opening surface and the variations of re/ri and rs/ri, and the radial and circumferential stresses within the plastic region for each data set. In the figures, SS indicates the results from similarity solution,
Fig. 4. Results for data Set 1. (a) Ground response curve. (b) Variation of re/ri. (c) Variation of rs/ri. (d) Variation of radial stress for pi ¼ 0. (e) Variation of circumferential stress for pi ¼0.
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
while FD indicates those from numerical stepwise procedure. From these figures, the following observations can be made: (1) Good agreement between results from similarity solution and numerical stepwise procedure can be seen for data Set 2 and Set 3 in both circular and spherical openings. (2) For data Set 1, good agreement of the ground response curve between two solutions is obtained. It is noticed that the result of ground response curve from similarity solution in this study is better than that of the previous study (Fig. 4 in [10]). However, the differences in the variations of rs/ri and circumferential stress for pi ¼ 0 from two solutions are shown in Fig. 4 (c) and (e) for both circular and spherical openings. As can be seen in Fig. 4(d) and (e), the variation of radial stresses in the softening region is soft, while the circumferential stresses in
157
the softening region is suddenly changed. The results of similarity solution are obtained using the initial values of Eqs. (54) and (55), while the numerical stepwise procedure shows the elastic-strain softening behavior. (3) Considering the data sets used in this study, the displacements induced by opening from the circular solution are about 4.3– 4.7 times larger than those from the spherical solution.
6.1. Effect of η* As mentioned in Section 4 and shown in Fig. 4(e), the numerical instability can arise in the similarity solution because of very small value of η* or sharp change of circumferential stress in the softening region. The effect of η* is investigated using the
Fig. 5. Results for data Set 2. (a) Ground response curve. (b) Variation of re/ri. (c) Variation of rs/ri. (d) Variation of radial stress for pi ¼ 0. (e) Variation of circumferential stress for pi ¼ 0.
158
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
Fig. 6. Results for data Set 3. (a) Ground response curve. (b) Variation of re/ri. (c) Variation of rs/ri. (d) Variation of radial stress for pi ¼ 0. (e) Variation of circumferential stress for pi ¼0.
data Set 1. Fig. 7 shows the results of ground response curve and circumferential stress for pi ¼0 with different values of η*. As η* decreases, the displacement at the opening increases and the slope of circumferential stress in the softening region is sharply changed. When η* ¼0.005, the results are obtained using the initial values of Eqs. (54) and (55).
7. Conclusions The similarity solution and numerical stepwise procedure for stresses and displacement have been presented for circular and spherical openings excavated in elastic-strain softening rock
masses compatible with M–C and H–B yield criteria. The accuracy and practical application of the proposed methods are illustrated by solving some examples. Although good agreement from similarity solution and numerical stepwise procedure can be obtained for the M–C yield criterion, the numerical instability can arise for the H–B yield criterion in the similarity solution. In order to overcome those instability problems, some measures are undertaken as follows: (1) The classical fourthorder R–K method without any adaptive step size is used to solve the system of ordinary differential equations. (2) Euler's method is used to consider the case of σr r0 in the residual region. (3) The initial values of σ~ θ ð1Þ and υ~ ð1Þ, such as Eqs. (54) and (55), are used to consider the case of very small value of η*.
K.H. Park / International Journal of Rock Mechanics & Mining Sciences 71 (2014) 151–159
159
Fig. 7. Effect of η* (data Set 1). (a) Ground response curve. (b) Variation of circumferential stress for pi ¼ 0.
Although the methods proposed in this paper are limited in scope, it appears to be useful for the preliminary design of a spherical or circular rock opening to obtain the ground response curve. The methods can also be used for the validation of numerical methods, such as the finite element method. References [1] Brown ET, Bray JW, Ladanyi B, Hoek E. Ground response curves for rock tunnels. J Geotech Eng ASCE 1983;109:15–39. [2] Detournay E. Elastoplastic model of a deep tunnel for a rock with variable dilatancy. Rock Mech Rock Eng 1986;19:99–108. [3] Wang Y. Ground response of circular tunnel in poorly consolidated rock. J Geotech Eng ASCE 1996;122:703–1008. [4] Carranza-Torres C, Fairhurst C. The elasto-plastic response of underground excavations in rock masses that satisfy the Hoek–Brown failure criterion. Int J Rock Mech Min Sci 1999;36:777–809. [5] Sharan SK. Elastic–brittle–plastic analysis of circular openings in Hoek–Brown media. Int J Rock Mech Min Sci 2003;40:817–24. [6] Alonso E, Alejano LR, Varas F, Fdez-Manin G, Carranza-Torres C. Ground response curves for rock masses exhibiting strain-softening behavior. Int J Numer Anal Methods Geomech 2003;27:1153–85. [7] Carranza-Torres C. Elasto-plastic solution of tunnel problems using the generalized form of the Hoek–Brown failure criterion. Int J Rock Mech Min Sci 2004;41:480–1.
[8] Sharan SK. Exact and approximate solutions for displacements around circular openings in elastic–brittle–plastic Hoek–Brown rock. Int J Rock Mech Min Sci 2005;42:542–9. [9] Park KH, Kim YJ. Analytical solution for a circular opening in an elastic–brittle– plastic rock. Int J Rock Mech Min Sci 2006;43:616–22. [10] Park KH, Tontavanich B, Lee JG. A simple procedure for ground response curve of circular tunnel in elastic-strain softening rock masses. Tunn Underg Space Technol 2008;23:151–9. [11] Lee YK, Pietruszczak S. A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass. Tunn Underg Space Technol 2008;23:588–99. [12] 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. [13] Wang SL, Yin SD. A closed-form solution for a spherical cavity in the elastic– brittle–plastic medium. Tunn Underg Space Technol 2011;26:236–41. [14] Wang SL, Zheng H, Li CG, Ge XR. A finite element implementation of strainsoftening rock mass. Int J Rock Mech Min Sci 2011;48:67–76. [15] Wang SL, Yin SD, Wu ZJ. Strain-softening analysis of a spherical cavity. Int J Numer Anal Methods Geomech 2012;36:182–202. [16] Chapra SC, Canale RP. Numerical methods for engineers. 4th ed.. New York: McGraw-Hill; 2002. [17] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran 77. Cambridge: Cambridge University Press; 1997. [18] Park KH. Similarity solution of a spherical or circular opening in elastic–brittle plastic rock mass. 2014. [in preparation].