Reliability analysis of underground excavation in elastic-strain-softening rock mass

Tunnelling and Underground Space Technology 60 (2016) 66–79

Contents lists available at ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Reliability analysis of underground excavation in elastic-strain-softening rock mass Li Song a, Hang-Zhou Li a,⇑, Chin Loong Chan b, Bak Kong Low b a b

Department of Civil Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore

a r t i c l e

i n f o

Article history: Received 19 January 2015 Received in revised form 12 April 2016 Accepted 22 June 2016

Keywords: Reliability analysis Elastic-strain-softening Tunnel Reliability-based design FORM SORM

a b s t r a c t This paper deals with the reliability analysis of a circular tunnel in elastic-strain-softening rock mass. Dilatancy angle which varies with softening parameter in different stress conditions is accounted for. Deterministic and probabilistic analyses of the circular tunnel in elastic-strain-softening rock mass are performed. Computational procedures for the first-order and second-order reliability methods (FORM/ SORM) are used in the reliability analyses of the elastic-strain-softening model. The results are in good agreement with those from Monte Carlo simulations incorporating importance sampling. Reliabilitybased design of the required support pressure for the circular tunnel is efficiently conducted. The effect of positive correlation between compressive strength and elastic modulus of the rock mass on the reliability of the tunnel is discussed. The influence of in situ field stress and support pressure as random variables on the probability of failure of the tunnel is investigated. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Convergence confinement method (CCM) is a widely used method in the design of support for underground excavations in rock masses. The method was initially developed in 1930s and refined by other researchers (e.g. Hoek and Brown, 1980; Brown et al., 1983; Hoek et al., 1995), as reviewed by Carranza-Torres and Fairhurst (2000). The method consists of three basic graphs: longitudinal deformation profile (LDP), support characteristic curve (SCC) and ground reaction curve (GRC). This research focuses on the GRC. Analytical solutions are often not available due to the complexity of engineering conditions. Circular tunnels are special cases which have analytical solutions (Ogawa and Lo, 1987; DuncanFama, 1993; Wang, 1996; Carranza-Torres and Fairhurst, 1999; Sharan, 2003, 2008; Park et al., 2008) and can be used for preliminary analysis of underground excavations and to guide design. It is more logical to regard the properties or input data of geomaterials as random variables rather than constant values because of their uncertainties. Reliability analyses that considered uncertainties in rock properties have been conducted by Mollon et al. (2009), Li and Low (2010), among others, using elasticperfectly-plastic model. Further investigations are conducted in ⇑ Corresponding author. E-mail address: [email protected] (H.-Z. Li). http://dx.doi.org/10.1016/j.tust.2016.06.015 0886-7798/Ó 2016 Elsevier Ltd. All rights reserved.

the present study involving more complex constitutive models, Lü and Low (2011), Lü et al. (2011, 2013) using elastic-perfectlyplastic model with first-order reliability method (FORM), secondorder reliability method (SORM) or response surface method (RSM). Zeng and Jimenez (2014) applied a method using FORM and SORM to evaluate the reliability of series geotechnical systems (a layered soil and a circular rock tunnel in a Hoek-Brown rock mass). A method proposed by Zhao et al. (2014) using Least squares support vector machines (LS-SVM) based RSM combined with FORM is applied in tunnel reliability analysis of elasticperfectly-plastic rock mass. Cho (2013) studied the reliability of a clayey soil slope considering multiple failure modes. Hoek and Brown (1997) suggested elastic-brittle-plastic, elastic-strain-softening and elastic-perfectly-plastic behaviors for very good quality hard rock masses, average quality rock masses and very poor quality soft rock masses, respectively. In fact most studies focus on the elastic-perfectly-plastic constitutive model for the theoretical and numerical analysis convenience. However, many rock masses belong to the average classification that is elastic-strain-softening model. Low et al. (2011) and Lü et al. (2013) discussed about the system reliability. Most studies at present mainly considered the unsatisfactory performance of the tunnel as individual failure modes. Due to the correlation of the failure modes and their sharing of common uncertain variables, the component failures are usually somewhat correlated and of different relative

67

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

2400 2200

βσ σ

2000

σσ

x2F

β=

R r

1400

μE σE

r

1200

βσ E

R 1000 Design Point

800

μσ

Failure domain

Limit State Surface (LSS) Limiting up/r i =1%

c

400 10

The matrix formulation of the Hasofer-Lind index for correlated normals is (e.g. Ditlevsen, 1981):

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ¼ min ðx  lÞT C1 ðx  lÞ

c

1600

600

2. Hasofer-Lind index and FORM algorithm based on an intuitive perspective

c

Safe domain

1800

E (MPa)

importance to the system failure. A yet more rational approach is to carry out a tunnel support design to comply with a target system probability of failure. Details can be found in Section 8 in this paper. In this paper, A VBA procedure of the elastic-strain-softening constitutive model is first created in Microsoft Excel to perform the iterative process and the reliability analysis of a circular tunnel in elastic-strain-softening rock mass subjected to a hydrostatic in situ stress field is conducted. Reliability index and failure probability with respect to plastic-zone radius and radial displacement of the circular tunnel are calculated. Probability density functions of plastic-zone radius and radial displacement are obtained using cubic spline interpolation method. System reliability-based design of support pressure to achieve an overall target reliability index is carried out. Other factors influencing the reliability analysis and design of the circular tunnel, such as the support pressure, the in situ field stress and the positive correlation between the compressive strength and the elastic modulus, are also investigated.

ð1Þ

15

20

25

30

35

40

45

50

σc (MPa) Fig. 1. Design point and normal dispersion ellipsoids illustrated in the original random variables’ space.

where x is a vector representing the set of random variables xi , l is the vector of mean values; C is the covariance matrix; and F is the failure domain. According to Eq. (1), the Hasofer-Lind index can be regarded as the minimum distance in units of directional standard deviations from the mean-value point of the random variables to the boundary of the limit state surface. An equivalent formulation (Low and Tang, 1997, 2004) for Eq. (1) is:

and the constrained optimization approach still apply in the original coordinate system, except that the non-normal distributions are replaced by an equivalent normal hyper-ellipsoid, centered not at the original mean of the non-normal distributions, but at the equivalent normal mean lN . The extension of the Hasofer–Lind index to correlated non-normals is known as the first-order reliability method (FORM). Eq. (3) can be rewritten as follows (Low and Tang, 2007):

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  T   xi  l i xi  li b ¼ min R1

b ¼ minx2F

x2F

ri

ri

ð2Þ

where R is the correlation matrix; and ri is the standard deviation of random variable xi . Eq. (2) was preferred by Low and Tang (1997) rather than Eq. (1) because the correlation matrix R is easier to set up than the covariance matrix C, and conveys the correlation structure more explicitly. For correlated non-normals, Eq. (2) can be rewritten as (Low and Tang, 2004):

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  T   N xi  lNi 1 xi  li b ¼ min R N N x2F

ri

ri

ð3Þ

where lNi and rNi are the equivalent normal mean and equivalent normal standard deviation of random variable xi , respectively. The values of lNi and rNi can be computed using the Rackwitz and Fiessler (1978) two-parameter equivalent normal transformation:

/fU1 ½FðxÞg f ðxÞ

ð4Þ

lN ¼ x  rN  U1 ½FðxÞ

ð5Þ

rN ¼

where x is the original non-normal variate, U1 ½ is the inverse of the standard normal cumulative distribution function (CDF), FðxÞ is the original non-normal CDF evaluated at x, /fg is the probability density functions (PDF) of the standard normal distribution, and f ðxÞ is the original non-normal probability density ordinate at x. For correlated non-normals, the ellipsoidal perspective (Fig. 1)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðnÞT R1 ðnÞ

ð6Þ

where n is a column vector of ni and ni ¼ ðxi  lNi Þ=rNi . When the value of ni is varied (automatically) during constrained optimization, the corresponding value of xi is automatically calculated as:

xi ¼ F 1 ½Uðni Þ

ð7Þ

The Low and Tang (2007) algorithm for FORM calculates the reliability index of Eq. (6) using Microsoft Excel’s built-in optimization routine Solver, subject to the constraint that the performance function gðxÞ ¼ 0 (where the x values are program-calculated from Eq. (7)), and by automatically changing the values of ni . Based on the reliability index, the probability of failure can be evaluated from:

pf  1  UðbÞ

ð8Þ

where UðÞ is the cumulative distribution function of the standard normal variable. 3. Problem description In this paper first a VBA procedure is created and verified in Microsoft Excel to perform the iterative process of a circular tunnel initially subjected to a hydrostatic in situ stress and then reliability analysis based on this procedure is carried out. The problem is as follows. When an underground opening is excavated in a stressed rock mass, the stresses in the vicinity of the new opening are

68

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

re-distributed (Hoek, 2007). Prediction of the response around excavation boundary of an opening is a very important component required in the design of the support. Although closed form solutions are only available for simple cases, they still have great value for conceptual understanding of the behavior of the excavations and for the testing and calibration of numerical models. For design purposes, however, these models are restricted to very simple geometries and material models. A circular tunnel initially subjected to a hydrostatic in situ stress field is analyzed, as shown in Fig. 2. The strain-softening model is used to describe the behavior of the rock mass which obeys the Hoek-Brown criterion. The model accounts for three different zones around the tunnel: (i) An elastic zone remote from the tunnel; (ii) an intermediate plastic zone in which the stresses and strains fall on the strain-softening portion; (iii) an inner plastic zone in which stresses and strains are limited by residual strength of the rock masses. The elastic-strain-softening model is shown in Fig. 3. It is required to describe the stresses and displacements in the plastic region to obtain the ground response curve. For simplicity, the Hoek-Brown criterion is used in this study instead of the generalized Hoek-Brown criterion. The parameter a is taken as 0.5. Using the symmetrical characteristic of the circular tunnel, the nonlinear Hoek-Brown criterion can be expressed as a function of radial and circumferential stresses rr and rh (Park et al., 2008):

rh rh rh

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ rr þ mp rc rr þ sp r2c for peak strength qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  rc rr þ sr2c for strain-softening strength ¼ rr þ m qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ rr þ mr rc rr þ sr r2c for residual strength

σ1 − σ 3 0 < η <η *

Strain Sof tening

η =η*

η =0

ε1 Elastic

Softening

Residual

− ε3

f 1

ð9Þ

− ε 3p

h

ð10Þ

1

ð11Þ

where rc is the uniaxial compressive strength of the intact rock, m  s, mr , sr are the and s are Hoek-Brown parameters, mp , sp , m, corresponding peak, strain-softening and residual values of the parameters m and s, respectively.

Fig. 2. A circular tunnel in a hydrostatic in situ stress field (ri is the inner radius of the tunnel, rs is the interface radius of strain softening region and residual region, rp is the interface radius of the elastic region and the strain softening region). The annular region rs < r < rp corresponds to the strain softening zone of Fig. 3, and the annular region ri < r < rs corresponds to the residual zone of Fig. 3.

ε 1e

− ε1p

ε1 αε 1e

Fig. 3. Elastic-strain-softening model (Brown et al., 1983).

The dilatancy angle is assumed to be constant and related to the friction angle / of the rock mass (Hoek and Brown, 1997), with w ¼ /=4, /=8 and 0 for very good quality hard rock, average quality rock, and very poor quality rock, respectively. The strain-softening model shown in Fig. 3 has a constant value of h, which is the gradient of ep1  ep3 line. However, Detournay (1986) noted some errors when using a constant dilatancy angle in elastic-perfectly-plastic Mohr-Coulomb medium and suggested that a variable dilatancy angle should be used instead. Alejano and Alonso (2005) also proposed the use of variable dilatancy by introducing a dilatancy factor. Generally, a strain-softening parameter is used to describe different parts of the stress-strain curves and can be defined in two different ways. According to the first method (Duncan-Fama et al., 1995) it is defined as a function of internal variables, whereas for the second method (Vermeer, 1997) it is defined in an incremental way. If a strain-softening model is defined, the softening parameter g will vary from zero to its peak value g corresponding to the three different stress-strain stages: elastic, strain-softening, and plastic stages, as shown in Fig. 3. In this study, the softening parameter g is used to consider the variation of the dilatancy angle and is defined as the major principal plastic strain (Duncan-Fama et al., 1995; Park et al., 2008). The following expression for g is from Case 3 in Park et al. (2008):

g ¼ eh  eeh

ð12Þ

where superscript e indicates elastic strain. The softening parameter g varies from 0 at the peak strength to g at the residual strength:

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

g ¼ eh  eeh ¼ aeeh  eeh ¼ ða  1Þeeh

ð13Þ

The variable dilatancy within the plastic region can be assumed  vary linearly from the peak values mp , sp and  s and w such that m, wp at g ¼ 0 to the residual values mr , sr and wr at g ¼ g , as shown in Appendix A, Eqs. (A1)–(A7). The equilibrium equation and strain displacement relation are given in Eqs. (A8) and (A9) according to plane strain state and axial symmetry around the circular tunnel (Brown et al., 1983). Because of algebraic complexity, it is impossible to obtain closed form solutions to the complete stress and strain distributions in this case (Brown et al., 1983). Numerical solutions can be obtained using a simple stepwise procedure that successively determines the stresses and strains on the boundaries of a number of annular rings into which the plastic zone is divided. The stepwise procedure begins at the interface of the elastic region and the plastic region with the radius r ¼ r 1 ¼ r p and stops at the internal face of the opening where r ¼ r i (Fig. 2). Details can be found in Eqs. (A10)–(A21). A VBA procedure is created in Microsoft Excel to perform the iterative process based on the sequence of calculation set out in the appendix of Park et al. (2008). In this study, two additional Eqs. (A22) and (A23) are included in the numerical procedure and the expression for g is replaced by Eq. (A21). In order to realize the stepwise calculation, a VBA procedure is created in Microsoft Excel to perform the iterative process based on the sequence of calculation set out in the appendix of Park et al. (2008). In this study, two additional Eqs. (A22) and (A23) are included in the numerical procedure and the expression for g is replaced by Eq. (A21).

4. Verification of the numerical procedure based on the strainsoftening model In order to verify the accuracy of the VBA iterative procedure, comparisons are made among three constitutive models. Constitutive model 1 is the strain-softening model, constitutive model 2 is the elastic-brittle-plastic model and constitutive model 3 is the elastic-plastic model. Solutions obtained by the numerical procedure developed for the strain-softening model in this study are compared with results reported in Brown et al. (1983). The parameters used in this deterministic analysis are from Brown et al. (1983): compressive strength rc ¼ 27:6 MPa, Young’s modulus E ¼ 1380 MPa, mp ¼ 0:5, mr ¼ 0:1, sp ¼ 0:001, sr ¼ 0, a ¼ 3:5, wp ¼ 19:47 calculated from the gradient h ¼ 2:0 and equation h ¼ ð1 þ sin wp Þ=ð1  sin wp Þ, wr ¼ 5:22 from f ¼ 1:2 and f ¼ ð1 þ sin wr Þ=ð1  sin wr Þ. The softening parameter at residual stage g can be obtained from Eq. (13) with the value g ¼ 0:004742, in which eeh can be calculated from Eq. (A18) at the first step of the iteration where rrð1Þ ¼ rre can be given by Eq. (A19), rhð1Þ ¼ rhe by Eq. (9) and shear strength G ¼ E=½2ð1 þ mÞ. Fig. 4 shows that the results obtained in this study are in good agreement with those from Brown et al. (1983). Given the same support pressure pi, rock masses with strain-softening behavior undergo larger displacement ui than elastic-plastic rock mass but smaller displacement than elastic-brittle-plastic rock mass. The elastic-brittle-plastic rock mass deforms significantly especially when the support pressure pi is low. Similarly, for a given support pressure pi, the extent of the plastic zone (ratio rp/ri) of the strainsoftening rock mass is larger than that of the elastic-plastic rock mass but smaller than that of the elastic-brittle-plastic rock mass. Although strain-softening rock masses are commonly encountered in practice, the deterministic and probabilistic studies are mainly focus on elastic-perfectly-plastic constitutive model for convenience. Thus the present study on the strain-softening

69

behavior of rock masses has its practical meaning and can provide some guidance in the design of tunnels and caverns. The verified iterative VBA procedure is used in the following reliability analysis. 5. Performance functions of the circular tunnel In general there are two kinds of failure mechanisms: structurally-controlled failure and stress-controlled failure. In this paper, only stress-controlled failure is considered. The induced stresses due to excavation can cause the failure of the tunnels or caverns if they exceed a certain level. Four parameters, namely compressive strength rc , Young’s modulus E, Hoek-Brown parameters mp and wp are regarded as normal (Gaussian) random variables with coefficient of variation equal to 15%, as shown in Table 1. Positive correlation between rc and E is considered and the correlation coefficient is assumed to be 0.5. The reliability procedure can deal with correlated nonnormal distributions, if desired. Two performance functions are investigated with respect to radius of plastic zone rp and radial displacement ui at the excavated boundary, as follows:

rp ri ui g 2 ðxÞ ¼ eL  ri

g 1 ðxÞ ¼ L 

ð14Þ ð15Þ

where L and eL are the limiting (i.e. permissible) values of r p =r i and ui =r i , respectively, and r i is the internal radius of the circular tunnel. When the performance function becomes negative, it means that a failure event has occurred. When Eq. (14) is negative, it means the extent of the plastic zone (rp/ri) exceeds the limiting (acceptable) value. When Eq. (15) becomes negative, the radial displacement ui at the excavated boundary (where r = ri, Fig. 2) is excessive. 6. Reliability analyses, results and discussions The reliability indices are calculated using the FORM algorithm of Low and Tang (2007) together with Microsoft Excel’s built-in constrained optimization routine Solver as shown in Fig. 5. Initially the column labeled ‘n⁄’ contain zeros. Solver is invoked to find the reliability index b by Setting Target Cell ‘b’ to minimum, via Changing Cells ‘n⁄’, Subject to the Constraint Cell ‘g1(x⁄)’ equal to zero. The original random variables x (on which the performance functions are formulated) are computed automatically from the probabilistic connections between n and x. The solution (b1 = 3.674) obtained by Solver is shown in Fig. 5 for the plastic zone performance function of Eq. (14). The solution for the radial displacement performance function of Eq. (15) (b2 = 1.752) is also shown in the same figure. Both solutions are for the case pi/p0 = 0.1. The corresponding probability of failure approximated by Eq. (8) are pf1 = 0.012% and pf2 = 3.99%, respectively. The solution in the column labeled ‘x⁄’ in Fig. 5 denotes the design point. The design point is the point on the limit state surface (LSS) where the performance function is zero. The LSS separates the safe domain from the unsafe domain in the space of the random variables. As illustrated in Fig. 1 for the two-randomvariable case, the design point is the point of tangency of the expanding ellipse (or the four-dimensional hyperellipsoid of the present case) with the limit state surface and is automatically searched during the optimization process. The reliability index b is the axis ratio (R/r, Fig. 1) of the ellipsoid touching the limit state surface and the one-standard-deviation dispersion ellipsoid. The axis ratio is the same along any radial directions for two concentric ellipsoids. The column labeled ‘n⁄’ indicates the sensitivity of the parameters. Among the four n⁄ values, the n⁄ value of wp has the

70

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

Fig. 4. Comparison of three constitutive models (elastic-strain-softening, elastic-brittle-plastic and elastic-plastic). (a) Relationship of internal support pressure pi and radius rp of plastic zone. (b) Relationship of internal support pressure and radial displacement ui of the tunnel.

Table 1 Random variables used in reliability analysis. Variables

Distribution

Para 1

Para 2

rc (MPa)

Lognormal Lognormal Lognormal Lognormal

27.6 1380 0.5 0.340

4.14 207 0.075 0.051

E (MPa) mp

wp

lowest absolute value in both sets of solutions. This shows that wp has the least influence on both performance functions of Eqs. (14) and (15). In cases where the analysis is time consuming, the reliability analysis can be simplified by omitting the less important random variables with little error. Monte Carlo simulation (MCS) is performed using the commercial software @Risk to compare the actual pf with the FORM pf as obtained above. The coefficient of variation (COV) for the MCS pf can be estimated from the following equation (Shooman 1968):

COV ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffi 1  pf Npf

ð16Þ

where N is the total number of simulation trials. Another use of Shooman’s equation is in estimating the number N needed to

achieve a certain accuracy. With an estimated FORM pf1 of 0.012%, Eq. (16) yields a 5% COV for N = 3 million trials, and a 2% COV for N = 21 million trials. On the other hand, corresponding to an estimated FORM pf2 of 3.99%, Eq. (16) yields a 5% COV for N = 10,000 trials, and a 2% error for N = 60,000 trials. For g1(x), the estimated pf1 is 0.019% from MCS (with 3 million Latin hypercube samples using @RISK program). For g2(x), the estimated pf2 is 4.84% from MCS (with 10,000 Latin hypercube samples using @RISK program). To reduce the computation cost involved in conducting MCS with small pf, an importance sampling technique can be used, such as that based on Melchers (1984). The importance sampling takes place in the rotated space of equivalent standard normal variables u, in the neighborhood of the design point u⁄. The MCS pf below are obtained by repeating five Monte Carlo simulations with importance sampling using the same N number of trials per simulation. For g1(x), the average MCS pf1 is equal to 0.018% with a COV of 7% using N = 1000 trials, and equal to 0.018% with a COV of 3% using N = 10,000 trials. For g2(x), the average MCS pf2 is equal to 4.78% with a COV of 6% using N = 1000 trials, and equal to 4.76% with a COV of 1% using N = 10,000 trials. This shows that importance sampling can improve the efficiency of Monte Carlo simulation significantly.

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

71

Fig. 5. Reliability analysis of circular tunnel with respect to two performance functions.

Yet another practical approach for approximating pf is the second-order reliability method (SORM), such as that based on Der Kiureghian et al. (1987). A paraboloid is constructed at the design point u⁄ in the rotated equivalent standard normal space. Chan and Low (2012) provide a simple procedure for fitting the paraboloid and computing its principal curvatures. The SORM pf is then calculated as a function of the principal curvatures and the FORM b. For the plastic zone performance function or g1(x), the solution for SORM pf1 is obtained as 0.018%. For the radial displacement performance function or g2(x), the solution for SORM pf2 is obtained as 4.74%. For the case in hand, the SORM results are much more accurate than the FORM results, in comparison with those obtained by MCS. Comparisons between FORM, SORM and MCS results are shown in Table 2 for the above case of pi/p0 = 0.1, as well as for

Table 2 Results of probabilistic analysis obtained using different methods. Support pressure ratio pi/p0

Low and Tang (2007)

Chan and Low (2012)

Monte Carlo simulation

FORM b

FORM pf (%)

SORM b

SORM pf (%)

MCS pf (%)

COV (%)

22.91 0.24 0.012 7.1  105 1.6  106

0.688 2.735 3.570 4.646 5.456

24.58 0.31 0.018 1.7  104 2.4  106

24.6 0.3 0.02 1  104 3  106

4 5 7 10 15

1.756 0.583 1.671 3.208 4.342

96.05 27.97 4.74 0.067 7.1  104

96.2 28.4 4.8 0.07 7  104

1 2 6 5 9

(a) Plastic-zone criterion 10–5 0.742 0.05 2.818 0.1 3.674 0.2 4.824 0.3 5.533 (b) Radial-displacement 10–5 1.704 0.05 0.651 0.1 1.752 0.2 3.313 0.3 4.470

criterion 95.58 25.75 3.99 0.046 3.9  104

Note: Five sets of Monte Carlo simulations with importance sampling of size 1000 for each case.

pi/p0 = 105, 0.05, 0.2 and 0.3. In each case, the MCS result is based on 5 simulations with importance sampling using N = 1000 trials. The COV is indicated beside the MCS pf. In all cases the SORM pf estimate is more consistent with the MCS pf, and hence more accurate than the FORM pf estimate. Table 2 also shows that the influence of the internal support pressure pi on the reliability index and the failure probability is significant for both plastic zone and radial displacement criteria. When the excavation is virtually unsupported, the probability that the plastic zone radius exceeds three times the internal radius of the tunnel is about 25% (see Table 2a, case pi/p0 = 105), and the probability that the radial displacement exceeds 1% of the internal radius of the tunnel is about 96% (see Table 2b, case pi/p0 = 105). On the other hand, when a support pressure equivalent to just 10% of the in situ field stress is applied to the excavation boundary, the probability that the tunnel will fail under the plastic zone criterion and the radial displacement criterion will be drastically reduced to about 0.019% and 4.78%, respectively. Therefore, the support pressure plays a very efficient role in restricting the plastic zone and the displacement of the tunnel. Notice that negative reliability indices appear in Table 2b. In using Eqs. (1)–(3) or (6) for computing the reliability index, one needs to distinguish negative from positive indices before using Eq. (8) to estimate the probability of failure. The computed b is positive if the performance function is positive at the mean-value point, i.e. the mean-value point is inside the safe domain. The computed b is negative if the performance function is negative at the mean-value point, i.e. the mean-value point is inside the failure domain.

7. Reliability analysis involving non-normal distributions Lognormal and beta distribution are also adopted here for comparison with the normal distribution. The random variables used in the analyses are summarized in Table 3.

72

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

Table 3 Parameters of lognormal and beta distribution. Variables

Distribution

Para 1

Para 2

rc (MPa)

Lognormal Lognormal Lognormal Lognormal

27.6 1380 0.5 0.340

4.14 207 0.075 0.051

BetaDist BetaDist BetaDist BetaDist

12 12 12 12

12 12 12 12

E (MPa) mp

wp

rc (MPa) E (MPa) mp

wp

Para 3

Para 4

6.9 345 0.125 0.085

48.3 2415 0.875 0.595

Note: Para 1 is mean value of the random variable, para 2 is the standard deviation of the random variable, para 3 and para 4 are the corresponding parameters k1 and k2, respectively.

The mean value l and standard deviation r of a beta distribution with parameters k1 , k2 , min and max are (Evans et al., 1983):

k1 k1 þ k2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi max  min k1 k2 r¼ k1 þ k2 k1 þ k2 þ 1

l ¼ min þðmax  minÞ 

mode ¼ min þðmax  minÞ 

k1  1 ; k1 þ k2  2

ð17Þ ð18Þ k1 > 1; k2 > 1

ð19Þ

The random variables for normal, lognormal and beta distribution are compared in Fig. 6. These three distributions of the four random variables have the same mean value and standard deviation. Good agreement of the probability density functions (PDF)

of the three random variables’ distributions is mainly due to low standard deviations which are assumed to be 15% of the corresponding mean values. Reliability indices for different limiting plastic zone radius rp and different limiting radial displacement ui for the assumed normal, lognormal and beta distributions for the underlying random variables (Tables 1 and 3 and Fig. 6) are calculated for the strainsoftening model of this paper, as shown in Fig. 7 and Table 4. It can be seen from Fig. 7 that the curves obtained based on normally distributed random variables are closer to the curves based on the 4-parameter beta distributions than those obtained from lognormally distributed random variables. This is consistent with the relationship among the three distributions of random variables shown in Fig. 6 where it may be observed that the PDF curves of normally distributed random variables is closer to the PDF curves of beta distributions than those of lognormal distributions. Table 4 shows that the reliability index decreases significantly when the support pressure decreases from pi =p0 ¼ 0:3 to pi =p0 ¼ 105 . In this study, numerical problem will occur for the case pi ¼ 0 when Microsoft built-in routine Solver is invoked. This may be attributed to the fact that the plastic part of the elastic-strain-softening model is discontinuous at the intersection of the softening part and the residual part, and the two points used in computing numerical derivatives may span different parts of such a discontinuity, causing numerical errors (Fylstra et al., 1998). It is found that when the ratio of the support pressure to in situ field stress is below 105 , the plastic zone radius and the radial displacement are already quite close to their final values (when pi =p0 ¼ 1011 ). Therefore in this study a support pressure ratio 105 instead of 0 is used to represent the scenario in which no support is applied.

Fig. 6. Comparison of the probability density functions (PDF) of random variables obeying the normal, lognormal and beta distributions of Tables 1 and 3.

73

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

Fig. 7. Comparison of FORM results for three different probability distributions of the underlying random variables. (a) Relationship between reliability index, internal support pressure pi, and limiting rp/ri. (b) Relationship between reliability index, internal support pressure pi, and limiting ui/ri.

Table 4 Probability of failure computed for 3 different distribution types. Support pressure ratio pi/p0

Probability of failure (%) obtained by assuming: Normal distribution

(a) Plastic-zone criterion 105 0.05 0.1 0.2 0.3

22.91 0.24 0.012 7.1  105 1.6  106

(b) Radial-displacement criterion 105 0.05 0.1 0.2 0.3

95.58 25.75 3.99 0.046 3.9  104

Lognormal distribution

Beta distribution

(24.58) (0.31) (0.018) (1.7  104) (2.4  106)

25.13 0.058 2.6  104 8.7  1010 0

(25.04) (0.058) (2.6  104) (8.7  1010) (0)

23.55 0.22 6.1  103 2.2  106 6.1  1011

(24.89) (0.23) (6.2  103) (1.9  106) (4.2  1011)

(96.05) (27.97) (4.74) (0.067) (7.1  104)

95.33 28.06 3.28 3.5  103 9.6  108

(95.49) (28.71) (3.44) (3.8  103) (1.1  107)

95.30 26.35 4.22 0.031 3.6  105

(96.02) (28.30) (4.69) (0.035) (3.6  105)

Note: Numbers without brackets are based on FORM. Numbers within brackets are based on SORM.

A series of reliability indices (and corresponding cumulative probability (CDF) values) can be obtained for different trial limiting ratios. Applying cubic spline interpolation (Kreyszig, 1993) to the CDF, the probability density function (PDF) curves can then be plotted. Details can be found in Low (2008). The PDF curves thus obtained for the normally distributed random variables of Table 1 are shown in Fig. 8. For comparison, the corresponding PDF curves obtained from Monte Carlo simulation using commercial software @Risk with 100,000 trials are also plotted.

Figs. 9 and 10 compare the probability of failure curves obtained from FORM/SORM and those from Monte Carlo simulations, for lognormal and beta distributions, respectively. In the case of pi/p0 = 105, results obtained from SORM are in better agreement with those from Monte Carlo simulations for certain limiting values of r p =r i or ui =r i for which the probabilities of failure obtained from FORM are slightly smaller than those obtained from Monte Carlo simulations. Theoretically, the probabilities of failure obtained from Monte Carlo will be identical to those inferred from FORM only when

74

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

Fig. 8. Comparison of results between FORM and Monte Carlo simulations, for normally distributed random variables and strain softening model. (a) Probability density function of rp/ri. (b) Probability density function of ui/ri.

the limit state surface is a hyperplane and the random variables are normally distributed and the correlation matrix is an adentity matrix. In this study, the limit state surfaces are not hyperplanes, as can be seen in Fig. 1 in the rc and E space. FORM algorithm only accounts for the corresponding domain behind the hyperplane and ignores the domain between the hyperplane and the limit state surface, resulting in slight inaccuracy for the case in hand. On the other hand, SORM builds on the FORM design point and FORM b value, and accounts to some extent for the curvatures of the limit state surface, leading to more accurate results.

8. Reliability-based design of support pressure It is more rational to design support according to a given probability of non-failure. This means the support pressure needs to be determined accounting for uncertainty. The given reliability index or probability of failure should be set to a target value, e.g. reliability index 2.5 (corresponding to a probability of failure of about 0.621%). The design support pressure which satisfies a target failure probability of 0.62% is shown in Fig. 11a and b as functions of the limiting plastic zone radius and the limiting radial displacement, respectively. Normal, lognormal and beta distributions are

considered in producing the reliability-based design chart in Fig. 11. It can be seen that for both the plastic zone criterion and the radial displacement criterion, the required support pressure decreases with the limiting plastic zone, and with the limiting radial displacement. The required support pressures in Fig. 11 correspond to a failure probability of 0.62% as determined by SORM instead of FORM for the following reasons. It was established in the earlier sections that SORM was more accurate than FORM. Furthermore, FORM was found to underestimate the failure probability for both failure modes in the present example. Therefore, it is more conservative to design the required support pressure based on SORM b indices than FORM b indices for the case in hand. For example, if a support pressure of pi = 0.475 MPa is applied to the excavation boundary, the FORM b is equal to 2.5 (and FORM yields pf = 0.62%) for the radial displacement failure mode. For the same support pressure, SORM yields pf = 0.8% which exceeds the permissible pf. A yet more rational approach is to carry out a tunnel support design to comply with a target system probability of failure. In other words, the tunnel is regarded as a series system, and the performance of the system is unsatisfactory if either plastic zone or radial displacement criterion is violated, or both criteria are simultaneously violated. The probability of failure of a series

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

75

Fig. 9. Comparison of results among FORM, SORM and Monte Carlo simulations, for lognormally distributed random variables and strain softening model. (a) Probability of failure versus limiting rp/ri. (b) Probability of failure versus limiting ui/ri.

system (pf,sys) can be assessed by computing the Kounias-Ditlevsen (Kounias, 1968; Ditlevsen, 1979) bi-modal bounds. The system failure probability bounds are calculated as functions of FORM/SORM b indices and the intermodal correlation coefficient between the two failure modes. The intermodal correlation coefficient is in turn dependent upon the direction vectors (linking the mean-value point and design point) of the two modes. Automatic and efficient implementation of the bi-modal bounds for systems with multiple failure modes is carried out by the simple program code in Low et al. (2011). Fig. 12 shows the support pressure (pi) required to achieve a system target failure probability of 0.62% (i.e. equivalent FORM bsys of 2.5) for the circular tunnel as a function of the limiting rp/ri ratio. The three curves correspond to limiting ui/ri ratios of 1%, 2% and 5%, respectively. Note that the three curves coincide for limiting rp/ri ratios not greater than 1.5. This indicates that, within this range, the value of pi depends only on the limiting rp/ri ratios. Earlier it has been established that a high value of pi is needed to maintain a reliability index of 2.5 for limiting rp/ri ratios lower than 1.5 (see Fig. 11a again). At such a high value of pi, the dominant failure mode is plastic zone failure (Mode 1). For example, when pi = 0.937 MPa, b1 = 2.50 (for limiting rp/ri = 1.25) whereas b2 = 4.17, 4.82 and 5.54 (for limiting ui/ri = 1%, 2% and 5%, respectively). In this case, the system probability of failure is largely

attributed to Mode 1, and therefore not sensitive to the limiting ui/ri ratio. Also note that the three curves in Fig. 12 shows decreasing trends with increasing limiting rp/ri ratios, and eventually approach constant pi values of about 0.5 MPa, 0.3 MPa and 0.1 MPa for limiting ui/ri ratios of 1%, 2% and 5%, respectively. At high limiting rp/ri ratios, the dominant failure mode is radial displacement failure (Mode 2). For example, when limiting rp/ri = 4, b2 = 2.50 (at pi = 0.5 MPa, 0.3 MPa and 0.1 MPa for limiting ui/ri = 1%, 2% and 5%, respectively) whereas b1 = 3.04, 4.05 and 4.68 (at pi = 0.5 MPa, 0.3 MPa and 0.1 MPa, respectively). In this case, the system probability of failure is largely attributed to Mode 2, and therefore not sensitive to the limiting rp/ri ratio. The pi values of 0.5 MPa, 0.3 MPa and 0.1 MPa are essentially the support pressures needed to limit the radial displacement to 1%, 2% and 5% of the original tunnel radius, respectively. The system reliability-based design chart in Fig. 12 is prepared for the case of normally-distributed random variables. If lognormal distribution (or any other non-normal distribution) is used in lieu of normal distribution, the same system reliability design approach can still be applied but resulting in a different design chart. Limiting rp/ri ratios as high as 4 and limiting ui/ri ratios as high as 5% are included in this study to provide insights under certain extreme situations, and should not be taken as design recommendations.

76

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

Fig. 10. Comparison of results among FORM, SORM and Monte Carlo simulations, for random variables obeying beta distributions and strain softening model. (a) Probability of failure versus limiting rp/ri. (b) Probability of failure versus limiting ui/ri.

The designer must specify his limiting ratios based on the project at hand.

9. Other factors influencing the reliability analysis In the above analysis, the compressive strength and elastic modulus of the rock mass is assumed to be positively correlated, with correlation coefficient equal to 0.5. The influence of the correlation between the compressive strength and elastic modulus on the reliability analysis is investigated. Fig. 13 shows the variation of the reliability index and corresponding probability of failure when the correlation coefficient varies from 0 to 0.9, for pi/p0 = 0.1 and normal random variables as shown in Table 1. The lower values of probability of failure correspond to the higher values of reliability index which means that if the correlation between the two random variables is not considered, the reliability will be overestimated. This can be understood intuitively in the random variables space. If random variables of resistance nature are positively correlated, the dispersion hyper-ellipsoid will have its major axis directed towards the limit state surface (Fig. 1) and will expand less to touch the limit state surface than uncorrelated dispersion hyper-ellipsoid. More can be found in Low (2008). Hoek (1998) pointed out that the uncertainty associated with estimating the properties of in situ rock masses has a significant

impact on the design of slopes and excavations in rock. Hoek (2007) also mentioned that the design of support is a very imprecise process, especially if shotcrete is involved. Based on earlier analyses of this study, it can be seen that support pressure has a significant effect on the reliability analysis results. A small change in support pressure may greatly influence the result of the reliability index or the failure probability of the tunnel. It is therefore more rational to account for the uncertainties in in situ field stress and support pressure and to treat them as random variables. Further studies are carried out, considering field stress and support pressure uncertainties. Results are shown in Tables 5a and 5b for plastic zone performance function and radial displacement performance function, respectively. Three scenarios are investigated. Scenario 1 involves the original four random variables i.e. rc, E, mp and wp. Scenario 2 involves five random variables with pi =p0 as the fifth random variable. Scenario 3 involves six random variables with pi =p0 and p0 as additional random variables. In Scenarios 2 and 3, the support pressure ratio pi/p0 is assumed to be normally distributed with mean values of 0.05. 0.1, 0.2 and 0.3 (see Tables 5a and 5b). In Scenario 3, the hydrostatic in situ field stress p0 is assumed to be a normal random variable with a mean value of 3.31 MPa. The coefficients of variation of the above two random variables are both equal to 15%. It can be seen that the probability of failure increases with the variation of in situ field stress and support pressure. Therefore it is recommended that in situ field stress

77

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

Required support pressure (MPa)

(a)

1.2 1 0.8 Normal Lognormal Beta

0.6 0.4 0.2 0

0

1

2

3

4

5

Limiting rp/ri

Required support pressure (MPa)

(b) 0.5 0.4 0.3

Normal Lognormal

0.2

Beta

0.1 0.0

0

1

2

3

4

5

6

Limiting ui/ri(%) Fig. 11. Reliability-based design of the required support pressure with respect to (a) plastic zone performance function and (b) radial displacement performance function, for b = 2.5 and random variables as shown in Tables 1 and 3.

Required support pressure pi (MPa)

1.2 1.0 0.8 0.6

Limiting ui/ri = 1%

0.4

Limiting ui/ri = 2%

0.2

Limiting ui/ri = 5%

0.0 1

2

3

4

Limiting rp/ri Fig. 12. Support pressure required for achieving an overall system reliability index of 2.5 (or system failure probability of 0.62%) for the case of normally-distributed random variables.

Fig. 13. Effect of correlation between compressive strength and elastic modulus on reliability index and probability of failure for pi/p0 = 0.1: (a) failure in plastic zone criterion and (b) failure in radial displacement criterion.

Table 5a Probability of failure computed for 3 different scenarios. (a) Plastic-zone criterion. Support pressure ratio pi/p0

Probability of failure (%) obtained for: Scenario 1

0.05 0.1 0.2 0.3

0.24 0.012 7.1  105 1.6  106

Scenario 2 (0.31) (0.018) (1.7  104) (2.4  106)

0.27 0.017 1.9  104 5.8  106

Scenario 3 (0.35) (0.025) (3.7  104) (1.3  105)

Note: Numbers without brackets are based on FORM. Numbers within brackets are based on SORM. Limiting rp/ri = 3.

0.75 0.062 8.2  104 2.2  105

(0.88) (0.080) (1.3  103) (4.1  105)

78

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79

Table 5b Probability of failure computed for 3 different scenarios. (b) Radial-displacement criterion. Support pressure ratio pi/p0

Probability of failure (%) obtained for:

0.05 0.1 0.2 0.3

25.75 3.99 0.046 3.9  104

Scenario 1

Scenario 2 (27.97) (4.74) (0.067) (7.1  104)

Scenario 3

26.15 4.62 0.11 3.4  103

(28.55) (5.51) (0.16) (5.6  103)

31.80 9.17 0.45 0.019

(33.32) (10.07) (0.55) (0.026)

Note: Numbers without brackets are based on FORM. Numbers within brackets are based on SORM. Limiting ui/ri = 0.01.

and support pressure should be considered as random variable in the reliability analysis and reliability-based design of the tunnels or caverns.

10. Summary and conclusions Due to the complexities of the problems, analytical solutions of underground excavations are commonly restricted to circular shapes with elastic-perfectly-plastic or elastic-brittle-plastic characteristics which are not always consistent with the behavior of practical rock masses. An elastic-strain-softening model is used in this paper. Dilatancy angle which varies with softening parameter in different stress conditions is taken into account. The solution for the elastic-strain-softening model is obtained through numerical iterative calculations. Reliability analysis was illustrated for a circular tunnel based on the elastic strain-softening model. First-order and second-order reliability methods (FORM/SORM) were used. Comparisons were made with Monte Carlo simulations incorporating importance sampling. Reliability-based design of the required support pressure for the underground excavation was demonstrated. The support pressure needed to achieve a target reliability index for individual failure modes, as well as to achieve an overall target reliability index for the system of multiple failure modes was evaluated. The probability of failure of a series system was assessed by computing the Kounias-Ditlevsen bi-modal bounds. The results can guide the analysis and design of actual tunnels. Other factors which influenced the reliability of the tunnel were investigated. Positive correlation between the compressive strength and the elastic modulus should be considered. It was shown that the results of reliability analysis were unconservative if the compressive strength and the elastic modulus were treated as uncorrelated random variables in cases that the failure probability is not very large. The effect of in situ field stress and support pressure on the reliability of the tunnel was also investigated by regarding them as additional random variables. It is recommended that in situ field stress and support pressure should be considered as random variables in reliability analysis and design.

Acknowledgement The authors would like to thank the Jurong Town Corporation of Singapore for providing research funding which enabled their research at Nanyang Technological University. This research was also supported by the Fundamental Research Funds for the Central Universities (xjj2013076, xjj2014127), the Nature Science Foundation of Shaanxi (2016JM5047) and the Open Research Fund of State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, under grant NO. 201105047.

Appendix A. Solution of strain-softening model Eqs. (A1)–(A20) below are as given in Park et al. (2008). Eqs. (A21)–(A23) are new expressions derived in this paper. Parameters  at the strain-softening stage are assumed to vary lin s and w m, early from the peak values mp , sp and wp at g ¼ 0 to the residual values mr , sr and wr at g ¼ g :

 ¼ mp  ðmp  mr Þ m s ¼ sp  ðsp  sr Þ

g g

g g

ðA1Þ ðA2Þ

 ¼ w  ðw  w Þ g w p p r 

g

ðA3Þ

For elastic-perfectly-plastic model, the parameters can be simplified (Brown et al., 1983):

 ¼ mp ¼ mr m s ¼ sp ¼ sr

ðA4Þ ðA5Þ

For elastic-brittle-plastic model, the parameters can be simplified (Brown et al., 1983):

 ¼ mr m s ¼ sr

ðA6Þ ðA7Þ

The equilibrium equation and strain displacement relation can be given as follows according to plane strain state and axial symmetry around the circular tunnel (Brown et al., 1983):

drr rr  rh þ ¼0 dr r du u er ¼  ; eh ¼  dr qffiffiffiffiffiffiffiffiffiffiffiffiffiffi r

rrðjÞ ¼ b  b2  a

ðA8Þ ðA9Þ ðA10Þ

where

  1  rc rrðj1Þ þ sr2c m 2  rc b ¼ rrðj1Þ þ km  2 r j1  r j k¼ r j1 þ r j

a ¼ r2rðj1Þ  4k

ðA11Þ ðA12Þ ðA13Þ

where the radii is given by

2ehðj1Þ  erðj1Þ  erðjÞ rj ¼ rj1 2ehðjÞ  erðj1Þ  erðjÞ

ðA14Þ

The circumferential strain increment can be given as

deh ¼ ehðjÞ  ehðj1Þ

ðA15Þ

The radial strain increment can be obtained as (Park et al., 2008):

derðjÞ ¼ deerðj1Þ  bðdehðjÞ  deehðj1Þ Þ

ðA16Þ

L. Song et al. / Tunnelling and Underground Space Technology 60 (2016) 66–79 

w e e where b ¼ 1þsin  , derðj1Þ and dehðj1Þ are radial and circumferential 1sin w

elastic strain increments at r j1 , respectively. The elastic strains at r j can be obtained as suggested in Park et al. (2008)

 1  ð1  mÞðrrðjÞ  p0 Þ  mðrhðjÞ  p0 Þ 2G  1  ð1  mÞðrhðjÞ  p0 Þ  mðrrðjÞ  p0 Þ ¼ 2G

eerðjÞ ¼

ðA17Þ

eehðjÞ

ðA18Þ

where G is the shear modulus and v the Poisson’s ratio of the rock mass. The radial stress at the interface is given in Brown et al. (1983):

rre ¼ p0  Mrc

ðA19Þ

where (Brown et al., 1983)

M¼

 12 1 mp 2 mp m0 mp þ þ sp  2 4 rc 8

ðA20Þ

The softening parameter g at rj can be calculated as:

g ¼ ehðj1Þ 

rhðj1Þ 2G

ðA21Þ

In the VBA procedure, the following initial values are used for the incremental strains:

Deh1 ¼ Deeh1 ¼ 0:01eh1 Der1 ¼ Deer1 ¼ 0:01er1

ðA22Þ ðA23Þ

References Alejano, L.R., Alonso, E., 2005. Considerations of the dilatancy angle in rocks and rock masses. Int. J. Rock Mech. Min. Sci. 42 (4), 481–507. Brown, E.T., Bray, J.W., Ladanyi, B., et al., 1983. Ground response curves for rock tunnels. J. Geotech. Eng. ASCE 109 (1), 15–39. Carranza-Torres, C., Fairhurst, C., 1999. The elasto-plastic response of underground excavations in rock masses that satisfy the Hoek-Brown failure criterion. Int. J. Rock Mech. Min. Sci. 36, 777–809. Carranza-Torres, C., Fairhurst, C., 2000. Application of convergence-confinement method of tunnel design to rockmasses that satisfy the Hoek-Brown failure criterion. Tunnel. Undergr. Space Technol. 15 (2), 187–213. Chan, C.L., Low, B.K., 2012. Practical second-order reliability analysis applied to foundation engineering. Int. J. Numer. Anal. Methods Geomech. 36 (11), 1387– 1409. Cho, S.E., 2013. First-order reliability analysis of slope considering multiple failure modes. Eng. Geol. 154, 98–105. Detournay, E., 1986. Elastoplastic model of a deep tunnel for a rock with variable dilatancy. Rock Mech. Rock Eng. 19, 99–108. Der Kiureghian, A., Lin, H.Z., Hwang, S.J., 1987. Second-order reliability approximations. J. Eng. Mech. 113 (8), 1208–1225. Ditlevsen, O., 1979. Narrow reliability bounds for structural systems. J. Struct. Mech. 7 (4), 453–472. Ditlevsen, O., 1981. Uncertainty Modeling: With Applications to Multidimensional Civil Engineering Systems. McGraw-Hill, New York. Duncan-Fama, M.E., 1993. Numerical modelling of yield zones in weak rocks. In: Hudson, J.A. (Ed.), Comprehensive Rock Engineering, vol. 2. Pergamon, Oxford, pp. 49–75.

79

Duncan-Fama, M.E., Trueman, R., Craig, M.S., 1995. Two and three dimensional elastoplastic analysis for coal pillar design and its application to highwallmining. Int. J. Rock Mech. Sci. Geomech. Abstr. 32 (3), 215–225. Evans, M., Hastings, N., Peacock, B., 1983. Statistical Distributions, second ed. Wadsworth, Belmont (CA). Fylstra, D., Lasdo, L., Watson, J., et al., 1998. Design and use of the Microsoft Excel Solver. Interfaces 28 (5), 29–55. Hoek, E., 1998. Reliability of Hoek-Brown estimates of rock mass properties and their impact on design. Int. J. Rock Mech. Min. Sci. 35 (1), 63–68. Hoek, E., 2007. Practical Rock Engineering. . Hoek, E., Brown, E.T., 1980. Underground Excavations in Rock. Chapman & Hall, London. Hoek, E., Kaiser, P.K., Bawden, W.F., 1995. Support of Underground Excavations in Hard Rock. Balkema, Rotterdam. Hoek, E., Brown, E.T., 1997. Practical estimates of rock mass strength. Int. J. Rock Mech. Min. Sci. 34 (8), 1165–1186. Kounias, E.G., 1968. Bounds for the probability of a union, with applications. Ann. Math. Stat. 39 (6), 2154–2158. Kreyszig, E., 1993. Advanced Engineering Mathematics, seventh ed. John Willey & Sons, New York. Li, H.Z., Low, B.K., 2010. Reliability analysis of circular tunnel under hydrostatic stress field. Comput. Geotech. 37, 50–58. Low, B.K., Tang, W.H., 1997. Reliability analysis of reinforced embankments on soft ground. Can. Geotech. J. 34 (5), 672–685. Low, B.K., Tang, W.H., 2004. Reliability analysis using object-oriented constrained optimization. Struct. Saf. 26 (1), 69–89. Low, B.K., Tang, W.H., 2007. Efficient spreadsheet algorithm for first-order reliability method. J. Eng. Mech. ASCE 133 (12), 1378–1387. Low, B.K., 2008. Chapter 3: Practical reliability approach using spreadsheet. In: Reliability-Based Design in Geotechnical Engineering: Computations and Applications. Taylor & Francis, Routledge, UK, pp. 134–168. Low, B.K., Zhang, J., Tang, W.H., 2011. Efficient system reliability analysis illustrated for a retaining wall and a soil slope. Comput. Geotech. 38 (2), 196–204. Lü, Q., Chan, C.L., Low, B.K., 2013. System reliability assessment for a rock tunnel with multiple failure modes. Rock Mech. Rock Eng. 46, 821–833. Lü, Q., Low, B.K., 2011. Probabilistic analysis of underground rock excavations using response surface method and SORM. Comput. Geotech. 38, 1008–1021. Lü, Q., Sun, H.Y., Low, B.K., 2011. Reliability analysis of ground-support interaction in circular tunnels using the response surface method. Int. J. Rock Mech. Min. Sci. 48, 1329–1343. Melchers, R.E., 1984. Efficient Monte-Carlo Probability Integration. Research Report No. 7. Department of Civil Engineering, Monash University, Clayton, Australia. Mollon, G., Dias, D., Soubra, A.H., 2009. Probabilistic analysis and design of circular tunnels against face stability. Int. J. Geotech. 9 (6), 237–249. Ogawa, T., Lo, K.Y., 1987. Effects of dilatancy and yield criteria on displacements around tunnels. Can. Geotech. J. 24, 100–113. Park, K.H., Tontavanich, B., Lee, J.G., 2008. A simple procedure for ground response curve of circular tunnel in elastic-strain softening rock masses. Tunn. Undergr. Space Technol. 23, 151–159. Rackwitz, R., Fiessle, B., 1978. Structural reliability under combined random load sequences. Comput. Struct. 9, 484–494. Sharan, S.K., 2003. Elastic-brittle-plastic analysis of circular openings in HoekBrown media. Int. J. Rock Mech. Min. Sci. 40, 817–824. Sharan, S.K., 2008. Analytical solutions for stresses and displacements around a circular opening in a generalized Hoek-Brown rock. Int. J. Rock Mech. Min. Sci. 45, 78–85. Vermeer, P.A., 1997. Non Associated Plasticity for Soils, Concrete and Rock. Physics of Dry Granular Media, Dordrecht, Netherlands, pp. 163–196. Wang, Y., 1996. Ground response of circular tunnel in poorly consolidated rock. J. Geotech. Eng. ASCE 122 (9), 703–708. Zeng, P., Jimenez, R., 2014. An approximation to the reliability of series geotechnical systems using a linearization approach. Comput. Geotech. 62, 304–309. Zhao, H., Ru, Chang.X., Yin, S., Li, S., 2014. Reliability analysis of tunnel using least square support vector machine. Tunn. Undergr. Space Technol. 41, 14–23.