Tunnel Probabilistic Structural Analysis Using the FORM

In this paper tunnel probabilistic structural analysis (TuPSA) was performed using the first order reliability method (FORM). In TuPSA, a tunnel performance function is defined according to the boundary between the structural stability and instability. Then the performance function is transformed from original space into the standard normal variable space to obtain the design point, reliability index, and also the probability of tunnel failure. In this method, it is possible to consider the design factors as the dependent or independent random parameters with arbitrary probability distributions. A software code is developed to perform the tunnel probabilistic structural analysis (TuPSA) using the FORM. For validation and verification of 
 TuPSA, a typical tunnel example with random joints orientations as well as mechanical properties has been studied. The results of TuPSA were compared with those obtained from Monte-Carlo simulation. 
The results show, in spite of deterministic analysis which indicates that the rock blocks are stable, that TuPSA resulted in key-blocks failure with certain probabilities. Comparison between probabilistic and deterministic analyses results indicates that probabilistic results, including the design point and probability of failure, are more rational than deterministic factor of safety.


Introduction
Structural instability is the most dominant failure mechanism of underground openings in moderately jointed rock masses.Traditional analyses in such openings have been largely based on rock mass classification methods.Rock load [1], RQD [2], RSR [3], RMR [4], and rock tunneling quality index () [5] are the most practical rock mass classifications.Despite the empirical effectiveness of these classifications, they largely ignored particular stability problem due to the formation of removable rock blocks around the tunnel walls.
Block instability is not necessarily dependent on any rock mass classifications.Block theory [6] provides a mathematical and geometrical procedure to define the stability of rock structures triggered by discontinuities geometry.For a number of joint intersections behind an excavation free face, block theory determines the combinations of joints half spaces which create a removable rock block, block geometry, probable sliding or falling direction, and the safety factor of block using a rigorous approach.Block theory method considers the rock mass and support properties as deterministic parameters.However, they improved our understanding of structural instability of tunnels, in spite of their assumptions, most of the geomechanical characteristics of rock masses such as discontinuities orientation and mechanical properties and also support properties vary in wide ranges.
Unlike the probabilistic stress-controlled instability analysis, few researches have been performed on probabilistic instability analysis of block failures in underground openings: [7][8][9].On the other hand, in these cases, the analysis is limited to the specific case studies with limited random dimension or specific probability distributions.So far, there is no direct study on the probabilistic analysis of structural instability in underground openings.Therefore, in this research a tunnel probabilistic structural stability analysis (TuPSA) was developed using first order reliability method.
For this purpose, at first, the reliability problem and first order reliability method (FORM) are explained, and then the performance function for tunnel structural stability is presented.By combining these two definitions, tunnel probabilistic structural stability analysis (TuPSA) will be explained.A computer code is developed for reliability analysis of structural failures and performing the TuPSA method.The reliability problem of a typical tunnel example is solved by developing computer code assuming different scenarios.Finally, the results are compared to those obtained from Monte-Carlo simulations.

First Order Reliability Method (FORM)
For a tunnel system, reliability could be defined as the ability of the tunnel or its components to perform the required functions under the specific uncertain conditions.A random uncertain parameters vector such as ( 1 ,  2 , . . .,   ), which could include rock material, support, and joints geometrical parameters, follows a multivariate probability density function ( 1 ,  2 , . . .,   ).For a specific designed tunnel system, there is a performance function ( 1 ,  2 , . . .,   ) that is defined to be zero at a limit state, less than zero at system failure domain, and greater than zero for safe state.The objective of reliability analysis is to calculate the probability of failure which can be expressed as [10,11] where   is failure probability and ( 1 ,  2 , . . .,   ) is uncertain parameters vector.The concepts of performance limit state, failure and safe regions, and probability of failure (domain of integration that is illustrated by shaded region) are illustrated in Figure 1.In Figure 1, shaded area indicates the failure domain and dashed lines are the contours that indicate the probability distributions of two variables (for sake of simplicity).
The general solution to (1) is the Monte-Carlo simulation that relies on repeated random sampling based on the multivariable density function to compute their performance results.The probability of failure is estimated from the ratio of number of failed cases,   , per total number of simulation iterations, , as follows: The sample size () must ensure coefficient of variation of calculated   , (CoV   ), becomes reasonably small.Coefficient of variation of   is estimated using the following equation [12]: The analytical solution of (1) is not available, unless the multivariable density function is normal and the performance function is linear or quadratic [10].These conditions led to the idea of FORM solution.In FORM solution, by replacing the physical nonnormal random variables, ( 1 ,  2 , . . .,   ), with the standard normal random variables (uncorrelated Gaussian random variables with zero mean and unit standard deviation), the performance function is transmitted to standard normal space spanned by standard normal variables ( 1 ,  2 , . . .,   ) (Figure 2).
Most likely failure point in standard space is called the design point or -point ( * in Figure 2).By replacement of the actual limit state ( = 0) with the approximate linear limit state function at design point (  = 0), the probability of failure could be approximated by [11]   ≈ Φ (−) .
Here,  is the distance between origin and design point in normal space and Φ is the cumulative density function (CDF) of standard normal variables.Using this transformation, the multiple integration of (1) reduces to a constrained nonlinear optimization problem as where  = ( 1 ,  2 , . . .,   ) is the vector of standard normal random variables,   is the transpose of ,  = ( 1 ,  2 , . . .,   ) is the physical nonnormal random variables vector, and  is the translation function used to convert the standard normal vector to the physical variables vector.For the noncorrelated and correlated physical variables, translation function can be defined as (6) [13] and (7), respectively: Each nonnormal component,   , can follow any arbitrary cumulative distribution function,   , and  is the Cholesky factor of correlation coefficients matrix of random variables, , [10].Cholesky factor can be determined using the following relation: Here,  is the correlation coefficients matrix of physical random variables and   is the transpose matrix of .The first order reliability method (FORM) could be performed on the systems with known performance function and specific random parameters by using ( 4)-( 8).

Tunnel Probabilistic Structural Stability Analysis (TuPSA)
In TuPSA method, by defining the tunnel performance function, it is possible to perform the FORM analysis using (4)- (8).Definition of tunnel structural performance function is one of the most important stages in TuPSA which is determined by the goals of tunnel design.The goal of geomechanical designs in dealing with structural instabilities is to prevent the blocks and wedges from sliding and falling into the underground opening.Therefore, the performance function is defined as follows: It is difficult to develop an analytical performance function which includes all (most important) design factors, but it is possible by using the nested functions.Details of the performance function calculation, used in this paper, are presented here.The main assumptions of TuPSA are as follows: (i) The blocks formed around the tunnel have maximum possible volume (considering the key blocks).
(ii) Specific falling/sliding blocks have no influence on other blocks' stability.
(iii) Bolts are considered as the simple forces applied to the blocks in specific directions.
(iv) The influence of field stress on blocks stability [14] is not considered.
(v) The distribution of uncertain parameters is available.
Many parameters are influencing the tunnel structural instabilities.The main parameters governing the overall stability of opening and their symbols used in analysis are as follows: (1) geometrical parameters: (a) tunnel axis orientation: trend and plunge, TO, (b) tunnel radius, , (c) joints orientation: joints dip and dip direction, JD, (2) mechanical parameters: (d) cohesion of joints,   , (e) friction angle of joints,   , (f) tensile strength of joints, Ten, (g) water pressure on joints face, , (h) unit weight of rock mass, , (i) unit weight of shotcrete,   , (j) shotcrete thickness, , (k) shotcrete shear strength,   , (l) bolts force pattern around the tunnel,   , (m) seismic coefficient, , (n) seismic force direction, Sd.
There are different consecutive stages to define the tunnel structural performance function.Some of these stages deal with the geometry of blocks and other ones deal with the limit-equilibrium analysis.Different software codes have been developed to compute these functions.By combining these codes using the procedure presented in Figure 3, the performance function for tunnel structural stability has been achieved.Some of these stages, which are used in this research, are presented here briefly and in the functional form.3.1.Removability Function.For 3 joint sets, considering all combinations of joints half spaces, there are eight tetrahedral block types codified 000, 001, 010, 011, 100, 101, 110, and 111.A 0 in the th place signifies that spherical triangle belongs to the upper half-space of th joint plane and conversely a 1 signifies that spherical triangle belongs to the lower half-space of th joint plane.As defined by Figure 3, by checking the kinematic instability of a block combined from joints half spaces, the removability of block is determined (from a geometrical point of view) based on block theory method (BTM).The inputs of this function are joints and tunnel axis orientations and also the code indicating the block and the binary output indicates the removability of the block: where Flag is a binary 0 or 1 that indicates the removability of the block (0: not removable, 1: removable) and BC is the block code.

Key-Block Geometry Function
. By determining the removable blocks with kinematic analysis, it is possible to check the kinetic instability of these blocks (Figure 3).block.In this study, the resultant active and passive force vectors are obtained from ( 12) and ( 13), respectively [15,16]: where  is active force vector,  is block weight vector,  is shotcrete weight vector,  is water force vector,  is seismic force vector,   is passive force vector,  is shotcrete shear resistance force vector, and  is bolt force vector.The probable sliding direction is computed by considering active forces ( vector).By knowing the sliding direction, it is possible to calculate the joints shear and tension strength using one of the joint behavior models such as Mohr-Coulomb model (Figure 3).Finally, a limit equilibrium analysis is performed to calculate the safety factor.The general formulation is as follows: Resisting forces (e.g., shear, tensile, and support forces) Driving forces (e.g., weight, water, and seismic forces) . ( In this study the performance function of each key block has been defined by block factor of safety as follows: The reliability problem of structural stability has been defined by placing of tunnel performance function (15) in general FORM problem (5).To solve resulting problem, a computer code has been developed.Due to high complexity of the tunnel performance function, it is impossible to introduce an analytical and or a closed form solution for the resulted optimization problem.To overcome this issue, a numerical optimization method based on "interior-point algorithm" [17] has been used.The optimization problem was solved using TuPSA code for a typical example which is presented in this work.The results show the great capability of FORM method to reliability analysis of structural instability with low failure probability instead of Monte-Carlo simulation.

Typical Example
Deterministic and probabilistic properties of a circular tunnel and surrounding rock mass are given in Table 1.Clear correlations are observed between joints mechanical properties as well as joints geometrical parameters.Correlation coefficient matrix of random parameters is illustrated in Table 2.The stability of the tunnel has been studied by considering 5 different scenarios from deterministic to fully probabilistic analysis.These scenarios are as follows: (S1) The parameters variability has no influence on the stability and only the mean values are effective in analysis.A deterministic safety factor and block performance were determined using the algorithm presented in Figure 3.
(S2) Only the joints strength parameters (, , Ten) variability has an influence on the stability and these parameters are assumed as the independent variables.
(S3) It is similar to scenario 2, but with considering the effects of correlations between random variables.
(S4) Only the joints orientation parameters ( 1 ,  1 ,  2 ,  2 ,  3 ,  3 ) variability has an influence on the stability and these parameters are assumed as the correlated random variables.
(S5) All the joints properties (mechanical and geometrical properties) variability has an influence on the stability and these parameters are assumed as correlated random variables.
Here, the reliability analysis steps for the general case (scenario 5) are expressed.The other solutions are simply derived from this case.
(1) Determine physical random vector using the standard normal variables.

Results and Discussion
By deterministic stability analysis using the mean values of joints properties (scenario 1), six removable blocks have been found and analyzed using the kinematic and kinetic instability analysis.The position of these blocks around the circular tunnel has been illustrated in Figure 5.The performance function determination algorithm which is illustrated in Figure 3 has been used in this analysis.The results of key blocks removability, geometry, safety factor, and performance are presented in Table 3.The probabilistic analysis using the TuPSA method has been carried out for scenarios 2-5 and the results are presented in Figures 6 and 7. Figure 6 illustrates the reliability index of each key block for different scenarios.Figure 7 illustrates the probability of each block failure for different scenarios.Note that the small probability of failures is illustrated in logarithmic scale.For validation, a Monte-Carlo simulation has been carried out for scenario 2. In the Monte-Carlo simulation, the number of iterations is selected so that, based on (3), the coefficient of variation of   (CoV   ) becomes less than 0.1.A comparison between results of TuPSA method and Monte-Carlo simulation for scenario 2 is presented in Table 4.
A comparison between reliability analysis using correlated and uncorrelated joint friction and cohesion strength (Figures 6-7) indicates that considering the negative correlation between joints friction and cohesion strength parameters reduces the probability of block failure and increases the tunnel reliability index.
As can be seen from Figures 6-7 (scenario 4), for the removable blocks with high safety factor (greater than 3),  the variability of joints orientation has very little impact on the event of failure.By performing a comparison between scenarios 2, 3, and 4, it can be concluded that for these blocks the variability of strength parameters has a more effective role in failure events.By comparing scenarios 2, 3, and 4 for block [1 1 1], it can be concluded that, for the marginal blocks which have a low safety factor (between 1 and 2), both of the geometrical and strength variabilities are important, separately as well as together, in the failure event.
Table 4 illustrates the results of FORM analysis and the Monte-Carlo simulation for the assumed case study.Considering the running time and accuracy of results, it is clear that FORM is more efficient and more precise than ordinary Monte-Carlo simulation (in finite iterations).This fact is more evident when the probability of failure is too small and consequently it requires too many simulation iterations to evaluate accurate probability based on (3).Table 4 illustrates that, on the contrary to FORM analysis, the simulation iterations and simulation time increase when the probability of failure decreases.

Conclusions
A reliability-based design methodology for tunnel structural instability (TuPSA) is developed.Many sources of uncertainties can be modeled and analyzed in evaluating the probabilistic tunnel stability.The physical random variables are considered to have arbitrary probability distribution with any types of correlations.A combination of tunnel structural stability models and first order reliability method (FORM) has been used for reliability analysis of structural instabilities around the circular tunnels.In this method, by transforming the performance functions from the space spanned by physical random variables into the noncorrelated standard normal random variables space, the reliability index and probability of failure are calculated via optimization techniques.
A case study with five scenarios from deterministic to full probabilistic state is presented to show the implementation of the TuPSA model and to obtain an estimate of the target value for each removable key block.It is possible to replace factor of safety with the performance of tunnel considering the tunnel design goals.Comparison between probabilistic and deterministic analyses results indicates that probabilistic results, including the design point and probability of failure, are more rational than deterministic factor of safety.Comparison between the failure probabilities obtained from FORM and those obtained from Monte-Carlo simulations indicates that FORM is the most efficient reliability method and its accuracy is very high and more than sufficient.
Deterministic methods are the base of design and analysis in rock engineering.However, the reliability analysis presented in this paper demonstrates the applicability of the statistical and probabilistic analysis in problems dealing with uncertainties and variability.However reliability analysis could consider the effects of uncertainties in stability analysis, faced with the same problem as the deterministic approaches.Both of the deterministic and probabilistic analyses use the same failure models that may have inherent weaknesses.

Figure 3 :
Figure 3: Procedure for TuPSA performance function definition for a specific block.

Figure 5 :
Figure 5: Removable key blocks positions around the tunnel.

Figure 7 :
Figure 7: Failure probability of each block considering different scenarios.

Table 1 :
At first, the geometry of biggest/critical blocks (key blocks) around the tunnel is determined based on BTM.Key-block geometry function determines the key-block volume,   , block face area,   , excavation face area of the block,   , trace length of joints in excavation rim,   , and angular intervals of block outcrops in the tunnel,  1 and  2 (Figure4).The key-block geometry function could be presented by [  ,   ,   ,   ,  1 ,  2 ] Properties of the circular tunnel.  : dip of th joint set; * *   : dip direction of th joint set.
= key-block-geometry (TO, , JD, BC) .(11)3.3.Factor of Safety Function.By defining a specific key-block geometry, it is possible to calculate the active and passive force vectors (in global ENU coordinate system) acting on the key *

Table 3 :
Results of deterministic analysis considering the mean values of random parameters (scenario 1).
[6]UL: key block created by combination of 1st and 2nd joints upper half spaces and 3rd joint lower half space (based on BTM; see[6]).

Table 4 :
Results of the reliability analysis of tunnel stability, considering the joints strength parameters (C, , Ten) as the independent random variables (scenario 2).* : it is not possible to calculate in reasonable time period with ordinary processing units.