# Safety Evaluation of Underground Caverns Based on Monte Carlo Method

Authors: Qifeng Guo; Zhihong Dong; Meifeng Cai; Fenhua Ren; Jiliang Pan
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In order to study the influence of joint fissures and rock parameters with random characteristics on the safety of underground caverns, several parameters affecting the stability of surrounding rock of underground caverns are selected. According to the Monte Carlo method, random numbers satisfying normal distribution characteristics are established. A three-dimensional model of underground caverns with random characteristics is established by discontinuous analysis software 3DEC and excavation simulations are carried out. The maximum displacement at the numerical monitoring points of arch and floor is the safety evaluation index of the cavern. The probability distribution and cumulative distribution function of the displacement at the top arch and floor are obtained, and the safety of a project is evaluated.


---

## Body

## 1. Introduction

More underground caverns with larger space are being developed in China and other countries. It is crucial to ensure the safety of underground cavern construction during the excavation and service process. Although the underground caverns should be excavated in stable rock mass, it is inevitable that the cavern will pass through weak structural zones of rock. It has been known that the strength of the rock will be greatly reduced by faults and joints in the rock. And the mechanical properties and weak joints are randomly distributed in the rock mass. Therefore, it is important to consider the randomness to evaluate the safety of the underground caverns.Qian and Zhou [1] carried out numerical simulation and field monitoring to investigate rock deformation and failure behaviors during excavation of underground cavern group under high in situ stresses in Jinping I Hydropower Station and found that the damage behaviors of rock masses subjected to high in situ stresses are characterized by the formation of multifractured zones. In consideration of the effect of the fault and joint on the stability of surrounding rock, Liu et al. [2] studied the formation and propagation mechanism of the surrounding rock fracture zone in the discontinuous jointed rock mass with different joint dip angles and densities. Leng et al. [3] developed a ubiquitous joint model in consideration of the nature of rock blocks and joints and studied the influence of the trend and tendency of weak joints on the deformation of surrounding rock. Zhou et al. [4] developed a micromechanics-based model to quantify the effects of various parameters on the complete stress-strain relationship and evaluated the complete stress-strain relationship and strength for the crack-weakened rock at the underground cavern complex of the Ertan Hydroelectric Project. Fracture mechanism of surrounding rock mass around underground caverns of Jinping I Hydropower Station was further investigated by stimulating growth and coalescence of cracks with the different scales by Zhou [5]. Niu et al. [6] analyzed the influence of fault occurrence on the stability of tunnel surrounding rock from two aspects of fault dip and angle. Zhou et al. [7] creates a novel meshless numerical method, called general particle dynamics (GPD), to simulate samples of rock-like brittle heterogeneous material containing four pre-existing flaws under uniaxial compressive loads. Lu et al. [8] studied the influence of the fault on the deformation and plastic zone of tunnel surrounding the rock. By comparing different inclination angles and different joint positions, Zhou et al. [9] evaluated the control effect of joints on surrounding rock stability from the aspects of shear stress distribution, displacement distribution and plastic flora are caused by underground excavation. Zhu et al. [10] invest the effects of the nonoverlapping length and flaw angle on the complete stress-strain curves, the stress of crack initiation, the peak strength, the peak strain, and the elastic modulus. Xi et al. [11–13] investigated the effect of heterogeneity on concrete cracking by cohesive crack models. Huang and Huang [14] studied the influence of natural in situ stress field on the stability of the surrounding rock of underground caverns with faults. Wang and Ren [15] studied the influence of the characteristics of joints on the stability of tunnel-surrounding rock in jointed fractured rock mass. Jia et al. [16] studied the influence of joints under different tunnel strikes on the stability of tunnel-surrounding rock.Although considerable achievements have been made in the study of the influence of fault and joint properties on the stability of surrounding rocks, the physical parameters of rocks and joints in most studies were deterministic and the majority of them were calculated by some specific parameters. However, for rocks and joints, the stability of their engineering systems is usually affected by a large number of uncertain factors. The Monte Carlo method is based on the “random number” calculation method, which is a method of calculation [17]. The basic idea is to understand a system through a large number of random samples and then get the value for calculation. The idea was born in the United States in 1940s “Manhattan project,” named after the gambling city of Monte Carlo, a symbol of probability. When the problem is the probability of a random event, or the expected value of a random variable, some “experimental” method is used to estimate the probability of the random event at the frequency of such an event, or to obtain some numerical characteristics of the random variable and take them as the solution of the problem. The Monte Carlo method theoretically requires a great deal of experimentation. The more the experiments, the more accurate the results. With the development of computer technology, the Monte Carlo method has been rapidly popularized in recent years. Modern Monte Carlo methods can take advantage of the high speed of a computer, which makes a laborious experimental process fast and easy. Liu et al. [18] established a probabilistic water quality prediction model based on Monte Carlo simulation. According to the calculation results, the range and probability of river water quality can be calculated by statistical analysis. Jiang et al. [19], based on the Monte Carlo simulation sampling method, constructed a method for evaluating the consequences of natural gas pipeline leakage and injection fire accident. Sun and Li [20] used the Monte Carlo simulation method to build an income risk assessment model for urban infrastructure projects based on the project internal rate of return and net present value. The physical quantities after Monte Carlo treatment can well reflect the randomness, which is more conducive to truly reflect the surrounding rock displacement during the excavation of underground caverns under the action of stochastic characteristics of rocks.In this paper, three parameters, elastic modulus of rock, joint orientation and axial angle of tunnel, and joint spacing are selected to calculate the maximum vertical displacement of the tunnel. Firstly, the parameters are randomized by the Monte Carlo method. Then, the probability distribution function of the maximum vertical displacement of the tunnel is obtained by quantitative numerical simulation. Finally, the obtained distribution function is used to evaluate the safety of an underground cavern during excavation.


## 2. Three-Dimensional Numerical Model with Stochastic Characteristics

### 2.1. Generating Random Parameters

The stochastic characteristics of rock mass and joint network can be obtained by field measurements. The distribution functions of the objective parameters are then built by field data obtained from actual measurement on the rock mass and joints. The objective parameters (e.g., joint orientation and spacing) in the numerical model must follow the probability distribution model.The physics of rock mass and joint usually satisfies the random variable of nonuniform distribution. For the simulation of the random characteristics of the rock and joint surface, we need to get the random number of corresponding distribution. According to the Monte Carlo method, the random number of uniform distribution in [0,1] interval is generated randomly, and then the random number of corresponding parameters is generated on this basis. For rocks and jointed faults, the values of their parameters are highly random, which can be regarded as random variables obeying normal distribution [21]. According to Box–Muller transform algorithm, the random number satisfying normal distribution is(1)xi=μ+σ−2lnrsin2πt,where σ and μ are the mean and standard deviation of the known normal sample data, respectively, and r and t are independent random numbers which obey uniform distribution in (0,1) range.


### 2.2. Numerical Model with Stochastic Characteristics for Underground Cavern Model

Discontinuous analysis software 3DEC is employed in this paper for the numerical simulations. The numerical model with stochastic characteristics for underground cavern model process is established as follows: firstly, blocks with a certain size are generated in the dimensions of the cavern. And blocks in the space of tunnels are defined as a group which could be removed to simulate the excavation process. Secondly, the joints with different values of tendency and inclination of joints are generated. The distribution function is to generate random numbers that obey uniform distribution in (0,1) and generate random numbers that meet normal distribution, respectively, according to equation (1). Finally, a number of stochastic numerical models following the probability distribution model are built. The random process is achieved by Fish programming in 3DEC.The normal distribution satisfied by the parameters after each randomization is shown in Table1. Figure 1 shows a typical example for the numerical model. The size of the numerical simulation sample is 100 m × 100 m × 100 m. The underground powerhouse is 59 m in height, 24 m in width, and 100 m in axial length. The number of elements of the numerical model is 1264,188, and the average size of elements is 0.8 m. Firstly, the Fish function of response is created by using def command and the Fish function of response is created by using random numbers. Using jset command, the function is invoked at the parameter assignment point. Fifteen joints through the model are generated according to random characteristics.Table 1
Normal distribution of random parameters.
Normal distribution characteristics
Modulus of elasticity (GPa)
Joint spacing (m)
Inclusion angle between joint strike and tunnel (°)
Mean value
13
2
35
Variance
2
0.5
3Figure 1
Underground cavern excavation model with stochastic characteristics.


## 2.1. Generating Random Parameters

The stochastic characteristics of rock mass and joint network can be obtained by field measurements. The distribution functions of the objective parameters are then built by field data obtained from actual measurement on the rock mass and joints. The objective parameters (e.g., joint orientation and spacing) in the numerical model must follow the probability distribution model.The physics of rock mass and joint usually satisfies the random variable of nonuniform distribution. For the simulation of the random characteristics of the rock and joint surface, we need to get the random number of corresponding distribution. According to the Monte Carlo method, the random number of uniform distribution in [0,1] interval is generated randomly, and then the random number of corresponding parameters is generated on this basis. For rocks and jointed faults, the values of their parameters are highly random, which can be regarded as random variables obeying normal distribution [21]. According to Box–Muller transform algorithm, the random number satisfying normal distribution is(1)xi=μ+σ−2lnrsin2πt,where σ and μ are the mean and standard deviation of the known normal sample data, respectively, and r and t are independent random numbers which obey uniform distribution in (0,1) range.


## 2.2. Numerical Model with Stochastic Characteristics for Underground Cavern Model

Discontinuous analysis software 3DEC is employed in this paper for the numerical simulations. The numerical model with stochastic characteristics for underground cavern model process is established as follows: firstly, blocks with a certain size are generated in the dimensions of the cavern. And blocks in the space of tunnels are defined as a group which could be removed to simulate the excavation process. Secondly, the joints with different values of tendency and inclination of joints are generated. The distribution function is to generate random numbers that obey uniform distribution in (0,1) and generate random numbers that meet normal distribution, respectively, according to equation (1). Finally, a number of stochastic numerical models following the probability distribution model are built. The random process is achieved by Fish programming in 3DEC.The normal distribution satisfied by the parameters after each randomization is shown in Table1. Figure 1 shows a typical example for the numerical model. The size of the numerical simulation sample is 100 m × 100 m × 100 m. The underground powerhouse is 59 m in height, 24 m in width, and 100 m in axial length. The number of elements of the numerical model is 1264,188, and the average size of elements is 0.8 m. Firstly, the Fish function of response is created by using def command and the Fish function of response is created by using random numbers. Using jset command, the function is invoked at the parameter assignment point. Fifteen joints through the model are generated according to random characteristics.Table 1
Normal distribution of random parameters.
Normal distribution characteristics
Modulus of elasticity (GPa)
Joint spacing (m)
Inclusion angle between joint strike and tunnel (°)
Mean value
13
2
35
Variance
2
0.5
3Figure 1
Underground cavern excavation model with stochastic characteristics.


## 3. Deformation Analysis Based on the Stochastic Models for Underground Caverns

The mechanical properties of the joints in the numerical model are shown in Table2. The bottom boundary of the whole model is fixed, and the normal displacement of the surrounding boundary is constrained. Gravity is applied to the whole model, and gradually compressive stress boundary is applied to the surrounding boundary. The minimum value of stress is 2 MPa, the maximum value is 10 MPa, and the lateral pressure coefficient of the stress boundary is 1; that is, the surrounding boundary is subject to stress of the same size. Due to the deep location of the underground powerhouse in the actual working condition and the homogeneous properties of surrounding rocks, it is assumed that the model has only one stratum with a density of 2700 kg/m3. According to the above parameters, the excavation of the underground powerhouse is simulated.Table 2
Table of remaining parameter values.
Joint surface cohesion (MPa)
Normal stiffness of joint surface (MPa)
Internal friction angle of joint surface (°)
0.5
0.1
35The maximum displacement of the top arch and the floor is taken as the safety evaluation index of the cavern. The numerical monitoring points are set at the corresponding positions. 100 simulations are carried out, and the vertical displacement contour from 5 typical simulation results is shown in Figures2–5.Figure 2
Vertical displacement of cavern from the typical simulation result 1.Figure 3
Vertical displacement of cavern from the typical simulation result 2.Figure 4
Vertical displacement of cavern from the typical simulation result 3.Figure 5
Vertical displacement of cavern from the typical simulation result 4.From the displacement contours, it can be seen that in the excavation process of underground caverns, the displacement of the top arch is generally larger than the displacement of the bottom plate and the displacement of the side wall of the cavern. The displacement value on both sides of the joint surface has a larger span. Some rock blocks will slide along the joint surface. The denser the joints pass through the cavern, the greater is the vertical displacement of the cavern. Therefore, during the excavation of underground caverns, the crown arch area with dense joint distribution should be paid more attention to.The displacement values at monitoring points in all calculation results are extracted and processed by statistics. The histogram of maximum displacement distribution of the cavern is obtained, as shown in Figures6 and 7.Figure 6
Histogram of top arch displacement distribution.Figure 7
Histogram of floor displacement distribution.We obtained the sampling calculation data of 100 groups in which the data for each group are independent and uniformly distributed to each other. When the sampling number is large enough, the arithmetic mean value of variables of each group converges to its expected value with probability 1. Therefore, the displacements of top arch and bottom plate converge to their expected value, respectively.It can be seen from Figures6 and 7 that the distribution of the displacement values of the top arch and the bottom plate is basically consistent with the normal distribution. After fitting the histogram, we can get the following results:The top arch displacement value satisfies normal distribution:N78,18The displacement of the floor satisfies the normal distribution:N22,10The probability density function of the maximum displacement of the crown arch is as follows:(2)fxtop=16πe−x−782/36.The probability density function of the maximum displacement of the floor is as follows:(3)fxfloor=1210πe−x−222/20.By integrating formulas (2) and (3), the cumulative distribution functions of the displacements of the crown and the floor are obtained, respectively:(4)Fxtop=∫−∞x16πe−x−782/36dx,(5)Fxfloor=∫−∞x1210πe−x−222/20dx.


## 4. Safety Evaluation of Underground Caverns

An underground powerhouse is located in the interior of a mountain with gentle topography on the top and steep bedrock walls on both sides with slopes over 50 degrees. The surface vegetation in the work area is more developed. The bedrock in the factory area is a huge porphyritic quartz two long porphyry in the late Yanshan period of the Songshan unit. In the workshop area, there are mainly three groups of joint fissures with NNE, NE, and NWW directions, mainly steep dip fissures, and a few gentle dip fissures locally.The discontinuous analysis software used in this paper is 3DEC. Through the simulation experiment of underground cavern excavation, the underground excavation space is simplified into workshop cavern, main transformer cavern, and four busbar caverns connecting the above two caverns. The overall length, width, and height of the model are 150 m × 150 m × 150 m. The geometry of the underground is imported as the building drawings. The bottom of the model is fixed, and the side boundary is constrained by the normal phase displacement. On the side of the model, the gradient stress boundary is applied with the maximum value 16 MPa, the minimum value 6 MPA, and the pressure coefficient is 1. The depth of the plant is about 150 m, the elastic modulus of surrounding rock is 12 MPa, and the density is 2700 kg/m3. The axonometric map of the underground powerhouse is shown in Figure 8.Figure 8
Three-dimensional axonometric map of the underground powerhouse model.The top, middle, and lower sections of the model are intercepted to monitor the displacement of the vault and the floor, respectively. The displacement contours of each face are shown in Figures9–11.Figure 9
Vertical Displacement of upstream workshop.Figure 10
Vertical displacement of the middle section of a factory building.Figure 11
Vertical displacement of the downstream powerhouse.The probability of response displacement is calculated according to equations (4) and (5). The results are shown in Table 3.Table 3
Probability of response displacement.
Monitoring value of crown displacement (mm)
Cumulative distribution probability (%)
Floor displacement monitoring value (mm)
Cumulative distribution probability (%)
90.9
96.64
4.9
0.19
153.3
99.99
35.9
70.77
50.3
34.78
36.6
70.79
90.9
96.64
4.9
0.19From the simulation results, it can be seen that the cumulative distribution probability of the maximum displacement of the roof arch in the upper and middle sections of the cavern is maximum when the underground powerhouse is excavated without support. It indicates that during the excavation process, the cumulative distribution probability of the maximum displacement of the roof arch in the upstream and middle sections is very large, the generation of the displacement value is relatively certain. Moreover, since the vertical displacement value of the roof arch in the upstream and middle sections has exceeded the safety range, it should be supported. The cumulative distribution probability of the vertical displacement value in the downstream is small, and the generation of displacement value is relatively random. The probability of the displacement of the chamber bottom plate in the middle and downstream is relatively large, and the displacement value is relatively certain and within the safe range. The probability of the displacement value of the upstream floor is extremely small, so the chamber bottom plate is generally safe in the process of excavation. Because the vertical displacement of the top arch in the upper and middle reaches has exceeded the safe range, it should be supported. The probability of cumulative distribution of vertical displacement value downstream is small, and the safety of the cavern roof arch in the excavation process is high. The displacement of the cavern floor is more likely to occur in the middle and lower reaches, and its value is within the safe range. The displacement of the upstream floor has a very small probability, so the cavern floor is safer in the whole excavation process.


## 5. Conclusion

In this paper, three parameters, elastic modulus of rock, angle between joint direction and tunnel axis, and joint spacing are randomized, and a model of underground cavern excavation with random characteristics is established. The maximum displacement of the roof arch and floor is taken as the safety evaluation criterion for numerical simulation calculation, and the cumulative distribution function of the maximum displacement of cavern is deduced. By evaluating the safety of underground cavern excavation in a project, the following conclusions are drawn:(1)
The parameters are randomized according to normal distribution. After 100 times of numerical simulation experiments, it is found that the maximum displacement of the roof arch and floor of the underground cavern is also approximately normal distribution. According to the calculation results, the cumulative distribution function of the maximum displacement of the cavern is obtained.(2)
The excavation process of an underground power station without support is simulated. The maximum displacement cumulative distribution function is used to evaluate the safety of the maximum displacement value of the cavern. It is concluded that the safety of the roof arch is lower and the floor is higher in the excavation process of the underground powerhouse. It is necessary to carry out the key support at the top arch of the upper and middle sections of the cavern.



---

*Source: 7214720-2020-01-13.xml*