Abstract

In this study, a combined anisotropic random field and GPU-accelerated Cholesky decomposition algorithm is proposed for diversion tunnel excavation. With the MATLAB programming control be used to combine the advantages of GPU and CPU computing, the efficient method of generating large numerical model random fields is realized. Based on the geological structure characteristics of red-layered soft rocks in central Yunnan, the anisotropic rock random field and tunnel excavation with different rotation angles are simulated. The results show that the variability of rotational anisotropy random fields and mechanical parameters of the red-bedded soft rocks in central Yunnan has a large influence on the deformation and stress state of the surrounding rocks at the excavation face of the tunnel. This study can provide theoretical and technical support for the design and construction of relevant rock engineering in the red-layered soft rock area in central Yunnan.

1. Introduction

With the rapid improvement of engineering construction level in China, the development and utilization of underground space are becoming more and more extensive, and the related engineering problems are gradually increasing [1, 2]. The surrounding rock body of a tunnel or underground chamber is a geological material with typical nonhomogeneous characteristics. However, in the traditional geotechnical engineering analysis, the rock is mostly assumed to be a continuous medium with homogeneous properties. This simplified treatment cannot consider the randomness and spatial variability of the actual engineering rock mechanical parameters [3, 4]. In the engineering design stage, due to the lack of profound understanding of the spatial variability of rock geological materials, engineers often cannot anticipate the potential hazards in underground engineering well. Therefore, how to analyze and evaluate the randomness of rock parameters reasonably is of great significance to the stability analysis of underground cavern engineering [47].

At this stage, there are three main methods for the stochastic analysis of mechanical parameters of geological materials. One is to consider the interval distribution of parameters in the calculation process, such as interval finite element [8, 9] and fuzzy finite element methods [10, 11]. These methods can input interval parameters directly and obtain a simulation result in a certain interval range by calculation. However, these methods are difficult to deal with the intrinsic model of nonlinear materials [12]. The second method is the mesoscale modeling of the heterogeneous media based on the statistical property of the microstructure [1319]. However, the analysis scale of this solution is not big enough for the large-scale geotechnic engineering.

The other approach is to characterize the spatial variability of the geological structural properties of the cavern envelope using random field theory [20]. Beacher and Ingra chose the Taylor expansion of the stochastic finite element method to analyze the reliability of ground settlement [21]. Vanmarcke and Grigoriu made pioneering work in the field of soil property parameters [22]. Suchomel and Mašín considered the spatial variable cohesion () and friction angle () of soil mechanical parameters to compare probabilistic methods for assessing slope stability and discussed the influence of parameters’ spatial variability on slope stability [23]. Li et al. proposed a new stochastic finite element method for analyzing slope reliability based on the consideration of spatial variability of slope soil coefficients, which improved the computational efficiency of slope reliability analysis [24]. The effect of spatial variability of infiltration parameters on slope stability was analyzed by Griffiths and Lane [25]. Cho [26] used random field theory to assess slope stability. Jiang et al. [27] conducted random field simulation of geological materials based on the Cholesky decomposition and investigated the effect of different autocorrelation functions on the reliability of slopes.

Based on the literature above, it is evident that the random field method is capable of describing the spatial variability of geotechnical material parameters more accurately. The random field method can be easily implemented in commercial or open-source programs. At present, the random field method for geotechnical materials is mainly applied to the stability analysis of slopes [28]. Less research has been done in the stability analysis of excavation in underground caverns or tunnels.

In this study, the modeling theory of anisotropic random field is optimized and an efficient generation algorithm for random field of underground caverns is developed. Due to a large number of model meshes, the covariance matrix decomposition method may exist the problem of slow computation. Based on the Cholesky decomposition theory and GPU acceleration technology, an efficient random field generation technique with GPU acceleration is proposed. An automated control calculation program was written using MATLAB. With the consideration of the geological structure characteristics of red-layered soft rocks in central Yunnan, an anisotropic rock random field model with different rotation angles was established to analyze the deformation and stress distribution state of the tunnel after excavation in central Yunnan. The influence of the randomness of the surrounding rock parameters on the secondary stress state of the surrounding rock was studied.

2. Methodology

2.1. Isotropic Random Field

There are methods to generate random fields. Currently, two main types of random field generation methods are commonly used [29]. In this study, a generator of spatially correlated random variables is combined with a discretization method to compute, for each node in the random field grid, a random variable that is correlated with other nodes in the random field. These spatially correlated random variables are assigned to a cell or integration point in the finite element model. The random field is generated by a combination of the covariance matrix decomposition method and the midpoint method. Using the covariance matrix decomposition method, a set of correlated random variables can be generated as follows: where is a vector that contains spatially correlated random variables, is the independent zero-mean vector representing the unit variance, normally distributed random variables, and is the Cholesky decomposition of the correlation matrix.

Therefore, the calculation of the correlation matrix is the key to this method. After loading the grid file to obtain the node coordinates, the correlation coefficient matrix for a set of random variables collected in the vector is defined as where is the correlation coefficient of and .

The correlation matrix is symmetric and positive definite, and the correlation coefficient matrix is decomposed into upper triangular matrix and lower triangular matrix by the Cholesky decomposition:

Currently, there are two main types of correlation functions, namely, exponential correlation functions (Exp) and squared exponential (SExp) correlation functions, as follows: where is the threshold value of the correlation function and and are the correlation lengths.

The correlation coefficient matrix is a symmetric matrix that can be obtained by a two-dimensional correlation function [30]. After obtaining the correlation matrix, it can be decomposed using the Cholesky decomposition method. When the SExp correlation function is chosen and the number of nodes is larger than the correlation length, the correlation matrix is nonpositive; i.e., all the eigenvalues of the matrix are no longer positive. In this case, it needs to be corrected. In the modified Cholesky decomposition algorithm, the parameters less than the tolerance value are set to zero to obtain stable results. This algorithm gives an “approximate” Cholesky decomposition. After the decomposition of the correlation matrix, the eigenvalues and eigenvectors are obtained. The decomposition of the correlation matrix is obtained according to where is the matrix containing the eigenvectors and is the diagonal matrix containing the eigenvalues.

Based on the derivation above, the calculation of random field CRV (correlated random variables) can be finally established as shown in where is the standard normal distribution mean, is the standard deviation of the standard normal distribution, and is the standard normal distribution random quantity.

2.2. Anisotropic Random Field

In order to study the modeling theory of anisotropic random fields in geotechnical bodies, it is necessary to recall the definition of correlation length. According to Equations (4) and (5), the value of correlation length is equal to the integration of the correlation function over the full domain. The phase correlation length is the key to measure the spatial variability of geological materials. The difference of its length in different orientations determines the anisotropic characteristics of the random field. The structural characteristics of the random field can be characterized by the correlation length contour. The points in space with the same coefficients associated with the center of the circle are connected to form the correlation length contour, as shown in Figure 1.

and denote the correlation lengths along the main directions of the two coordinate axes. is the correlation length in a direction in coordinate space. is the rotation angle of the coordinate axes with respect to the horizontal datum. It follows that when the correlation lengths along the principal directions of the two coordinate axes are different, the contour changes from an isotropic circle to an ellipse. This describes the random field with horizontal anisotropy characteristics. When the axes are rotated by an angle based on the horizontal line, the length of the projection on the axes also changes. This situation characterizes the random field with rotational anisotropy.

If the anisotropy of rocks is going to be characterized by random fields, suitable modifications according to the definition of random field generation have to be made. In combination with the random field theory described above, this study proposes a method for generating random fields applicable to the anisotropy of rock materials. Along a certain principal axis direction where anisotropic characteristics need to be reflected, the correlation length of that direction is appropriately increased. For example, for the horizontal anisotropy shown in Figure 1(b), the desired anisotropic random field can be obtained by increasing the correlation length of the axial direction. In contrast, the rotational anisotropy shown in Figure 1(c) is not along the principal axis direction. Based on Equations (8) and (9), the new modification and of the correlation length projection along the -direction and -direction can be obtained by changing the angle of the reference coordinate, instead of taking the standard coordinate axis as reference. Taking the exponential correlation function as an example, the random field of rotational anisotropy can be obtained by rotating counterclockwise.

2.3. Anisotropic Random Field Simulation Based on GPU Acceleration

Hardware acceleration refers to that in the process of scientific computing; due to the excessive amount of data to be computed, the computer assigns computational tasks exclusively to a specific piece of hardware to process, reducing the load on the core processor and thus improving computational efficiency. Hardware acceleration technology, widely used in image processing, is commonly known as the graphics processing unit (GPU). Before the advent of the GPU, the CPU had been responsible for a large amount of data processing tasks in the computer. In terms of design thinking, CPUs are suited to complete a single task as fast as possible, usually a line of computation from start to finish. Multimedia data processing is a task that computers are often faced with and usually requires computer hardware that can provide great computational density and efficient parallel computing. In fact, CPUs have difficulty in meeting their huge multimedia data processing requirements in their hardware itself. Software optimization alone could not deal with this problem, until the advent of GPUs.

GPUs are tasked with composing and displaying images of millions of pixels on the screen, requiring parallel processing of millions of tasks. CPUs and GPUs differ greatly in architecture, as illustrated in Figure 2. It can be seen that the GPU has more logical algorithm units (ALUs) and is also multithreaded and parallel, without excessive control and storage functions. Therefore, GPUs have super high computational efficiency and good processing power for dense and large amount of data and complex algorithms. However, the CPU needs to control the whole computation process and store the data at the same time, and it has a small number of logical operation units and can only process one data line in one time period, which is inefficient for large and complex array computation.

CUDA is a general-purpose parallel computing architecture introduced by NVIDIA that contains calling instructions and a GPU computing engine that enables the GPU to solve complex computational problems. Program code developed based on CUDA separates out the complex array processing part of the entire computational processing flow and calls the GPU to accelerate it. Each discrete block of computational requirements is further subdivided into a new set of instructions. Finally, all of the subdivided computational tasks are executed simultaneously on the GPU processor. This is the principle of efficient parallel computing GPU hardware acceleration.

The traditional method of improving Cholesky decomposition to generate random fields is good, but the computational efficiency is very low and there are few practical applications. In order to put the random field generated by this method into practical application, the random field generation efficiency needs to be improved significantly. In the random field generation method proposed in this study, the decomposition of the correlation matrix is calculated using MATLAB programming. Among them, the modified Cholesky decomposition can be found in MATLAB with the corresponding eigenvalue decomposition operation, while in the case of larger and more complex models with more coordinate data. The correlation matrix array generated by the length of the coordinate data is huge. The decomposition process of the correlation matrix array is extremely time-consuming. The cycle efficiency of simulation calculation is seriously affected. To address this difficulty, this study uses GPU acceleration technology to improve the simulation efficiency.

In the key step of correlation matrix decomposition, the efficiency of GPU parallel computation increases significantly when the matrix size gradually increases, and the efficiency of cyclic simulation increases dramatically. However, the entire process of correlation matrix decomposition remains in the hands of the CPU. In the 2D random field simulations carried out in this study, MATLAB was used for the conversion of the simple code. When proceeding the modified Cholesky decomposition, the gpuArray command is called to load the relevant matrix into the GPU. After completing the decomposition of the correlation matrix with efficient parallel operations in the GPU, the required eigenvector matrix and eigenvalue matrix are called back to the CPU for the next computation by the gather command. Table 1 compares the decomposition of the correlation matrix with different sizes when computers with CPU and GPU computation schemes are used, respectively. As can be seen, the advantage of GPU acceleration gradually emerges as the array of operational matrices becomes larger, and its computational speed is much faster than that of the CPU, which saves a lot of time when doing multiple round-robin operations to compute probability problems.

3. Simulation of the Red Layer in Central Yunnan

3.1. Geological Background

The tectonics of rock stratum in soft rock body refers to the spatial distribution and arrangement among the sedimentary layers in the rock stratum. For example, the thick-layered block-thin laminated structure rock body refers to the sandwich-like and interlayer-like structures of the spatial combination arrangement of each stratigraphic rock unit (which can be divided according to lithology, layer thickness, etc.). It can be seen that the original rock stratigraphic structure of the weak geotechnical body is closely related to the engineering rock structure type.

In this study, the stratigraphic structure characteristics of the red-layered weak rock body in the Chuxiong section of the Yunnan Central Water Diversion Project are analyzed as an example. Through on-site engineering geological investigation and geological data analysis, the spatial variation law of the geological structure of the red-layered soft rock is analyzed. In this study, the mudstone content and the spatial distribution and arrangement of other rock layers were selected as the analysis index. Based on this, the red-bedded soft rock stratigraphy was classified as mudstone-free type, interbedded mudstone type, interbedded type, mudstone-dominated type, and all-mudstone type. Obviously, the mud content of the rock strata of different tectonic types in the above order is gradually increasing. In a sense, the engineering geological properties of the red-bedded soft rocks will gradually diminish.

Figure 3 shows the mud-free type of chondrites, which is characterized by a rather heterogeneous spatial distribution of strata. The distribution of different rocks is random and loose, and there is no prominent anisotropy in a certain spatial direction. The main rocks are gravels, and the corresponding random field model can be considered an isotropic random field.

Figure 4 shows the structure of the intercalated mudstone-type stratigraphy, which shows the intercalation of soft mudstone layers between hard sandstone layers. There is obvious anisotropy in specific spatial directions and spatial variability along the strike of the interbedded layers. The anisotropic random field model can be used to simulate this type of interbedded mudstone formation.

Figure 5 shows the anisotropic characteristics of the spatial structure of interstratified strata. In the red layer region of central Yunnan, the interstratified strata are mainly quartz sandstone, calcareous siltstone, and calcareous mudstone and siltstone layered with each other, forming a spatial structure of one layer of different rock-soil bodies. In the spatial direction, the variation of interstratified strata is obvious, and the anisotropic random field is applicable to describe their geological and structural characteristics.

Figure 6 shows the geological and tectonic features of the all-mudstone-type red layer. The rock types of the all-mudstone red layer are mainly calcareous and siltstone mudstone. In terms of geological structure, there is no obvious anisotropic feature in the all-mudstone red layer. However, because of the large content of soft rocks, the deformation of the all-mudstone red layer is larger under the excavation and unloading conditions. Therefore, its engineering risk should be considered as a key issue.

3.2. Modeling and Analysis Process for Red-Layered Soft Rock Tunnels in Central Yunnan

In the water diversion project in central Yunnan, a section of horseshoe-shaped hydraulic tunnel was selected for the study. The burial depth of this section of hydraulic tunnel is about 600 m, and the surface topography is gently sloping. The numerical model is shown in Figure 7. The model has a length of 80 m, a height of 80 m, a bottom wall of 8.72 m, and a height of 9.43 m. The rock and tunnel areas are divided in the model, and a hexahedral structured mesh is used to divide the model mesh. The number of meshes in this numerical model is 7225, and the number of nodes is 14594.

This section of hydraulic tunnel is in the softly interbedded cavern section of the Lion Rock tuff. The initial ground stress in this section of the hydraulic tunnel can be found by the depth of burial and the weight of the rock mass. The fluctuation value of the cohesion of the rock strength parameter is suggested to be between 0.1 and 0.4 MPa, and the fluctuation value of the internal friction angle is suggested to be between 25 and 35 degrees. The site engineering geological survey report provides suggested values for other rock mechanics parameters, as shown in Table 2.

Based on the random field theory, the hydraulic tunnel excavation simulation requires a preassigned model. According to the working conditions, a numerical model of the tunnel excavation area as well as the surrounding rock is established with the grid be divided. The preprocessing of the dissected mesh is done by a written program. The coordinates of the grid centroids are exported, and the coordinate matrix is generated in MATLAB. Then, the correlation matrix is generated from the coordinate matrix, and the random field matrix is generated using the covariance matrix decomposition method. Based on the Mohr-Coulomb principal structure model, the attributes to be assigned are selected and the coefficient-based equations are obtained using quadratic curve fitting. The random field matrix is substituted through code programming control to generate an initial random assigned model of the cavity based on random field theory. Numerical simulation analysis is performed on the produced numerical model. The maximum displacement of the surrounding rock in the excavation area of the hydraulic tunnel and the corresponding coordinate positions are recorded for a variety of random field conditions. The whole process can be accomplished by executing operating system commands with MATLAB. The specific flow chart is shown in Figure 8.

4. Study on the Influence Law of Anisotropy with Field

4.1. Computing Conditions

To simplify the analysis, this study assumes that the tunnel excavation is completed in one go. When the numerical simulation of cavern excavation is implemented in FLAC software, if the material properties of the rock mass in the tunnel excavation area are directly changed to null, it may cause a sharp release of ground stress, resulting in instantaneous collapse of the cavern. In order to reflect the progressive excavation process of the tunnel, the elastic modulus of the rock body in the excavation area of the cavern chamber can be reduced to 5% of the initial one by gradual weakening. Then, the material properties are changed to null, and the numerical simulation of the tunnel excavation process is finally completed. By the method above, the excavation numerical simulation of anisotropic random field cavern chamber with different rotation angles is done. The simulation results are analyzed as follows.

4.2. Simulation Results and Analysis

Figure 9(a) shows the enlarged view of the local model when the random field is horizontally anisotropic. After the excavation, the deformation results of the cavern are shown in Figure 9(b). It can be seen that the deformation of the right side wall of this hydraulic tunnel is larger than that of the left side wall under the condition of horizontal anisotropic random field. The maximum deformation is 81.3 mm, which is located in the middle of the right side wall, about 3.6 m from the bottom edge. The random field parameter at this location is the minimum value of the cavern envelope. The main stress distribution of the cavern chamber is shown in Figures 9(c) and 9(d). Due to the excavation, a significant ground stress release was generated in the critical area of the tunnel envelope. The major principal stress decreased from an initial 14.50 MPa to 0.34 MPa. The minor principal stress decreased from 14.50 MPa (compressive stress) to a tensile stress of 0.02 MPa. Due to the variability of mechanical parameters caused by the rock random field, the stresses in the cavern envelope exhibit an asymmetric distribution. However, the general distribution pattern is similar. The further away from the tunnel excavation critical surface, the less the stress release and shows an irregular circular divergence.

An enlarged local model of the anisotropic random field with 15 degrees of rotation is shown in Figure 10(a). After the end of excavation, the deformation of the tunnel is shown in Figure 10(b). It can be seen that in the anisotropic random field with 15 degrees of rotation, the deformation of the tunnel surrounding rock shows a phenomenon of large right and small left. The maximum deformation is 85.7 m, which is located in the middle of the right side wall, about 4 m from the bottom edge. A strip of rock with weak parameters is distributed at this location. The main stress distribution of the cavern is shown in Figures 10(c) and 10(d). It can be seen that the surrounding rock of the tunnel has a significant stress release due to excavation. The major principal stress decreases from 14.50 MPa to 0.19 MPa, and the minor principal stress decreases from 14.50 MPa in compression to 0.02 MPa in tension, while the stress distribution around the chamber is asymmetric due to the randomness of the mechanical parameters. The values of the random mechanical parameters of the rock mass in the direction to the left of the vault are larger. Therefore, the large principal stress release in the internal surrounding rock at this location is not very obvious compared with that in the surrounding rock at other locations.

The local model enlargement of the anisotropic random field with 30 degrees of rotation is shown in Figure 11(a). After the end of excavation, the cavern excavation deformation is shown in Figure 11(b). In general, it shows a large deformation of the right side wall and a small deformation of the left side wall. The direction of the most value of displacement is orthogonal to the direction of random field anisotropy. The maximum deformation is 65.2 mm, which occurs in the middle of the sidewall. The mechanical parameters of the tunnel envelope are higher in this simulation; therefore, the maximum value of the displacement is smaller than that of the random field simulation in the direction of 15 degrees of rotation. The principal stress distribution in the tunnel is shown in Figures 11(c) and 11(d). The significant stress relief in the tunnel envelope due to excavation resulted in a decrease of the major principal stress from 14.50 MPa to 0.27 MPa. Due to the randomness of the parameters, the stresses in the tunnel surrounding rock show asymmetric distribution characteristics.

The enlarged local model of anisotropic random field with 45 degrees of rotation is shown in Figure 12(a). After excavation, the deformation of tunnel surrounding rock is shown in Figure 12(b). Under this condition, the tunnel deformation shows the distribution pattern that the left side wall is large while the right side wall is small. The maximum deformation occurs in the middle of the left sidewall with a value of 77.9 mm, and the main stress distribution in the chamber is shown in Figures 12(c) and 12(d). The excavation resulted in a significant release of the ground stress in the tunnel surrounding rock. The major principal stress decreased from 14.50 MPa to 0.44 MPa. The minor principal stress decreased from 14.50 MPa in compressed condition to 0.02 MPa in tensile stress. Although the tunnel shape is symmetrical, the randomness of the parameters makes the stresses in the tunnel envelope show asymmetric distribution.

An enlarged local model of the anisotropic random field rotated by 60 degrees is shown in Figure 13(a). Figure 13(b) demonstrates the deformation of the tunnel excavation. Overall, the deformation of the left side wall is larger than that of the right side wall. The maximum displacement is 76.2 mm, which occurs in the middle of the left side wall at about 3.6 m from the bottom edge. Figures 13(c) and 13(d) show the main stress distribution in the tunnel envelope. Due to the excavation, the tunnel surrounding rock has produced obvious stress release.

Figure 14(a) shows the local model of anisotropic random field with 75 degrees of rotation. Figure 14(b) shows the simulation results of the deformation of the tunnel envelope after excavation. In general, the displacement of the left side wall of the tunnel enclosure is larger than that of the right side wall. The distribution direction of the most value of displacement is orthogonal to the random field anisotropy direction, and the maximum deformation value of 81.6 mm appears in the direction of 165 degrees on the left side of the vault. The main stress distribution of the cavern chamber, after excavation, and the ground stress of the tunnel surrounding rock occurred a significant stress release. The simulation results of large and small principal stresses are shown in Figures 14(c) and 14(d). The major principal stress decreases from 14.50 MPa to 0.36 MPa and remains in the compressive state. The randomness of the parameters makes the stresses in the cavity surrounding rock show asymmetric distribution.

Figure 15(a) shows an enlarged view of the local model of the rotated 90 degrees of anisotropic random field. The deformation of the surrounding rock after tunnel excavation is shown in Figure 15(b). It can be seen that, in general, the surrounding rock of the tunnel shows a pattern that the deformation of the right side wall is larger than that of the left side wall. The distribution direction of the displacement maximum is orthogonal to the random field anisotropy direction. A strip of rock with weaker mechanical parameters appears on the right side of the tunnel envelope. This increases the displacement of the surrounding rock in this area. The maximum deformation of the tunnel envelope is 76.9 mm, which occurs in the middle of the right sidewall about 4 m from the bottom edge. Figures 15(c) and 15(d) show the main stress distribution of the cavern. It can be seen that there is a significant stress release in the tunnel surrounding rock due to excavation. The major principal stress decreases from 14.50 MPa to 0.30 MPa.

For comparison between random fields with different rotation angle, the information is collected and summarized in Table 3. It can be seen that, due to the variability of mechanical parameters caused, the stresses in the cavern envelope exhibit an asymmetric distribution. There is a significant stress release in the tunnel surrounding rock due to excavation in each simulation, and the general distribution pattern is similar. In addition, when the rotation angle is 45°, 60°, and 75°, the maximum deformation occurs in the middle of the left side wall and has a slightly larger major principal than other conditions.

5. Conclusion

In this study, numerical simulation of excavation deformation of hydraulic tunnels was carried out for the anisotropic random field of red-bedded soft rocks in central Yunnan. The Mohr-Coulomb model was chosen for the rock material, and the internal friction angle and cohesion were set as random parameters, while other parameters were defined as constants. The effects of different correlation functions and different correlation lengths on the cavern excavation are compared in the isotropic random field model. The deformation characteristics and stress distribution laws of the anisotropic random field with different rotation angles are also compared.

The stress distribution patterns of anisotropic random field with different rotation angles are compared, because the randomness of the parameters causes the displacement and stress field of the cavern surrounding rock after excavation to show an asymmetric distribution law. The distribution law of maximum displacement is mainly related to the high and low values of the random parameters of the excavated cavern envelope. The maximum displacement value generally appears at the lowest parameter value. For the deeply buried tunnel in this study, the horizontal ground stress is large. The maximum displacement of the tunnel envelope appears at the side walls on both sides and not at the vault position. At the same time, there is an obvious stress release in the tunnel envelope. The excavation-induced ground stress release is weaker in the areas where the stochastic parameters are taken to be higher than elsewhere. Large deformation also does not occur in these domains.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Acknowledgments

This research is supported by the Fundamental Research Funds for the Central Universities (B200201059) and the Natural Science Foundation of China (Grant Nos. 51709089, 51939004, 12062026, and 11772116).