A Prediction Method Based on Monte Carlo Simulations for Finite Element Analysis of Soil Medium considering Spatial Variability in Soil Parameters
With the Stochastic Finite Element Method (SFEM), the spatial variability of soil properties can be incorporated into the analysis of geotechnical structures. Although this method is significantly superior in principle to the homogeneous analysis of soil parameters, generalizing the method in engineering practice is difficult due to its computational inefficiency. In this paper, we propose a new method for the fast calculation of convergence results. The proposed method introduces a distance space to the Monte Carlo Method (MCM) random field instances and, considering the importance of a safety margin in structures, uses selected spatial interpolation to predict the MCM instances to be solved. Two case study simulations are presented. The results show that compared to the full Monte Carlo Simulation, the fast calculation method proposed in this paper can achieve very accurate convergence results while substantially reducing the computational cost, and the simulation errors for the structure are on the safer side.
The spatial variability of soil properties has been a major aspect of geotechnical design and investigation in recent decades. Phoon and Kulhawy  pointed out that geotechnical variability is derived from three aspects: inherent variability, measurement error, and transformation uncertainty. Elkateb et al.  noted that inherent variability is generated by aspects such as the parent material properties and the weathering erosion of a rock mass. Therefore, soil properties do not vary randomly in space, rather such variation is gradual and can be viewed as a superposition of smoothly varying trends and fluctuating components. As traditional deterministic analyses do not account for spatial variability in soil properties, the corresponding results may not reflect the true failure mechanisms and response of geostructures to loading . Vanmarcke , the first person to propose a random field model of soil profiles, used the numerical characteristics of the random field to describe the natural variability of soil properties. After that, Vanmarcke et al. [5, 6], Kiureghian and Ke , Spanos and Ghanem , Degroot and Baecher , and Zhang and Ellingwood  conducted in-depth investigations of the spatial variability of soil properties and improved the random field theory.
The Stochastic Finite Element Method (SFEM) combines random field analysis, the finite element method, and the Monte Carlo Method (MCM) . In recent years, the SFEM has been widely used to study the impact of the spatial variability of soil properties on geotechnical structures. El-Ramly et al.  and Srivastava and Babu  studied the effect of the spatial variability of soil parameters on slope stability and deformation; Griffiths and Fenton  and Jimene and Sitar  analysed foundation settlement; Huang et al.  investigated the differential settlement of tunnels given the spatial variability of soil. Wang and Yeh , Zakian et al. [18, 19], and Yeh and Rahman  studied seismic responses considering soil spatial variability. Zakian and Khaji [21, 22] used a stochastic spectral finite element method to investigate wave propagation in a stochastic medium. The MCM is both simple and direct. In theory, this method can be used to analyse the spatial variability of geotechnical structures; in practice, however, Monte Carlo Simulation (MCS) has high computational costs when used to analyse complex structural systems. Huang et al.  employed the SFEM to conduct calculations for a specific engineering case and noted that at least 600 MCM instances are required to obtain convergence results for the probability characteristics. Cho  used the Latin hypercube sampling technique to generate random soil properties; compared to simple random sampling, Latin hypercube sampling reduces the number of calculations, but the number of calculations required by the fixed value method is still large. Ghanem and Spanos  proposed the spectral expansion method of nodal random variables, but for complex structural systems, this method involves many equations and involves difficult calculations. To reduce the computational cost of finite element analysis using Monte Carlo Simulation, Dyson and Tolooiyan  introduced a method for random field similarity measures to predict the convergence value of the factor of safety by comparing the similarity between MCM random field instances, which greatly reduced the computational cost, and the proposed method was compared with the full MCS to verify its effectiveness.
On the basis of Dyson and Tolooiyan’s idea, a new method for predicting MCM random field instances based on spatial interpolation is proposed in this paper. This method is expected to popularize the analysis of the spatial variability of soil in practical engineering.
2. Stochastic Finite Element Method
2.1. Numerical Characteristics of Random Fields
The random field model proposed by Vanmarcke  essentially uses a continuous stationary random field to simulate the soil profile. The core idea of this theory is to use the variance reduction function to combine the point variability and spatial variability in a soil body to transform the point properties of the soil into the average spatial properties. Indicators that describe the spatial variability in soil parameters include the mean μ, standard deviation s, and scale of fluctuation δ. The coefficient of variation (COV) is defined as the ratio of s to μ (formula (1)), which can be used to measure the inherent geotechnical variability:
The scale of fluctuation δ is an intrinsic soil parameter that indicates the soil parameters are basically correlated within this range and essentially uncorrelated outside this range. Therefore, the scale of fluctuation δ reflects the correlation of soil parameters: small δ indicates large variation and weak correlation of the soil parameters while large δ indicates smooth variation and strong correlation of the soil properties. Discretization of a random field requires that the autocorrelation function (ACF) be determined. The ACF is used to describe the magnitude of the correlation between points in space , and there are a number of alternative correlation functions [27, 28], such as the Markovian spatial correlation function, the square exponential correlation function, the second-order autoregressive model, and the cosine exponential correlation function. Considering the generality and simplicity of the Markovian spatial correlation function , this function (formula (2)) is used as the ACF in this study to discretize the random field. The proposed method in this paper is also applicable to other types of correlation functions:where δx, δy, and δz are the scales of fluctuation in the x, y, and z directions, respectively; Δx, Δy, and Δz are the distances between the elements along the x, y, and z directions, respectively.
Related studies [2, 29] have shown that the soil parameters in nature mostly follow a normal or lognormal distribution, and lognormal distributions are more common because they avoid the generation of negative values. Similar to the existing random field generation method , the lognormal random field is obtained by first simulating a normal random field. The lognormal distribution of E means that ln E is normally distributed. Therefore, the standard deviation s and mean u of the normal distribution of lnE are given by
2.2. Random Field Generation Method
With reference to studies by other researchers [30–32], the generation of random fields in the present study is based on covariance matrix decomposition. The finite difference software FLAC3D is employed in this research. FLAC3D provides an embedded FISH language, which can be used for the development of FLAC3D. The FISH language is used to generate random fields, and the specific steps are as follows:(1)Build the numerical model with n elements and obtain the coordinates of the centre point of each element.(2)Use formula (2) to calculate of the ith element relative to all elements (including the ith element itself) and obtain an n-order column vector. Obtain n column vectors by varying i from 1 to n, and combine them using i as the row number to form a matrix , which is the autocorrelation covariance matrix of the model:(3)Decompose by Cholesky decomposition:(4)Y is a randomly generated n-dimensional column vector, and the components of which are independent of one another and follow the standard normal distribution as a whole:(5)Allocate Z to corresponding elements in the numerical model to generate a random field model in FLAC.
Repeat the previous steps M times to generate all the required MCM random field instances.
3. Existing Methods
To shorten the SFEM computation time, Dyson and Tolooiyan  proposed a method to predict individual MCM random field instances instead of carrying out the conventional finite element simulation. This method is based on a similarity measurement technique (the Frobenius norm of the difference between the two matrices). Dyson and Tolooiyan used case studies to illustrate that similar random fields produce similar simulation results. According to this phenomenon, the results of each set of MCM instances to be calculated can be predicted by simply simulating the first 50 sets of MCM instances and then making similarity judgements between the remaining MCM random field instances to be calculated and the first 50 sets.
Dyson and Tolooiyan’s algorithm greatly reduces the computational cost of the MCM, and the prediction results are close to those of full MCS. Nevertheless, there is still room for improvement. First, the highest similarity between the two sets of random fields does not necessarily mean that the two sets of MCM instances have the closest results. There could be other methods for further improving the prediction accuracy. Second, it is not clear if the convergence results obtained by prediction methods reduce the structural safety margin. This issue determines if it is safe and reliable to apply the prediction results to actual engineering projects and poses the most concern to engineers. Based on these two considerations, we improved the specific prediction method and proposed a new method for the fast calculation of convergence results in models considering the spatial variability of soil.
4. New Fast Calculation Method
4.1. Definition of Distance
The fast calculation method proposed in this paper is based on a distance space. The normal random vector Z generated by formula (6) is a dimensional column vector, for which the distance is provided aswhere , .
It is easy to see the following:(1).(2)The necessary and sufficient condition of is .(3).(4)Suppose , , and are three points in the n-dimensional distance space defined above. Physically, these are three sets of randomly generated parameters and the corresponding spatial states. By using the Cauchy inequality, it can be proved that satisfies the triangle inequality, that is, .
The above observations show that the normally distributed random vector Z is an n-dimensional distance space, denoted as , of the provided distance .
4.2. Distribution Characteristics
The origin is defined in the distance space :
To investigate the distribution of MCM instances, SFEM analyses were performed on two case study simulations based on the aforementioned theory. The scatter plot showed that there was an obvious banded distribution. Therefore, a pair of squeeze parallel lines was set to envelop the points. The method is described as follows:(1)Perform calculations for the first N sets of MCM instances, and plot the scatter plot, with the distance from the MCS instances to the origin as the horizontal coordinate and the simulation result of the MCM instances as the vertical coordinate.(2)Let two parallel lines be and such that the points appear in the region between the squeeze parallel lines.(3)Solve for , , and as follows. First, set the possible range of values for ; in this case, set . Second, solve for and to minimize the distance between the two parallel lines when . Third, find , , and values that correspond to the minimum value of as takes each integer between −2000 and 2000;(4)To obtain stable squeeze parallel lines, some boundary points are allowed to be excluded. It is worth noting that each point in the figure represents a set of random field simulation results. To reduce the distortion caused by the excluded boundary points, the exclusion of at most N×2% boundary points is recommended to obtain the corrected , , and .
Figure 1 is the plot of the actual calculation results of the one set of model processes described above, with numbers of known sets, N, of 30, 50, 75, 100, 200, 300, 400, and 500. The figure shows that when N is small, the slope and intercept of the squeeze parallel lines change greatly as N increases; when N ≥ 50, the slope and intercept of the squeeze parallel lines basically do not change as N increases, indicating that when N = 50, the range in which and possibly appear is similar to that of a full MCS. Considering the balance between the calculation economy and the prediction accuracy, the calculation results of the first 50 sets of MCM instances are used for a fast calculation of the convergence results of the full MCS.
The distance between the two sets of MCM instances, , being the minimum does not necessarily mean that the calculation results of the two sets are the closest; nevertheless, in general, there are usually close calculation results when the distance between sets of MCM instances is small. To ensure that the convergence result obtained by the prediction method is on the safe side, a selection mechanism is introduced, the principle of which is determined by the type of the target result; for example, larger displacement and stress simulation results or smaller structural resistance simulation results increase the structure’s safety margin.
Based on the aforementioned research, the prediction method in this study can be performed using the following procedure, which is summarized in Figure 2:(1)Select the first N sets of MCM instances as the known sets. For each known set of MCM instances, obtain its distance from the origin, , and its model calculation result, .(2)Stop calculating and instead predict the result for the (N + 1)th set and beyond. Supposing that the set y is the MCM instances to be predicted, calculate the distance between each known xth set of MCM instances and the yth set of MCM instances to be predicted. Find the five known sets with the smallest . According to the aforementioned selection mechanism, if the displacement and stress are to be predicted, retain the four sets with the largest values out of the five sets; if the structural resistance is to be predicted, retain the four sets with the smallest values out of the five sets. It should be noted that the use of five known groups for prediction is not required; the number of known groups can be adjusted according to the needs of the study. However, if the number of known groups is too small, the prediction results may contain a large deviation.(3)Predict the result of the yth set of MCM instances, , using the inverse distance weighted (IDW) method:(4)To exclude possible outliers from the prediction results, a discriminating mechanism is set up according to Section 4.2. More than 98% of the points are located between the squeeze parallel lines; the remaining less than 2% of the points that are outside the squeeze parallel lines generally deviate within 10% from the straight line on the same side. Therefore, if the result of the MCM instances predicted by the above algorithm deviates from the squeeze parallel lines by more than 10%, this prediction result is considered abnormal, and the prediction result is corrected to the value of the straight line on the same side under the same horizontal coordinate.(5)Perform a statistical analysis after the prediction.
5. The Convergence of a Tunnel
5.1. Numerical Modelling
To study the effect of spatial variability in the elastic moduli of geotechnical structures, a two-dimensional finite difference model is established using a metro shield tunnel in Shanghai as a prototype based on the case study of Huang et al. . The model has a free boundary at the top (ground surface), normal displacement constraint boundaries on all sides, and a fixed boundary at the bottom, as shown in Figure 3.
The numerical analysis software FLAC3D is used to establish the model with a total of 9344 elements, and D is taken as 6.2 m. To facilitate the calculations with the model, the model adopts the elastic constitutive law and only considers the condition of complete drying. Poisson’s ratio is 0.33, and the density is 1800 kg/m3. The elastic modulus is assumed to follow a lognormal distribution with a mean of 20 MPa, a of 0.3, and a scale of fluctuation of . Using formulas (3) and (4), the lognormal distribution with mean and variance is transformed to the standard normal distribution with mean and variance . Therefore, 600 sets of elastic moduli are generated using the method in Section 2.2. Figure 4 shows the spatial distribution of the elastic moduli generated by a set of typical random fields, with the different colors of the elements representing different values of elastic moduli.
The original in situ stress field is formed by calculating the equilibrium under the dead weight, and then the horizontal and vertical displacements, and , of the profile after cavern excavation are calculated (as shown in Figure 5).
5.2. Algorithm Performance
Figure 6 is a comparison of the convergence results obtained using the existing method , the proposed method, and the full MCS, using the number of sets of random fields as the horizontal coordinate and the average of the horizontal and vertical displacements as the vertical coordinate. Note that, for more than 600 instances, the mean value of the simulation results no longer fluctuates significantly as the number of instances increases, and the fluctuation range is only approximately 0.2%, indicating that the results have converged, which is consistent with the results of Huang et al. . The clearly converged results for 800 instances shown in the figure were taken as the full-MCS result.
A workstation equipped with an Intel (R) Core(TM) i9-7980XE CPU was used to calculate the numerical models of this case, and the calculation time of a single instance was approximately 19 minutes. The full MCS calculated 800 sets of instances for a total of approximately 253 hours. The method in this paper only needs to calculate 50 sets of instances, which takes approximately 16 hours. The running time of the prediction process is less than 1 minute. Therefore, compared with the full MCS, the proposed method reduces the computation time by approximately 93%.
In this case study, the convergence results obtained by the existing prediction method  and the proposed method are both larger than those obtained by the full MCS. Therefore, the safety margin of the structure is not reduced, and both of the simulation errors in the structure are on the safer side. Moreover, the convergence results of the method proposed in the present study are closer to those obtained by the full MCS. In the calculation of , the relative error between the convergence results of the proposed method and the convergence results of the full MCS is 0.39%, which is 66.7% less than the existing method. In the calculation of , the relative error between the convergence results of the proposed method and the convergence results of the full MCS is 0.55%, which is 8.3% less than the existing method. A detailed comparison is presented in Table 1.
6. The Face Stability of Tunnels
6.1. Numerical Modelling
When tunnelling, normal stress is usually applied to the working face as the supporting force, and the magnitude of this stress affects the displacement of the centre of the working face. There is a critical minimum supporting force that ensures the stability of the working face. Once the supporting force is less than this critical value, unstable displacement of the working face occurs. To investigate the magnitude of the critical supporting force, a tunnel model is established, as shown in Figure 7, in reference to the study of Cheng et al. . The shallow buried circular tunnel has a diameter of 10 m, a burial depth of 10 m, and lengths of 20 m, 35 m, and 30 m in the x, y, and z directions, respectively. Except for a free top surface, the rest of the boundaries are constrained in the normal direction. The structure is discretized into 22,152 elements.
The numerical model only considers dry conditions. The Mohr–Coulomb yield criterion commonly used in geotechnical engineering is used in the numerical model, and the basic parameters include an elastic modulus of 24 MPa, Poisson’s ratio of 0.3, a density of 1800 kg/m3, and a tensile strength of 1 kPa.
The original in situ stress field is established by calculating the equilibrium under the dead weight. After tunnel excavation, a normal supporting stress P is applied on the working face, and the displacement of the centre of the working face along the Y direction is calculated. The value of P is gradually decreased so that a series of P values and the corresponding can be obtained. The stress-displacement curve, with P as the horizontal coordinate and as the vertical coordinate, is plotted in Figure 8. After the supporting stress exceeds the minimum pressure point, the horizontal displacement increases sharply as the supporting stress P decreases, so the supporting stress corresponding to the minimum pressure point is considered as the minimum supporting stress.
In this case study, the effects of the spatial variability in the friction angle and the cohesion on the minimum supporting stress of the working face are investigated. Note that the calculation results of the minimum supporting force are related to the step size of the supporting stress in the calculation process. To minimize the effect of the step size on the results, the step size is set to a very small value of 40 Pa.
6.2. Cohesion Variation
The cohesion is assumed to follow a lognormal distribution with a mean of 10 kPa, a COV of 0.3, and a friction angle of 15° and with the rest of the basic parameters the same as above. The scale of fluctuation is , and 500 sets of random cohesion fields are generated using the method in Section 2.2. Figure 9 shows the typical spatial distribution of the cohesion formed by a set of random numbers, with the different colors of the elements representing different cohesion values.
6.3. Friction Angle Variation
The friction angle is assumed to follow a lognormal distribution with a mean of 15°, a COV of 0.3, and a cohesion of 10 kPa, with the rest of the basic parameters the same as above. The scale of fluctuation is , and 500 sets of random friction angle fields are generated using the method in Section 2.2. Figure 10 shows the typical spatial distribution of the friction angle formed by a set of random numbers, with the different colors of the elements representing different friction angle values.
6.4. Algorithm Performance
Figure 11 is a comparison of the convergence results obtained using the existing method , the proposed method, and the full MCS, where the number of sets of random fields is taken as the horizontal coordinate and the average of the minimum supporting force is taken as the vertical coordinate. As in the previous case, the calculated result for 800 sets of random fields was taken as the full-MCS result.
The same workstation was used to calculate the numerical models of this case, and the calculation time of a single instance was approximately 45 minutes. The full MCS calculated 800 sets of instances for a total of approximately 600 hours. The method in this paper only needs to calculate 50 sets of instances, which takes approximately 38 hours. The running time of the prediction process is less than 1 minute. Therefore, compared with the full MCS, the proposed method reduces the computation time by approximately 93%.
Figure 11 shows that under both conditions of this case study, the minimum supporting force obtained by the existing prediction method  is smaller than the convergence result of the full MCS. The application of this prediction result in engineering projects may cause accidents. The minimum supporting force obtained by the proposed method is higher than the convergence result of the full MCS, and the deviation is small. In the cohesion variation group, the relative error between the convergence results of the proposed method and the convergence results of the full MCS is 0.03%, which is 98.1% less than the existing method. In the friction angle variation group, the convergence result of the proposed method is the same as that of the full MCS. A detailed comparison is given in Table 2.
7. Discussion and Conclusion
A fast calculation method for SFEM convergence results based on spatial interpolation was presented in this study. Two case study simulations were provided. In the first case study, the horizontal and vertical displacements of an underground cavern were calculated considering the spatial variability of the soil elastic modulus. In the second case study, the minimum tunnel supporting force was calculated considering the spatial variability of the cohesion and the friction angle. A comparison with full MCS was carried out. In case study 1, the convergence results obtained by the existing prediction method and the proposed method do not reduce the safety margin of the structure, and the deviation due to the fact that the proposed method is smaller. In case study 2, the convergence results obtained by the existing prediction method for the structure deviate greatly towards the more unsafe side while the proposed method produces smaller errors, and the results deviate towards the safer side.
This method can be applied to the analysis of geotechnical structures considering spatial variability. Compared with the traditional full MCS, this method can obtain convergence results on the safer side under the premise of substantially reducing the calculation burden, thereby facilitating the application of spatial variability analyses of soil parameters in practical engineering.
As a final remark, this study only considered the effect of the spatial variability in the parameters of homogeneous soil. The analysis of multilayered soil considering spatial variability is a topic deserving further investigation.
Some or all data, models, or codes generated or used during the study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
E. H. Vanmarcke, “Probabilistic modeling of soil profiles,” Journal of the Geotechnical Engineering Division-Asce, vol. 103, no. 11, pp. 1227–1246, 1977.View at: Google Scholar
E. H. Vanmarcke, Random Fields: Analysis and Synthesis, MIT Press, Cambridge, England, 1983.
H. Huang, S. W. Gong, and L. ZhangJuang, “Simplified procedure for finite element analysis of the longitudinal performance of shield tunnels considering spatial soil variability in longitudinal direction,” Computers and Geotechnics, vol. 64, pp. 132–145, 2015.View at: Publisher Site | Google Scholar
Q. Cheng, S. X. Luo, and X. Q. Gao, “Analysis and discuss of calculation of scale of fluctuation using correlation function method,” Rock and Soil Mechanics, vol. 21, no. 03, pp. 281–283, 2000, In Chinese.View at: Google Scholar
G. A. Fenton and D. V. Griffiths, Risk Assessment in Geotechnical Engineering, John Wiley & Sons, Hoboken, NJ, USA, 2008.
G. B. Beacher and T. S. Ingra, “Stochastic FEM in settlement predictions,” Journal of the Geotechnical Engineering Division, vol. 107, no. 4, pp. 449–463, 1981.View at: Google Scholar
H. Z. Cheng, J. Chen, Z. F. Hu, and J. H. Huang, “Face stability analysis for a shield tunnel considering spatial variability of shear strength in sand,” Rock and Soil Mechanics, vol. 39, no. 8, pp. 3047–3054, 2018.View at: Google Scholar