An iterative algorithm for random upper bound kinematical analysis

Stochastic analyses are currently used extensively in a wide range of geotechnical applications. This is the reason of soil spatial variability caused by the geological processes that form the soil layer (e.g., Ferreira et al., 2015). The influence of soil spatial variability is examined via numerical methods based on the finite element method, finite difference method, or finite limit analysis. These methods are used to examine the impact of soil strength spatial variability on engineering structures such as foundations, slopes, or sheet pile walls (e.g., Griffiths et al., 2002; Fenton and Griffiths, 2002, 2008; Griffiths and Fenton, 2004; Srivastava et al., 2010; Kasama and Whittle, 2011; Huang et al., 2013; Simoes et al., 2014; Zhu et al, 2017; Kawa et al., 2019; Halder and Chakraborty, 2020; Pramanik et al., 2020; Żyliński et al., 2020; Kawa and Puła 2020; Viviescas et al., 2020; Pramanik et al., 2021). However, the numerical efficiency of most existing probabilistic methods is not sufficient for practical applications where three-dimensional analyses are crucial. This is true especially for three-dimensional issues for which the three-dimensional nature of soil strength spatial variability suggests that there are difficulties with its two-dimensional simplifications (e.g. plane strain conditions). This is the reason researchers investigated the possibility of providing more efficient 3-D approaches for bearing capacity analysis in the case of soil spatial variability, e.g., Chwała (2019; 2020), Li et al. (2021), Li et al., (2021). In the earlier papers by Puła and Chwała (2018) and Chwała (2019), the algorithm for using the upper bound limit theorem and Vanmarcke's spatial averaging (Vanmarcke, 1977a,b; 1983) is presented for two-dimensional and three-dimensional bearing capacity evaluation, respectively. However, the numerical efficiency achieved by Chwała (2019) can be further improved by using the so-called constant covariance matrix approach (Chwała, 2020). The constant covariance matrix approach means that the matrix is obtained using the expected values of soil strength parameters; the concept was first proposed by Puła (2004; 2007). Generally speaking, the above approaches are based on random field discretisation to a single random variable (each random variable corresponds to a specific dissipation region in a failure mechanism). These variables are correlated with the covariance matrix, and this matrix is the basis for generating the average soil strength parameters (e.g., Fenton and Griffiths, 2008; Puła and Chwała, 2015). However, the solution suggested by Chwała (2019) for three-dimensional issues is based on one iteration of the covariance matrix (for more details see the next section). The same approach is used for two-dimensional analyses by Puła and Chwała (2018). Therefore, there is still a need of verify this assumption. To study the impact of the iteration number of covariance matrix and to extend the use of the proposed method to other geotechnical applications, an iterative algorithm was developed and is described in this study. The main objective is to provide a general procedure that can be utilised for a variety of geotechnical applications, where the upper bound limit theorem (e.g., Shield and Drucker, 1953) and Vanmarcke's spatial averaging are used. Another important objective of this study is to compare the simplified approach, which is based on the constant covariance matrix with the iterative procedure that is proposed here. The iterative procedure is based on redetermining the covariance matrix in the subsequent iterations; thus, the covariance matrix corresponds to the failure geometry. As a numerical example, the proposed algorithm is used to analyse the three-dimensional bearing capacity of a shallow foundation under undrained conditions and the bearing capacity of two-layered soils. The comparison of the iterative approach with the constant matrix approach is crucial for further practical applications of the simplified algorithm for three-dimensional issues like soil sounding location optimization (Chwała, 2020; 2021). This comparison is important due to the enormous reduction in computation time. Namely, three-dimensional analyses that were performed in this study for a constant covariance matrix took about 0.5 s using a standard notebook. Therefore, there is a great potential to improve the numerical efficiency in the three-dimensional probabilistic analyses by using the constant covariance matrix approach. However, the iterative procedure is needed to examine the impact of this assumption on the resulting bearing capacity probabilistic characteristics.

The proposed procedure is devoted to random analyses of geotechnical problems in the framework of the upper bound limit theorem. Its engineering applications are directed toward foundation bearing capacity evaluation and slope stability evaluation for spatialy variable soil. The procedure is described in as general manner as possible; however, some information directed to the examples used in the study are given. Two problems are investigated employing the proposed iterative algorithm, i.e., the algorithm proposed by Chwała (2019) for 3-D bearing capacity analysis, and the algorithm proposed by Chwała and Puła (2020) and Chwała and Kawa (2021) for two-layered soil bearing capacity calculation. As mentioned in the previous section, the procedure is based on the determination of the covariance matrix, which corresponds to the actual geometry of failure. The components of the covariance matrix are established via Vanmarcke spatial averaging, and the necessary formulae have to be derived independently for each problem (those for the three-dimensional bearing capacity are given by Chwała (2019), and for two-layered soil by Chwała and Puła (2020)). As indicated earlier, the above-mentioned algorithms use one iteration of the covariance matrix, meaning that after determination of the covariance matrix (which is based on the failure geometry that was obtained for uncorrelated soil strength parameters, see details in the algorithm description below), the correlated soil strength parameters are generated. Then, they are used to find the optimal failure geometry and the corresponding final bearing capacity. However, there is an inconsistency in the above approach, which is that the optimal failure geometry differs from the failure geometry for which the covariance matrix was established. This may introduce additional uncertainty in the observed results, which can be eliminated by the iterative procedure that is proposed in this study. In the algorithm detailed below, the covariance matrix is determined as many times as the assumed number of iterations N. The difference between the subsequent covariance matrices decreases when the number of iterations increases. To ensure the generalisability of the algorithm (i.e., the possibility that it can be used for various geotechnical problems), its layout is as universal as possible. For further discussion, the following algorithm is denoted as ‘A’.

Step 1. Preliminary assumptions. At the beginning stage, information about the initial random field characteristics of the soil spatial variability, is required. This information consists of a mean value μ, the point variance σ², the type of probability density function, and the scale of fluctuation θ (or scales of fluctuation if they are distinguished in the vertical and horizontal directions). Moreover, the random field correlation structure is assumed; the most commonly used covariance functions are Gaussian, see Eq. (1), and Markov, see Eq. (2) (see Fenton and Griffiths, 2008). The stationary and separable random fields are assumed in the algorithm. (1) R(Δx,Δy,Δz)=σ2exp{ −[ (πΔxθx)2+(πΔyθy)2+(πΔzθz)2 ] } R\left( {\Delta x,\Delta y,\Delta z} \right) = {\sigma ^2}\exp \left\{ { - \left[ {{{\left( {{{\sqrt \pi \Delta x} \over {{\theta _x}}}} \right)}^2} + {{\left( {{{\sqrt \pi \Delta y} \over {{\theta _y}}}} \right)}^2} + {{\left( {{{\sqrt \pi \Delta z} \over {{\theta _z}}}} \right)}^2}} \right]} \right\} (2) R(Δx,Δy,Δz)=σ2exp{ −[ 2| Δx |θx+2| Δy |θy+2| Δz |θz ] } R\left( {\Delta x,\Delta y,\Delta z} \right) = {\sigma ^2}\exp \left\{ { - \left[ {{{2\left| {\Delta x} \right|} \over {{\theta _x}}} + {{2\left| {\Delta y} \right|} \over {{\theta _y}}} + {{2\left| {\Delta z} \right|} \over {{\theta _z}}}} \right]} \right\} where, θ_x, θ_y, and θ_z are scales of fluctuation in the specified directions, and Δx, Δy, Δz are distances along x, y and z directions. If two soil strength parameters are described by random fields (e.g. angle of internal friction and cohesion), the correlation ratio ρ between them can be included in the analyses (in this case, the size of the covariance matrix from step 4 is enlarged; for more details see Puła and Chwała, 2015). Next, the corresponding failure mechanism has to be chosen or established. The failure geometry has to be kinematically admissible according to the upper bound theorem.

Step 2. Generation of independent soil strength parameters. According to the initial random field characteristics (μ, σ₂, and the type of probability density function), an appropriate number of independent soil strength parameters γ₀=γⁱ₀ are generated using a random number generator. This number is equal to the number of dissipation regions in the considered failure mechanism n (i=1,…,n). If two soil strength parameters are considered, both values are generated for each dissipation region.

Step 3. Determination of the optimal failure geometry. For soil parameters γ₀, the value of the corresponding upper bound limit load p₀ (bearing capacity) is computed based on the total energy dissipation in the failure mechanism. However, it is necessary to find the minimum upper bound limit load, and the optimal failure geometry is searched using a dedicated optimisation procedure. Thus, the optimal failure geometry Δ₀ and the corresponding upper bound limit load p₀ are determined. The optimisation procedure should be chosen individually for the considered problem, and there are no specific restrictions on this issue.

Step 4. Determination of covariance matrix and generation of correlated soil strength parameters. The covariance matrix C₀ is determined based on the optimal failure geometry Δ₀. The size of the matrix is n × n (or 2n × 2n if two correlated soil strength parameters are considered), and the matrix is symmetrical and positive definite (thus, the Cholesky decomposition exists, see e.g., Horn, 1985). The components of the covariance matrix between two dissipation regions can be derived from the following formula (for more details see Puła 2004; 2007 or Puła and Chwała, 2015): (3) Cov(Xi,Xj)=1| V1 || V2 |∫V1∫V2R(xi,yi,zi,xj,yj,zj)dVi(xi,yi,zi)dVj(xj,yj,zj) \matrix{{{\rm{Cov}}\left( {{X_i},{X_j}} \right) = {1 \over {\left| {{V_1}} \right|\left| {{V_2}} \right|}} {\int_{{V_1}}} {\int_{{V_2}}R\left( {{x_i},{y_i},{z_i},{x_j},{y_j},{z_j}} \right)} } \cr {d{V_i}\left( {{x_i},{y_i},{z_i}} \right)d{V_j}\left( {{x_j},{y_j},{z_j}} \right)} \cr } Where indexes i and j run from 1 to the number of dissipation regions n in the failure mechanism, and V₁ and V₂ are the integration regions (section, surface, or volume). Based on C₀, the independent soil strength parameter vector γ₀ is transformed into the correlated parameter vector γ^C₁ by calculating the product of standardized γ₀ vector and triangular matrix which is the result of Cholesky decomposition of the covariance matrix C₀ (for more details, see Puła and Chwała, 2015 or Chwała 2019).

Step 5. Beginning of the iterative procedure. Set k=1, and proceed to step 6.

Step 6. Set k=k+1. Using the optimisation procedure from step 3 for the correlated soil parameters γ^C_k, the optimal failure geometry Δ_k and corresponding upper bound limit load p_k are calculated.

Step 7. Determine the covariance matrix C_k for the optimal failure geometry Δ_k. Based on C_k, the independent soil strength parameter γ₀ is transformed into the correlated parameter γ^C_k+1 using the covariance matrix C_k.

Step 8. Using the optimisation procedure from step 3 for the soil parameters γ^C_k+1, the optimal failure geometry Δ_k+1 and corresponding upper bound limit load p_k+1 are calculated.

Step 9. If k≥k_max, proceed to step 10; otherwise, proceed to step 6. The maximum number of iterations k_max can be assumed arbitrarily after investigation of the convergence of the results. Additionally, different conditions for stopping the procedure are possible.

Step 10. For the end of the numerical procedure, the last values of the vector Δ_kmax and p_kmax are the optimal failure geometry and the corresponding upper bound limit load, respectively.

First, the above steps of the algorithm can be developed and adjusted in an individual pattern to be suitable for specified issues that are examined, such as slope stability or foundation bearing capacity.

Steps 2 to 10 are repeated N times, where N denotes the number of Monte Carlo realisations. Generally, the choice of the number N determines the accuracy that is obtained, i.e., a greater N provides more accurate estimations. Thus, a reasonable compromise between accuracy and computation time must be established. The overall computational efficiency can be improved using a framework that is characterised by faster convergence, such as subset simulation (Au and Beck, 2001).

As discussed above, the difference between subsequent covariance matrices C_k and C_k+1 is expected to decrease when the iteration number k increases. This suggests that there is a decreasing tendency in the differences between subsequent bearing capacities p_k and p_k+1 and failure geometry parameters Δ_k,i and Δ_k+1,i. The tendency can be expressed by Eq. (4). (4) limk→∞pkpk+1=limk→∞Δk,iΔk+1,i=1.0 \mathop {\lim }\limits_{k \to \infty } {{{p_k}} \over {{p_{k + 1}}}} = \mathop {\lim }\limits_{k \to \infty } {{{\Delta _{k,i}}} \over {{\Delta _{k + 1,i}}}} = 1.0 The values resulting from Eq. (4) are influenced by numerical errors that are caused by finite numerical accuracy and the procedure that ensuring the positivity of the covariance matrix. However, a specified iteration number k can be found, for which the difference between two subsequent values from Eq. (4) becomes negligible or equal to the numerical accuracy. Moreover, due to the results described in the next sections, the stabilisation is observed for a relatively small value of k. An increase in the iteration number k_max causes a proportional decrease in computational efficiency. However, for some engineering problems, the differences between p₁ and p_k (where k→∞) should be relatively small or even negligible. Therefore, in such cases, the computational efficiency can be significantly improved using only one or two covariance matrix iterations. One covariance matrix can also be determined for all Monte Carlo simulations, and the natural choice is to determine it for the failure geometry, which corresponds to the expected values of the soil strength parameters. The idea was first proposed by Puła (2004) and discussed later for a two-dimensional Prandtl failure mechanism by Puła and Chwała (2015). In this study, a simple general algorithm for using a constant covariance matrix is proposed and used for three-dimensional bearing capacity evaluation and bearing capacity of two-layered soil. In the numerical example section, results for both approaches are discussed. For each case, the algorithm that is based on the constant covariance matrix may provide sufficient accuracy for engineering practice, and it provides a dramatic improvement in computational efficiency. However, a previous investigation is required if the constant covariance matrix is to be used in numerical analyses. To illustrate the algorithm layout, the flowcharts are shown in Fig. 1. The path denoted by ‘A’ represents the ten steps described above, and the path ‘B’ represents a slight modification of ‘A’, where only Step 3 differs. The difference is a distinct way to determine the optimal failure geometry, i.e., in the path ‘B’ the optimal failure geometry is established for the expected values of the soil strength parameters. Path ‘B’ can be utilised to speed up the convergence of the results when the iterative procedure is used. The distinct path ‘C’ is the algorithm in which a constant covariance matrix is utilised (the covariance matrix in Step 3 in the path ‘C’ is determined only once for all Monte Carlo simulations).

Figure 1

The results obtained in the two previous sections are for different applications of the random failure mechanism method (RFMM). The first concerns the bearing capacity of two-layered soil when the failure mechanism is based on two-dimensional simplification, and the second concerns the bearing capacity of the rectangular foundation in the fully three-dimensional case. In both cases, soil spatial variability is considered three-dimensional. According to the obtained results, the differences between the approach with constant covariance matrix ‘C’ (k=0) and the iterative approach ‘A’ (k=1, k=2, k=3, k=4, k=5, and k=6) are generally in a narrow range. As expected, the greatest differences are observed between the constant covariance matrix approach and the iterative one for the first two iterations of the covariance matrix. However, despite the stabilisation of bearing capacity mean values and standard deviation for two and more iterations (k), some fluctuations can be observed, as shown in Fig. 4, Fig. 5, Fig. 7, and Fig. 8. The range of this fluctuation is not important for this study and does not impact the conclusions that can be drawn based on the shown examples. In the author's opinion there are two main sources for this fluctuation in the case of iteration indexes greater than 2, the first possibility is the numerical accuracy, especially in the calculation of high-order integrals, which is calculated to obtain the covariance matrix, another issue is the accuracy of the optimization procedure, which finds approximation to the global minimum of the objective function, not the exact global minimum. The second source can be caused by small modifications of the covariance matrix that are introduced to ensure the positive definiteness of the matrix (the modification occurs if the matrix is not positive definite). In the author's opinion, those factors combined may cause the observed small fluctuations in bearing capacity, mean values and standard deviation. The final accuracy can be influenced also, by the way, the described algorithms work. It can be observed especially for the longer foundations for which the averaging region can be significantly greater than the scale of fluctuation. Therefore, one additional assumption that is not discussed earlier has to be introduced here, that the foundation is rigid. Consequently, it is assumed by default in the numerical examples described above that the wall rests on the foundation and provides infinite foundation stiffness. As a result, the soil strength parameters can be determined as averages of large regions. Certainly, all above mentioned uncerntenties may overlap in the iterative procedure. In Fig. 9 a dependency of bearing capacity mean values, standard deviations and variation coefficient on undrained shear strength coefficinat of variation are shown. The case of two-layered soil is considered. To maintain the clarity of the plot only results for constant covariance matrix (k=0) and the last result for iterative procedure (k=6) are shown in Fig. 9. Three values of c_u coefficient of variation are analysed; namely, v_cu=0.25, v_cu=0.5 and v_cu=1.0. A very slow increase in the differences between mean values and standard deviations of bearing capacity in both cases (k=0 and k=6) is observed with the increase of v_cu (see Fig. 9a and Fig. 9b). Despite those differances, the resulting bearing capacity coefficient of variation is almost the same for each considered v_cu (see Fig. 9c).

Figure 9

All numerical results shown in this study were performed for relatively small number of Monte Carlo realisations. However, as mentioned earlier in this study the main goal is to examine the impact of using more iterations of the covariance matrix. For this particular approach the accurate determination of bearing capacity mean value and standard deviations is not necessary. Despite this, in Fig. 10 an example of such analyses is shown to demonstrate that even relatively low number of Monte Carlo realisations (N = 200) provides quite good estimation of bearing capacity mean values and standard deviations. The conclusions from the study are presented in the next section.

Figure 10

The study propose an original iterative algorithm (path ‘A’ or ‘B’ in Fig. 2) for the upper bound analyses of geotechnical problems with consideration of soil strength spatial variability. The resulting limit loads are realisations of a random variable for which probabilistic characteristics must be determined using the Monte Carlo method. The use of the iterative algorithm ensures the consistency between the failure geometry and the covariance matrix, from which the average soil strength parameters are determined - the averaging procedure is following Vanmarcke's spatial averaging (Vanmarcke, 1977a, 1977b, 1983). Moreover, the general algorithm for using a constant covariance matrix is presented in the study. Performing an analysis via the algorithm based on the constant covariance matrix (see algorithm ‘C’ in Fig. 2) provides a dramatic improvement in computational efficiency. Therefore, the possibility of using a constant covariance matrix is very promising for practical applications and three-dimensional analyses. The main objective of the study is to investigate weather these algorithms give similar results. According to the delivered examples it is shown that, the constant covariance matrix approach provides satisfactory results. However, the investigation required a comparison of the results provided by both approaches. This was done through the iterative algorithm proposed in this study. However, despite good agreement between the two considered approaches for the analysesd scenarios, i.e., bearing capacity of rectangular foundation and two-layered soil under plane strain assumption, other geotechnical applications need to be verified individually. Therefore, the conclusions are only limited to those two scenarios.

The proposed iterative algorithm was utilised in this study to evaluate the bearing capacity of shallow foundations. Two cases were considered, i.e., the bearing capacity of two-layered soil, and the three-dimensional analyses, both were performed for a variety of foundation shapes and scales of fluctuation. To summarize, the following conclusions can be drawn:

The study presents an original iterative algorithm for upper bound analyses of geotechnical problems with the inclusion of soil strength spatial variability. The proposed algorithm can be applied not only to bearing capacity evaluation, as shown in the study, but also for other applications like slope stability or retaining wall analyses. The use of the iterative procedure ensures consistency between the failure geometry and the covariance matrix; therefore, the algorithm allows the recognition and control of the impact of this issue on the resulting estimates. Basing the iterative algorithm on the upper bound theorem provides the opportunity to utilise it as a reference for other probabilistic methods.