#### Abstract

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 rock-support interaction analysis. The proposed method for estimating the reliability index and determining the tunnel support parameters was introduced in detail. To illustrate the proposed method, it was applied to two tunnels with rock-bolt and combined support systems. The results indicated that the method could obtain accurate solutions with a dramatically improved computing efficiency, guaranteeing the tunnel stability at the same time. The 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.

#### 1. Introduction

There 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–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. Therefore, 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–13]. The 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].

The optimization of engineering structure is essential to structure design in various fields. The 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–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]. The 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–37]. The CCM, which simulates three-dimensional 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 rock-bolts 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]. The 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. The 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.

The objective of this study was to present a practical approach to performing reliability-based designs of rock-support systems. In this study, an RBD was applied to a rock-support interaction design in a tunnel. First, the rock-support interaction analysis was reviewed. Second, a first-order reliability method (FORM) and an RBD were presented. Then, 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.

#### 2. Rock-Support Interaction Designs for Rock Tunnels

The 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 rock-support interactions.

The rock-support interaction analyses, which are based on CCM, require the rock response curves and the support reaction curves, as illustrated in Figure 1. The 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. The surrounding rock was assumed to behave as an elastic-perfectly plastic material. The 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. The 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. The 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 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. The 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.

##### 2.1. Rock Response Curve

The 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. The solutions are summarized in this study’s Appendix section.

##### 2.2. Support Reaction Line

For the support of the reaction line, the support stiffness and maximum support pressure need to be determined. The 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.

The 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; 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.

The maximum support pressure of the support system *p*_{totmax} was determined by the following equations [6]:

##### 2.3. 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). The 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.

#### 3. 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. This 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]. Then, chaotic particle swarm optimization (CPSO) was adapted as an optimization method to search the design variables.

##### 3.1. First-Order Reliability Method (FORM)

The determination of the reliability index is important to a reliability-based design (RBD). In this study, a first-order reliability analysis method (FORM) was used to calculate the reliability index. The 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}. Then, 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). Then, the procedure of determining the reliability index was complemented in Excel VBA with Solver (Figure 4). The code of reliability analysis is listed in Figure 4.

##### 3.2. Chaotic Particle Swarm Optimization

Kennedy and Eberhart originally designed particle swarm optimization (PSO) for large search spaces [46]. The 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;

*V*_{i}= () is its current speed; and

*P*_{i}= (

*p*

_{i1},

*p*

_{i2}, …,

*p*

_{in}) is its best position. Then, 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:where 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 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 can be constrained to the range [, ] in order to control excessive roaming of particles outside the search space. Then, the particles will fly toward a new position according to equation (16). This 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 , *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. Therefore, 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 into the PSO [47]. This well-known logistic equation was incorporated into a simple PSO. The logistic equation was defined as follows [48]:where *μ* is the control parameter; *x* is a variable; and *n* = 0, 1, 2, …. Then, 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 , *r*_{1}, and *r*_{2} of equation (15) were controlled by the following equations:where and *r*_{i}(0) were not 0, 0.25, 0.5, 0.75, and 1.

##### 3.3. Reliability-Based Design Using CPSO

###### 3.3.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. Then, 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. The 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}. This 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.

###### 3.3.2. 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 *i*th failure mode was not less than the target reliability index (), then **d**_{i}, and represented the *i*th element of **d** and the lower and upper permissible values for **d**_{i}, respectively. The 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: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 *i*th failure model, respectively; *d*_{i} represents the design variables; and are the lower and upper limitation of the *i*th 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. The procedure of determining the reliability index using Excel Solver is described in Section 3.1.

###### 3.3.3. Objection Function

An objective function is essential to any optimization problem. Therefore, 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. Then, 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. The 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:

###### 3.3.4. CPSO-Based RBD

The procedure of an RBD includes two loops. The 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. Then, the CPSO was used to search the design variables. The flowchart of this study’s CPSO-based RBD is shown in Figure 5.

##### 3.4. 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. Then, 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. The 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

#### 4. Validations

This 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). The optimal design variables of the supports were determined, and the rock-support interaction of the tunnels was evaluated.

##### 4.1. 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. The in situ stress was *p*_{o} = 10 MPa. The rock mass properties were uncertain and therefore were regarded as random variables with normal distributions. The statistics of the random variables are listed in Table 1. Anchored rock-bolt was used to stabilize the tunnel. The diameter of the rock-bolt (*d*_{b}) was 25 mm. The Young modulus (*E*_{b}) was 200 GPa. The ultimate failure load (*T*_{bf}) and load-deformation constant (*Q*) of the rock-bolt 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.

The spreadsheet of the reliability-based design is shown in Figure 6. The parameters of the CPSO and the searching range of the design variables are shown in Figure 6. The 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). The value of the constraint of reliability was negative (Figure 6), which indicated that the results had met the reliability criterion. The 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. The 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 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.

##### 4.2. Example 2: Combined Support System

A 10.7 m diameter highway tunnel was excavated in frail-quality gneiss, at a depth of 122 m below the surface. The in situ stress was *p*_{o} = 3.31 MPa. The Poisson ratio and the unit weight of the rock were 0.2 and 0.02 MN/m^{3}, respectively. The 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. The 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. The diameter of the rock-bolt (*d*_{b}) was 25 mm. The Young modulus (*E*_{b}) was 207 GPa. The ultimate failure load (*T*_{bf}) and load-deformation constant (*Q*) of the rock-bolt were 0.285 MN and 0.143 m/MN, respectively. The 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.

The spreadsheet of the reliability-based design is shown in Figure 10. The 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). The value of the constraint of reliability was negative (Figure 10). The maximum cycle of the CPSO was 50. The time of computing was approximately 1,725 seconds. Figure 11 shows the curve of the rock-support interaction. The 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.

#### 5. Conclusions

In this study, a novel reliability-based design approach that took the uncertainties into account was proposed to design support systems in tunnels. This new design was applied to two circular tunnels, one with a rock-bolt support system and one with a combined support system. The results showed that the proposed method was able to obtain accurate solutions with a dramatically improved computing efficiency. The proposed method could also be generalized and used for any shape tunnel with a numerical solution. The 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. The 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.

#### Appendix

#### Solutions for the Circular Tunnel Using M-C and Empirical Yield Criterion

In the solution based on Mohr–Coulomb (M-C) yield criterion, it was assumed that the circular tunnel was excavated in a continuous, homogeneous, isotropic, initially elastic rock mass, which had been subjected to hydrostatic far-field stress *p*_{0} and uniform support pressure *p*_{i}. Therefore, if *p*_{i} was less than the critical pressure *p*_{cr}, a plastic zone existed. According to the Mohr–Coulomb criterion, the plastic zone radius *r*_{p} and the inward displacement of the tunnel wall *u*_{ip} could be obtained by the following equations [43]:where *E* is the elastic modulus and is the Poisson ratio. The values of *p*_{cr}, *k*, and *s* were obtained by the following equations:where *φ* denotes the friction angle and *c* is the cohesion.

For the solution based on the empirical yield criterion, if the plastic failure occurred around the tunnel, then the radius of the plastic zone could be obtained by the following equations:

The deformation around the tunnel was obtained by the following equation:where *σ*_{c} represents the uniaxial compressive strength of the intact rock; *m* and *s* are the material constants for original rock mass; and *m*_{r} and *s*_{r} are the material constants for the broken rock mass.

If the deformation around the tunnel was elastic, then the deformation around the tunnel could be obtained by the following equation:

#### Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by the Open Research Fund of the State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, China (Grant no. Z020006).