The Reliability-Based Design Optimization of considering Rock-Support Interaction for Rock Tunnels

Uncertainty is critical to the tunnel design. In this study, a novel reliability-based design (RBD) method was developed by integrating the rock-support interaction and its uncertainties for rock tunnels. In this method, the rock-support interactions were analyzed based on a convergence-confinement method. Uncertainties were estimated by reliability analysis using Excel Solver. Chaotic particle swarm optimization (CPSO) was adopted to search the optimal tunnel design parameters based on the rocksupport interaction analysis. )e proposed method for estimating the reliability index and determining the tunnel support parameters was introduced in detail. To illustrate the proposedmethod, it was applied to two tunnels with rock-bolt and combined support systems. )e results indicated that the method could obtain accurate solutions with a dramatically improved computing efficiency, guaranteeing the tunnel stability at the same time. )e developed method provides an excellent way to deal with the uncertainty in tunnel design. It was proved to be an efficient and effective method for the support designs of rock tunnels.


Introduction
ere are lots of uncertainty problems existing in various engineering systems. To improve the performance of engineering systems, reliability-based design (RBD) has been widely used in various engineering fields [1][2][3][4][5]. RBD involves the finding of balanced designs that are economical and reliable in the presence of uncertainties. Rock-support interaction analyses and designs are of critical importance to the success of any tunneling project [6]. Rock-support interaction analyses are usually affected by some uncertainties which cannot be considered using classic deterministic analysis approaches such as convergence-confinement methods (CCM). However, the overwhelming majority of the previous rock-support interaction studies were based on deterministic models which did not consider the uncertainties of the rock mechanical properties or loading of surrounding rock mass [7]. It has become evident that the precise knowledge of the rock-support interactions is impossible, even if an impractical amount of in situ measurements and laboratory tests are performed. erefore, to consider the uncertainties such as rock mechanical parameters, property of support, and in situ stress, probabilistic approaches have been developed for the designs and analyses of tunnels [8][9][10][11][12][13]. e stability of the tunnel face in heterogeneous soil was investigated using a probabilistic approach in the southern extension of Line 3 of the Tehran subway [14]. A probabilistic approach was proposed for face stability analysis of subway tunnels by consolidating the multiple failure mechanisms into the response surface method [15]. Zhou et al. proposed a new method to analyze the stability of a tunnel using P-wave seismic velocity [16]. Cheng et al. compared the random variables method and the random fields method in reliability analysis of tunnel face [17]. Lv et al. analyzed the reliability of rock-support interactions using a response surface method [18]. Boon et al. developed a concept of interaction diagram to make preliminary design decisions of joint rock tunnel given a target maximum allowable convergence based on discrete element method [19]. e optimization of engineering structure is essential to structure design in various fields. e support design for rock tunnels involves a complex problem and a tradeoff between safety and economy because of various uncertainties present in rock mass properties [20]. RBD can consider the uncertainties of structures through the addition of probabilistic constraints (reliability index) [21]. RBD can overcome some limitations and ambiguities in the load and resistance factor design (LRFD) approach and automatically reveal critical design scenarios and obtain design values of loads and resistance [22]. Recently, RBDs have attracted increasing amount of attention in various fields [23][24][25][26]. Phoon et al. proposed a simplified RBD approach based on the load and resistance factor design and multiple resistance factor design for practical implementation [27]. Chalermyanont and Benson developed a set of RBD charts for mechanically stabilized earth walls [28]. Oreste considered the probabilistic distributions of rock properties and support materials to obtain the probabilistic distribution of the support safety factors and then designed supports [29]. Chalermyanont and Benson developed a two-phase approach based on the RBD method for mechanically stabilized earth walls [30]. Wang et al. applied RBDs to spread foundations and drilled shafts by combing Monte Carlo simulations or expanding the RBDs [31]. To improve efficiency, an updated RBD was proposed, which was then combined with the improved reliability methods [32]. e RBD was performed to compute probabilistic face pressures under different COVs at a target safety level based on the polynomial chaos expansion method [33]. Ji et al. developed the RBD method for geotechnical engineering based on an inverse FORM algorithm [34]. RBD, which can consider the uncertainties in design, has some evident advantages over deterministic optimization designs, and its results are rational and meet the engineering reality. However, optimization and reliability assessments require a great deal of repeated computation, which reduces computational efficiency and limits its implementation in practical rock engineering.
Many researchers have developed various methods to estimate the rock-support interaction for determining the support system [6,[35][36][37]. e CCM, which simulates threedimensional problems as the rock-support interaction based on a two-dimensional simplified approach, is well known and widely used to estimate the tunnel's required load capacity. Gschwandtner and Galler developed a CCM approach considering the time-dependent material of the support by investigating different support scenarios of rockbolts and shotcrete [38]. Paraskevopoulou and Diederrichs estimated the time-dependent deformation of the tunnel using the CCM [39]. Fuente et al. applied the CCM to the full-face excavation of circular tunnels with a stiff support system [40]. Su et al. evaluated the stability of the tunnel in the weak rock using the CCM [41]. e CCM provides an excellent way to determine the support system for the circular tunnel in the preprimary design phase. However, the traditional CCM method cannot deal with the uncertainty that influences the analysis of rock-support interaction. In this study, a reliability-based design method was developed to overcome the above issues by combining the CCM, FORM, and RBD. e optimal approach is essential to the reliability analysis and RBD. Various bioinspired soft computing approaches have been developed for engineering systems. In this study, the Excel Solver and chaotic PSO were adopted in reliability analysis and RBO to avoid the local minimum of optimal problem. e objective of this study was to present a practical approach to performing reliability-based designs of rocksupport systems. In this study, an RBD was applied to a rocksupport interaction design in a tunnel. First, the rocksupport interaction analysis was reviewed. Second, a firstorder reliability method (FORM) and an RBD were presented. en, the procedure reliability-based design of the tunnel was presented in detail. Finally, the proposed method was applied to two tunnels related to the rock-bolt system, the combined support system of rock-bolt and shotcrete, and some conclusions were drawn.

Rock-Support Interaction Designs for
Rock Tunnels e potential for instability in the surrounding rock in tunnels is an ever-present threat to both the safety of the workers and equipment in the tunnel. Support structures (rock-bolt, shotcrete, steel set, and so on) are significant to maintaining the stability of the surrounding rock in tunnels. Specific support designs for tunnels may be required in order to reduce the potential instability to an absolute minimum [42]. Rock-support interaction analyses are critical parts of the support designing of tunnels. Convergence-confinement methods (CCM) provide a procedure to evaluate rocksupport interactions. e rock-support interaction analyses, which are based on CCM, require the rock response curves and the support reaction curves, as illustrated in Figure 1. e two curves in a circular tunnel with hydrostatic stress are presented in the following paragraphs.
In this study, in order to present the concepts of the rock-support interactions, a circular tunnel was analyzed in a hydrostatic stress field, in which the horizontal and vertical stresses were equal. e surrounding rock was assumed to behave as an elastic-perfectly plastic material. e failure was assumed to occur with zero plastic volume change. Since the support was modeled as an equivalent internal pressure, the reinforcement provided by the support could not be taken into account in this type of simple model. A typical plot of the rock response curve is presented in Figure 1. e support pressure p i in Figure 1 is that which was provided by the rock mass in the tunnel face through which the tunnel was being advanced. At a distance of approximately one diameter ahead of the tunnel, the rock mass was not influenced by the presence of the tunnel, and the support pressure p i equaled the in situ stress p o , which corresponded to point A on the rock response curve. As the tunnel advanced, the support provided by the rock mass diminished, and the rock mass responded elastically up to point B, at which point a plastic failure of the rock mass was initiated. e radius r p of the plastic zone and the radial convergence u both increased as the support pressure decreased, as illustrated in Figure 1. Eventually, at approximately two tunnel diameters behind the face, the support pressure p i provided by the face had 2 Advances in Civil Engineering decreased to zero, and the radial convergence u reached its final value (point D). Once the support had been installed, and in full and effective contact with the rock, the support started to deform elastically, as shown in Figure 1. e displacement u o by which the support is installed is shown in Figure 1 before the support became effective. Depending upon the characteristics of a support system, it will deform elastically in response to the closure of the tunnel as the face advances away from the point under consideration. Equilibrium will be achieved (point C) if the support reaction curve intersects with the rock response curve.

Rock Response Curve.
e rock response curve is defined as the relationship between the fictitious internal pressure p i and the radial displacement of the wall u r , as shown in Figure 1. For a circular tunnel, the relationship can be constructed from the elastoplastic analysis of the stresses and strains developing in the rock mass with an analytical solution. In this study, a Mohr-Coulomb (M-C) yield criterion [43] and an empirical yield criterion [6] were adopted to illustrate the application of CCM in rock-support interaction analyses. e solutions are summarized in this study's Appendix section.

Support Reaction Line.
For the support of the reaction line, the support stiffness and maximum support pressure need to be determined. e support stiffness and maximum support pressure of anchored rock-bolt can be determined according to the following equations [6]: where L is the length of the rock-bolt; d b is the diameter of the rock-bolt; E b is the elastic modulus of the rock-bolt; Q denotes the load-deformation constant for the rock-bolt; T bf represents the ultimate failure load from pull-out test; and S c and S L are the circumferential and longitudinal spacing of the rock-bolt, respectively. Figure 2 shows the spreadsheet for computing the support stiffness and maximum support pressure of the rock-bolt. e support stiffness and maximum support pressure of the shotcrete are given by the following equations [6]: where E c is the elastic modulus of the concrete or shotcrete; v c is the Poisson ratio of the shotcrete; t c denotes the thickness of the shotcrete; and σ c represents the uniaxial compressive strength of the shotcrete. For the combined support system, the stiffness of the support system k total was determined by the following equation: where k 1 and k 2 are the stiffness of the support systems 1 and 2, respectively. e maximum support pressure of the support system p totmax was determined by the following equations [6]:

Rock-Support Interaction Analysis.
Once the rock response curve and support reaction curve were determined, a rock-support interaction analysis was completed to determine the equilibrium point (point C in Figure 1). e code for computing the equilibrium point is shown in Figure 3. Figure 3 shows the spreadsheet of the rock-support interaction analysis based on an empirical yield criterion.

Reliability-Based Design for a Rock Tunnel
For reliability-based design, the optimization and determination of the reliability index require a great deal of repeated computation. is computation time decreases the efficiency and limits the application of RBDs in practice. In this study, Low and Tang introduced the algorithm to compute the reliability index using Excel Solver software [44]. en, chaotic particle swarm optimization (CPSO) was adapted as an optimization method to search the design variables.

First-Order Reliability Method (FORM).
e determination of the reliability index is important to a reliabilitybased design (RBD). In this study, a first-order reliability analysis method (FORM) was used to calculate the reliability index. e Hasofer-Lind index has been previously widely used [45], for which the matrix formulation for a correlated normal is as follows: where X is a vector representing the set of random variables x i ; μ is the vector of the mean values; C is the covariance matrix; and F is the failure domain. Based on the Hasofer-Lind index, Low and Tang presented a practical and transparent FORM procedure, in which the constrained optimization was based on the perspective of an expanding ellipsoid in the original space of the basic random variables [45]. Next, in order to obviate the computation of the equivalent normal means and equivalent normal standard deviations, they proposed a new efficient algorithm for FORM by varying the dimensionless number n i . en, equation (10) could be rewritten as follows: where R is the correlation matrix and n is a column vector of n. When the value of n i varied during the strained optimization, the corresponding value of x i was calculated as follows: In this study, equation (11) was used to determine the reliability index in Excel Solver software. Excel Solver was used to solve the optimization in equation (11). en, the procedure of determining the reliability index was complemented in Excel VBA with Solver ( Figure 4). e code of reliability analysis is listed in Figure 4.

Chaotic Particle Swarm Optimization.
Kennedy and Eberhart originally designed particle swarm optimization (PSO) for large search spaces [46]. e technique involves simulating social behavior among individuals (particles) which are "flying" through a multidimensional search space.
In regard to the minimum problem, by supposing that f(X) is the objective function, then X i � (x i1 , x i2 , . . ., x in ) will be the current position of the particle; en, the best position of particle i can be computed according to the following equation: If the population is s, and is the global best position among all the particles, then According to the PSO theory, the following equations represent the process of the evolution:

Advances in Civil Engineering
where v i is the velocity of particle i for the distance to be traveled from its current position; x ij is its position; p ij is its best previous position; and p g is the best position among all the particles in the population; r 1 and r 2 represent two independently uniformly distributed random variables with range [0, 1]; c 1 and c 2 denote the positive constant parameters called the acceleration coefficients and control the maximum step size, respectively; and w is the inertia weight that controls the impact of the previous velocity of a particle on its current velocity.
In a standard PSO, equation (15) was used to calculate the new velocity according to its previous velocity, as well as to the distance of its current position from both its own best historical position and its neighbors' best position. Generally speaking, the value of each component in v i can be constrained to the range [, v max ] in order to control excessive roaming of particles outside the search space. en, the particles will fly toward a new position according to equation (16). is process was repeated until a user-defined stopping criterion was reached.
In PSO, the parameters (for example, the inertia weight factor) are crucial for the efficient identification of the optimum solution. Many researchers believe that the parameters w, r 1 , and r 2 in equation (15) are the key factors affecting the PSO convergence. In a simple PSO, the inertia weight cannot ensure that the optimization is entirely ergodic within the phase space due to its randomness. Moreover, parameters r 1 and r 2 cannot ensure the entire ergodicity within the search space. erefore, to enrich the search behavior and avoid entrapment at a local optimum, Zhao et al. proposed a novel approach that combined chaotic mapping with certainty, ergodicity, and stochastic property   (2) For i = 1 To N P(i) = Pi(i) u(i) = ui(i) Next y = 100 i = 1 Do While (y > 0) x0 = u(i) y0 = P(i) x1 = u(i + 1) y1 = P(i + 1) y = (Psm * (x0 -uso) / usm -y0) * (Psm * (x1 -uso) / usm -y1) Advances in Civil Engineering into the PSO [47]. is well-known logistic equation was incorporated into a simple PSO. e logistic equation was defined as follows [48]: where μ is the control parameter; x is a variable; and n � 0, 1, 2, . . .. en, although the logistic equation was deterministic, it exhibited chaotic dynamics when μ � 4, and x 0 was not 0, 0.25, 0.5, 0.75, and 1. In this study, the parameters w, r 1 , and r 2 of equation (15) were controlled by the following equations: where w(0) and r i (0) were not 0, 0.25, 0.5, 0.75, and 1.

Performance Function of the Reliability Analysis.
In this study, in order to investigate the reliability analysis of a tunnel and its support system, two failure modes of a tunnel (for example, the deformation exceedance criterion and the support system failure) were considered. en, the performance function was obtained by the displacement solutions of the tunnel wall, and the support pressure of the support system is as follows: where u r is the inward displacement of the tunnel wall; u max denotes the maximum deformation of the tunnel; p i represents the support pressure of the support system; and p max is the maximum capacity of the support system. e performance function became negative when the inward displacement of the tunnel wall was u r ≥ u max , or the support pressure was p i ≥ p max . is indicated that the tunnel or support system would experience failure. u max and p max were calculated by using equations (1)-(9) in Section 2. u r and p i could also be calculated using the equation in this study's Appendix section.

Reliability-Based Design.
In this study, d denoted the design variables in the engineering system, and C(d) denoted the objective function of the design. By assuming the design requirement was that the reliability index for the ith failure mode was not less than the target reliability index (), then d i , and represented the ith element of d and the lower and upper permissible values for d i , respectively. e purpose of the RBD was to determine a set of design variables that could minimize the cost function (C(d)) without violating any of the design requirements. Mathematically, this problem could be expressed as follows:  Advances in Civil Engineering where C(d) is the objective function, such as the cost; β i (d) and β i T are the reliability index and reliability constraint of the ith failure model, respectively; d i represents the design variables; and are the lower and upper limitation of the ith design variables, respectively; and n m and n d denote the number of reliability constraints and design variables, respectively. An RBD differs from the determination optimization designs and involves reliability constraints. In this study, FORM (equation (8)) was used to compute the reliability index. e procedure of determining the reliability index using Excel Solver is described in Section 3.1.

Objection
Function. An objective function is essential to any optimization problem. erefore, in order to address the problems of a reliability-based design, the design variables need to search for the optimal value based on the objective function using CPSO. en, the problems can be converted into the following nonconstrained optimization form using a penalty method: where f i (d) is the penalty function and M is the penalty coefficient. e choice of the penalty coefficient M is crucial for improving the convergence and accuracy of equation (23). In this study, f i (d) was defined as follows:

CPSO-Based RBD.
e procedure of an RBD includes two loops. e inner loop is used to determine the reliability index based on FORM, and the outer loop is used to search the design variables. In this study, Excel Solver software was adopted to calculate the reliability index. en, the CPSO was used to search the design variables. e flowchart of this study's CPSO-based RBD is shown in Figure 5.

Rock-Support Interaction Analysis Using a Reliability-Based Design.
In this study, Excel Solver software was used to compute the reliability index using a FORM method in the inner loop.
en, CPSO was adopted to search the optimal design variables in the outer loop. A rock-support interaction model was used to determine the equilibrium point, and the support pressure and deformation of the rock mass were calculated to determine the performance function in the reliability analysis. e stepwise procedure was as follows: Step 1. Determine the design variables and their ranges according to the tunnel problem Step 2. Determine the random variables and statistical parameters Step 3. Determine the parameters of the CPSO Step 4. Activate Excel Solver to compute the reliability index using a FORM method Step 5. Compute the objective function value and search the optimal design variables using CPSO

Validations
is study illustrated and analyzed the performance of the approach mentioned above, with one tunnel characterized by a rock-bolt support system and one tunnel containing a combined support system (rock-bolt and shotcrete). e optimal design variables of the supports were determined, and the rock-support interaction of the tunnels was evaluated.

Example 1: Tunnel with a Rock-Bolt Support System.
A six-meter diameter tunnel (r � 3 m), which had been excavated in a fair-quality, blocky sandstone rock mass, was designed with a rock-bolt support system. e in situ stress was p o � 10 MPa. e rock mass properties were uncertain and therefore were regarded as random variables with normal distributions. e statistics of the random variables are listed in Table 1. Anchored rock-bolt was used to stabilize the tunnel. e diameter of the rock-bolt (d b ) was 25 mm. e Young modulus (E b ) was 200 GPa. e ultimate failure load (T bf ) and load-deformation constant (Q) of the rockbolt were 0.254 MN and 0.143 m/MN, respectively. In this study, the initial deformation before the rock-bolt was installed; the length of rock-bolt and the circumferential and longitudinal space of the rock-bolt were successfully determined.
e spreadsheet of the reliability-based design is shown in Figure 6. e parameters of the CPSO and the searching range of the design variables are shown in Figure 6. e initial deformation before the rock-bolt was installed, length of the rock-bolt, and circumferential and longitudinal spaces of the rock-bolt were 5 mm, 2 m, 0.5 m, and 1.49 m, respectively ( Figure 6). e value of the constraint of reliability was negative (Figure 6), which indicated that the results had met the reliability criterion. e time of computing was approximately 1,724 seconds (less than half an hour), which was found to be acceptable for the rock tunnel design. Figure 7 shows the rock-support interaction curve for the tunnel based on the determined support parameters. As shown in Figure 7, the deformation of the tunnel met the demands of the tunnel, and the tension force was equal to the maximum support pressure. In order to improve the safety of the tunnel, a performance function of a support system was introduced (equation (21)). Figures 8 and 9 detail the converging process of the proposed method. Figure 8 shows the variation of the objective function with the cycle. e results indicated that the proposed method displayed a better computing efficiency and the ability to determine the global solution. It located the solution at approximately cycle 50, as shown in Figure 8. In Advances in Civil Engineering other words, it only required half of the 1,724 seconds (less than 15 minutes) to determine the support parameters. It demonstrated a remarkably improved computing efficiency for the RBD. Figure 9 shows the evolutionary process of each individual, which were the same results as those detailed in Figure 8. It was observed that the value of the objective function almost coincided with cycle 50 and cycle 100 in this study. e Young modulus (E), material constants of the original rock mass (m, s), and material constants of the broken rock mass (m r , s r ) were found to be normal distributions of the random variables. e statistical parameters of the random variables are listed in Table 2. A combined support system with anchored rock-bolt and shotcrete was used to stabilize the tunnel. e diameter of the rock-bolt (d b ) was 25 mm. e Young modulus (E b ) was 207 GPa. e ultimate failure load (T bf ) and load-deformation constant (Q) of the rock-bolt were 0.285 MN and 0.143 m/MN, respectively. e elastic modulus (E c ), uniaxial compressive strength (σ c ), and Poisson's ratio (]) of the shotcrete were 20.7 GPa, 34.5 MPa, and 0.25, respectively. In this study, the initial deformation before the rock-bolt was installed, as well as the thickness of the shotcrete, length of the rock-bolt, and the circumferential and longitudinal spaces of the rock-bolt were determined. e spreadsheet of the reliability-based design is shown in Figure 10. e initial deformation before the rock-bolt was installed, along with the length, circumferential space, longitudinal space, and thickness of the lining of the rock-bolt which were 5 mm, 2 m, 0.5 m, 1.49 m, and 50 mm, respectively ( Figure 10). e value of the constraint of reliability was negative ( Figure 10). e maximum cycle of the CPSO was 50. e time of computing was approximately 1,725 seconds. Figure 11 shows the curve of the rock-support interaction. e deformations of the rock mass, rock-bolt, and lining were determined to be safe and met the requirements of the tunnel. Figure 12 shows the variation process of the design variables.

Conclusions
In this study, a novel reliability-based design approach that took the uncertainties into account was proposed to design support systems in tunnels. is new design was applied to two circular tunnels, one with a rock-bolt support system and one with a combined support system. e results showed that the proposed method was able to obtain accurate solutions with a dramatically improved computing efficiency. e proposed method could also be generalized and used for any shape tunnel with a numerical solution. e reliability index was calculated using various reliability methods, with the exception of FORM. A CPSO method was utilized in the demonstration cases in this study. However, any search optimization method could be applied, such as a genetic algorithm, gradient-based methods, and so on. e reliability-based design of the rock-support interaction design was found to conform to the practical requirements of tunnels. It provided a rational and reliable way to conduct the analyses of tunnel stability and support designs in tunnel projects. Rock response curve Support system Figure 11: Rock-support interaction curve for the tunnel using the optimal combined support system. D �