1. Introduction

Bearing capacity evaluations of shallow foundations in the case of spatially variable soil are currently used primarily for two-dimensional analysis. One of the reasons for this situation is numerical efficiency, which limits the use of three-dimensional finite element formulations for random analyses. However, recent attempts have been made to solve such issues, e.g., [1]. Moreover, most existing probabilistic methods have been applied to single-layered soil. However, from an engineering perspective, two-layered soil is also a highly probable situation. An especially important situation is when the bottom layer is a soft soil and the top layer is a strong soil. This scenario has been extensively examined in the case of deterministic analyses (e.g., [2–5]); however, there have been very few studies on random analyses of two-layered soil (e.g., [6]). It is worth mentioning that it is possible to perform these types of two-dimensional analyses by using OptumCE software.

However, the lack of random approaches that are efficient and capable of introducing three-dimensional analyses was the motivation to propose an original approach that utilizes a kinematic approach in conjunction with Vanmarcke spatial averaging [7–9]. An analogous algorithm for plane strain conditions and a multiblock failure mechanism for a one-layered soil was proposed by Puła and Chwała [10]. The present study is also motivated by the engineering application to the problem of the bearing capacity of working platforms, which are usually constructed over soft subsoil to ensure safe traffic of heavy trucks. A working platform usually consists of compacted cohesionless material. Due to the controlled compaction procedure, the relatively homogenous soil layer is obtained. However, the spatial variability of the natural soil layer situated below the man-made layer can play a significant role in bearing capacity estimation. This situation is reflected in the numerical analyses. The importance of the problem can also be illustrated by investigating failures of working platforms, e.g., Channel Tunnel Rail Link Contract 310 [11] or Lublin (eastern Poland) road bypass (e.g., [12]).

In this study, as a deterministic background, a failure mechanism proposed by Michałowski and Shi [3] is assumed. The selected failure mechanism is dedicated for two-layered soil, i.e., the top layer is medium or dense sand, and the bottom layer is soft clay. According to the above described assumption, the spatial variability in soil strength is considered only in the bottom layer; thus, this approach characterizes the spatial variability in the undrained shear strength. The problem of designing working platforms is coded in BR 470 [13]. This study provides an efficient method for analysing such engineering problems and demonstrates that the random bearing capacity of two-layered soil in the case of spatial variability in the bottom layer can be efficiently evaluated by the developed kinematic approach. Moreover, despite the use of a two-dimensional mechanism, the spatial variability of the undrained shear strength in the soft clay layer is considered in three dimensions. As a result, an assessment of the impact of foundation length on random bearing capacity can be estimated. The importance of performing three-dimensional analyses is noted in this study. This is especially true if soil spatial variability is considered; thus, the spatial variability has a three-dimensional structure, and using two-dimensional simplifications can significantly influence the final results. A numerical algorithm is proposed in this study. According to this numerical procedure, numerous two-layered soil scenarios are analysed. Moreover, a simple approach to separate the impact of the thickness of the homogenous layer on random bearing capacity estimates from the impact of the spatial variability of the bottom layer and foundation length is provided in the study. The obtained results are shown and discussed in section 5.

In summary, the objective of this study was to propose an efficient approach for three-dimensional bearing capacity analysis of two-layered soil. The considered situation reflects the working platform scenario in which the top layer is a man-made layer of compacted cohesionless soil (in this study assumed to be homogenous) and the bottom layer is a natural, spatially variable, cohesive soil.

2. Algorithm

2.1. Failure mechanism description for random analyses

The deterministic failure mechanism proposed by Michałowski and Shi [3] is used in this study. The geometry of the failure mechanism was modified to be suitable for random analyses. As a result, a new formula for bearing capacity is introduced. A detailed description of the random failure mechanism is provided below. The first assumption concerns the homogeneity of the top soil layer (note that in this study, a sand layer is considered to be the top layer). In the case of a working platform, the upper layer is practically a man-made material. Therefore, the spatial variability of this layer is neglected, and only a deterministic case is provided. As a result, the spatial variability in soil properties is considered for the bottom layer (purely cohesive soil). Moreover, for the bearing capacity estimations performed in this study, the bottom layer is assumed to be weightless. In the case of undrained conditions, the weight of the soil can be neglected. Moreover, as shown in Chwała [14] and also in the case of random analyses, the consideration of soil weight has a very limited impact on the final bearing capacity estimations. As in the deterministic version of the failure mechanism [3], the random failure mechanism proposed in this study consists of rigid blocks of soil that move with respect to each other. This failure mechanism is constructed in a particular manner so that there is no velocity discontinuity along the line that separates the top soil layer from the bottom layer. For the numerical analyses performed in this study, the failure mechanism of the 11th block is assumed. An example geometry and the rigid block numbering convention are shown in Fig 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Rigid block numbering convention for the random version of the failure mechanism for two-layered soil.

Note that t is the thickness of the top layer (homogenous sand) and that the bottom layer is assumed to be spatially variable clay.

https://doi.org/10.1371/journal.pone.0231992.g001

Note that the block shown in Fig 1 is a single rigid block, i.e., the block denoted by O₁P₁A₁A₂P₂O₁ defines a rigid block. Due to assumption of the homogeneity of the sand layer (top layer), all velocity jump vectors within this layer are inclined to the slip lines at the same angle that is equal to the assumed internal friction angle of sand φ (this angle is shown in Fig 2). Note that the velocities shown in Fig 2 can be obtained from the velocity hodograph shown in Fig 2. The velocity hodograph is constructed according to the geometrical relations between the velocity jump vectors shown in Fig 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Detailed geometry of the failure mechanism.

Note that if spatial variability of the bottom layer is assumed, the failure mechanism is non-symmetrical.

https://doi.org/10.1371/journal.pone.0231992.g002

The first noticeable difference when comparing the failure mechanism shown in Fig 2 with the deterministic failure mechanism is that in the random case, the failure geometry is no longer symmetric. Due to the different soil strength properties for each velocity discontinuity, the failure geometry is asymmetrical with respect to the vertical line passing through the foundation centre. An example random failure geometry is shown in detail in Fig 2. Based on Fig 2 and the upper bound theorem, the formula for bearing capacity can be derived. As an illustrative example, some parts of the bearing capacity formula are derived below. Based on the upper bound theorem, the upper limit of the bearing capacity can be found by comparison with the dissipation energy along velocity discontinuities and gravitational forces. Note that no overburden pressure is considered (a foundation is placed on the soil surface, as shown in Figs 1 and 2). According to results presented by Chwała [14], the self-weight of soil in the case of purely cohesive soil can be neglected when the spatial variability of soil strength is considered. As an example, the energy dissipated along the velocity discontinuity line P₁A₁ can be calculated as shown in Eq (1). (1) where c₁ and l₁ are the cohesion and dissipation line length, respectively, and v₁ is the velocity jump vector. All three values are shown in Fig 2. By using an analogous method, the remaining formulas for energy dissipation within the bottom layer can be derived. Moreover, because of the cohesionless character of the top layer, there is no energy dissipation within it (c = 0). Therefore, in the considered failure mechanism, the energy dissipates along sixteen lines within the bottom layer (see Fig 2). Another factor that influences the final bearing capacity estimation is the power of the gravitational forces on the top layer. The formula for region O₁P₅C₁ is demonstrated in Eq (2). (2) where F₆ is the area of the rigid block indexed by 6 (see Fig 1). Note that in the case of block 5, the corresponding area F₅ is the area of block 5 within the top layer. The velocity vector in brackets [v₁₀] denotes the vertical component of the velocity jump vector v₁₀. The unit weight of the sand layer is denoted by γ. By using an analogous method, the remaining formulas for the power of gravitational forces within the top layer can be derived. The power of the limit load that the foundation can be subjected to can be calculated by multiplying the bearing capacity p by the velocity v₀ of the rigid block beneath the foundation (the foundation moves with the same velocity that is vertically oriented). Note that the formulas provided above assume by default that the calculations are performed for a = 1.0 m (see Fig 2); however, if other foundation lengths are considered, the corresponding value of a must be included. By summing all terms that relate to energy dissipation and gravitational forces and by transforming the resulting equation due to p, the following formula can be obtained: (3) where, as indicated earlier, v_i is the magnitude of the velocity jump vector along the slip line i, [v_i] denotes the vertical component of velocity v_i (negative if upward), and l_i is the length of the slip surface i. In Eq (3), the undrained shear strengths are denoted simply by c_i, where the index i is the same as that for the slip line lengths l_i (see Fig 2).

Design charts for calculating bearing capacity were provided in a paper by Michałowski and Shi [3]. As shown in Fig 2, when increasing the depth of the bottom layer (t), the critical value of t can be found. This critical value means that the failure mechanism is no longer deep (spread out in both soil layers) but becomes much shallower (spread out only in sand). This issue was described in detail by Michałowski and Shi [3], who showed that the critical depth of the weak layer depends on the angle of internal friction φ of the sand and the ratio c_u/γb. Nevertheless, this issue is not considered in this study; thus, the examined scenarios relate to the case in which a deep failure mechanism occurs (as shown in Figs 1 or 2). The investigated situation is adequate for working platforms in which the top sand layer is relatively thin.

The formula shown in Eq (3) is crucial for determining the bearing capacity in the case of spatial variability in the strength properties of the bottom layer by applying different undrained shear strengths on each dissipation region. However, it is essential to properly calculate these undrained shear strengths. The method used in this study for that purpose is detailed in the next section. Moreover, Eq (3) provides the bearing capacity if the failure geometry is known (see Fig 2). Therefore, an optimization procedure is necessary to find the minimum value of the limit load because the method is based on the upper bound approach.

2.2. Spatial averaging

As indicated in the previous section, the method for determining average undrained shear strengths is a crucial element of the proposed approach. Because the sand layer is assumed to be non-random, the spatial averaging is subject only to the random field X that describes the spatial variability of the bottom layer. For this purpose, the method proposed by Puła and Chwała [15] is used here. Note that the approach described in the mentioned paper concerns a one-layered soil and Prandtl failure mechanism. However, with some modifications and adaptations, this approach can be utilized for the scenario considered in this study. The modified version of the spatial averaging procedure is described below. Additionally, to maintain the clarity of this paper, some general information is also provided below.

The approach is based on Vanmarcke’s spatial averaging concept [7–9]. Generally, the procedure replaces the need to generate a random field X by generating a set of correlated random variables (random vector X_V). Those random variables (X_V,1, X_V,2, …, X_V,n) or random vector components correspond to the dissipation regions resulting from the considered failure mechanism. Due to the assumed stationarity of the random field, after the averaging procedure, the mean values of the random vector components are preserved; however, their variances are subject to reduction. The average is calculated as shown in Eq (4): (4) where V_i is the averaging domain that corresponds to the dissipation region and X is the initial random field. Because spatial averaging is performed in three dimensions, the dissipation regions in this study are rectangles (see Fig 2). The variance of the new random variable X_Vi is given by: (5) where is the variance of the random field X. According to earlier experience (e.g., [16, 17]), in this study, the field X is assumed to be a lognormal random field. However, both isotropic and anisotropic correlation structures are examined. As a correlation function within the random field, the Gaussian covariance function was selected and used in the numerical analyses, as shown in Eq (6). (6) where θ_x, θ_y and θ_z are fluctuation scales [7, 17]. To obtain the covariance function from Eq (5), the right-hand side of the equation must be multiplied by , i.e., . Note that in this study, . Eq (5) is applied to the three-dimensional scenario. Despite the assumed plane strain conditions, spatial averaging is performed in three dimensions. The coordinate system is shown in Fig 2. However, in further analyses, both horizontal fluctuation scales are assumed to be equal. Moreover, they are not distinguished and are denoted as θ_h (θ_h = θ_x = θ_y). The fluctuation scale along the z direction is called the vertical fluctuation scale θ_v (θ_v = θ_z). According to the failure mechanism shown in Fig 2, there are 16 regions in which spatial averaging is necessary. To determine these regions, a covariance matrix that describes the mutual correlations between each pair of regions is needed. To calculate the components of the covariance matrix, dedicated formulas must be derived. Derivation is possible due to the equations shown in Eq (7) (see [15, 18]).

(7)

The above equation determines the covariance between two random variables X_Vi and X_Vj resulting from the averaging (Eq (4)) corresponding to the dissipation regions V_i and V_j (here, rectangles are shown in Fig 2; e.g., P₁A₁P’₁A’₁). Note that i, j = 1, …, 16. As a result, the final size of the covariance matrix is 16×16. The covariance matrix is symmetric and positively definite. Below, some formulas for the covariance matrix components are derived. First, consider the dissipation region P₁A₁P’₁A’₁ (points P’₁ and A’₁ are not shown in Fig 2, but their coordinates are the same as those of P₁ and A₁, respectively, except that the y coordinates are shifted by α). Fig 2 shows the angle α (coloured red), which is used for parametric representation of the region P₁A₁P’₁A’₁. Indeed, the considered region can be expressed by using the following parametrization (Eq (8)):

Eq.d, by using the following parametrization the considered region can be expresses to be derived. is needed.

(8)

By substituting Eqs (8) and (6) into Eq (7), the following formula for the new variance of region P₁A₁P’₁A’₁ is obtained: (9)

Note that in Eq (9), the region P₁A₁P’₁A’₁ is denoted more simply as P₁A₁; this convention is also used later. Other variances of regions described by points P_i and A_i can be derived analogously.

In the case of region A₁A₂, the transformation shown in Eq (10) can be utilized as follows: (10) where the angle δ is shown in Fig 2 (coloured red).

(11)

Note that in Eqs (9) and (11), the coordinates along each axis are subtracted from each other. As a result, the coordinates of points P₁ and A₁ are reduced. However, the situation is different when the covariance between two different regions is calculated. To illustrate this situation, the formula for the covariance between regions P₁A₁ and A₁A₂ is given in Eq (12).

(12)

Based on an analogous procedure, the other covariance matrix elements were derived. To enable presentation of the covariance matrix, a simplification in the naming convention is introduced. First, start from the components considered above, i.e., and . Namely, the variances are denoted by w_i, where i is the index of the considered region (see Fig 2); the covariances are denoted by k_ij, where i ≠ j. As a result, the covariance matrix is given in Eq (13).

(13)

In the presented method, the covariance matrix size depends on the number of dissipation regions resulting from the failure mechanism. The size given in Eq (13) is valid for the failure mechanism of the 11th block shown in Fig 2. However, if a more rigid block is considered, the size of the covariance matrix is correspondingly larger.

The method for determining the covariance matrix for the considered failure mechanism is provided above. Next, based on C_X, the average set of undrained shear strengths is computed by the method described in Appendix A based on the approach proposed by Puła and Chwała [15]. Through the algorithm, the sixteen independent values of undrained shear strengths (c_u,1, …, C_u,₁₆), each dedicated to a specific dissipation region, are transformed to the correlated values (). The transformation is based on the Cholesky decomposition of the covariance matrix (Eq (13)). More details on this issue are also provided in the numerical algorithm section.

2.3. Simulated annealing optimization scheme

As mentioned earlier, the procedure for calculating bearing capacity described in section 2.1 is applicable for a known failure geometry. However, before the covariance matrix is determined, the optimal failure geometry must be found. Here, optimal means the failure geometry that provides the lowest possible bearing capacity for the specified soil strength parameters. Therefore, the optimization procedure is crucial to use the multiblock failure mechanism shown in Fig 2. As a basic optimization procedure, a simulated annealing scheme is utilized in this study [19, 20]. Moreover, earlier authors’ experiences with the application of simulated annealing to optimize the geometry of kinematic failure mechanisms noted that this approach could be used here [10, 21]. As a result, the approach described below is created. The flow chart of the procedure is shown in Fig 3; to facilitate navigation, the step numbers have been added in the below description. Generally, the optimization method is based on slight changes in the failure geometry in subsequent simulations. The method starts with the initial geometry parameters (step 1) from which the first bearing capacity is calculated (step 2). After that, the control parameter values are set (step 3). For each simulation, the corresponding bearing capacity is calculated using Eq (3). One bearing capacity is called the current value p_c, and this value is compared with those from the subsequent simulations p_n. This method allows for accepting a worse solution than the current one (p_c < p_n); however, the probability of such a scenario P_a decreases during the simulation process, and at the end of the simulation, this probability is nearly zero (only better solutions are accepted). Note that here, the better solution means the failure geometry that provides lower bearing capacity estimation. The concept of accepting worse solutions is due to overcoming local extrema; therefore, the final result is an approximation of the global extremum. In the simulated annealing scheme, there are some controlling parameters that control the simulation process (step 3). The formula used to calculate the acceptance probability is shown in Eq (14).

(14)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Flow chart for the optimization procedure.

A detailed description is in the text.

https://doi.org/10.1371/journal.pone.0231992.g003

The parameter T is responsible for decreasing the acceptance probability P_a during the simulation process. Therefore, T also decreases during simulation. Note that Eq (14) is used if p_c < p_n; however, if p_n < p_c, the formula given in Eq (14) reaches a value greater than one, and as a result, a better solution is always accepted (step 8). According to the literature (e.g., [22]), to effectively set the optimization method, the average initial value of P_a should be approximately 0.5. The crucial element of the described optimization procedure is a procedure for generating a slightly different failure geometry based on the current one. This is a so-called neighbouring set of geometric parameters (step 6). This procedure is used to generate each failure geometry considered by the simulation process (omitting the first geometry that is determined arbitrarily). The geometry of the failure mechanism shown in Fig 2 has 18 degrees of freedom. Therefore, those 18 parameters are the arguments of the objective function that is intended to minimize Eq (3). As a result, the arguments of the procedure for determining a neighbouring failure geometry are those 18 parameters. The following parameters are selected for this purpose: and coordinates; coordinates, i = 1, …, 10; and lengths l₂, l₃, l₄, l₁₀, l₁₁ and l₁₂. Each geometric parameter is changed in a random order by increasing its initial value by an amount generated from a uniform distribution U[−0.01 m, 0.01 m]. Moreover, the procedure takes into account additional limitations on the values of the changed parameters, for example, to prevent overlapping of adjacent blocks, namely, to ensure that, e.g., (see Fig 2). As mentioned earlier, for a new set of geometric parameters, the bearing capacity p_n is calculated using Eq (3) (step 7) and compared with the current value p_c (step 8). The acceptance probability of a worse solution is modified by the parameter T; however, some simulations are performed under constant P_a. In this study, for each P_a level, z = 2000 simulations are carried out (step 5). Then, P_a is decreased by modifying the value of T, namely, by multiplying T by an additional parameter: αT (step 9). The parameter α is kept constant during the simulation process and based on the literature and earlier experiences is assumed to be equal to 0.9. The simulation process continues as long as T > T_min (step 4). In this study, T_min = 10⁻⁸. When T reaches T_min, the algorithm ends, and the current bearing capacity p_c is treated as the optimal bearing capacity (step 10). Note that the values of the parameters that control the simulation process have to be selected individually for the considered issue.

According to the performed tests, the parameters chosen for the purpose of this study are definitely appropriate in terms of the achieved accuracy. The geometry of the failure mechanism obtained by using the optimization procedure described above for a non-random case was identical to the geometry obtained by Michałowski and Shi [3]. One example of this comparison is shown in Fig 4. For the initial geometry, the bearing capacity is approximately 210 kN/m. However, at the beginning stage of the simulation process, the value reaches over 280 kN/m. Despite this, during a simulation process with decreasing P_a, the current bearing capacity p_c moves towards the optimal value.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Value of the current bearing capacity p_c as a function of the number of simulations in the simulated annealing procedure.

The horizontal dashed line shows the final bearing capacity that is equal to that obtained by Michałowski and Shi [3]. The example was calculated for the following parameters to make the comparison possible: b = 1.0 m, γ = 20 kN/m³, c_u/γb = 1.0, t = 10 m and φ = 35°.

https://doi.org/10.1371/journal.pone.0231992.g004

The results obtained for the optimization procedure were also compared with the results obtained using finite element limit analysis implemented in OptumG2 [23] software for different numbers of finite elements and calculation methods. The obtained results are juxtaposed in Table 1. The results were obtained under plane strain conditions. To enable their comparison to the three-dimensional case, additional deterministic analyses were carried out by using ZSoil software [24]. In the case of t = 0.5 m bearing capacities of 142.5 kN/m, 140.6 kN/m and 136.8 kN/m were obtained for foundation lengths 6 m, 10 m and 16 m, respectively. In the case of t = 1.0 m bearing capacities of 215.6 kN/m, 198.7 kN/m and 191.2 kN/m were obtained for foundation lengths 6 m, 10 m and 16 m, respectively. As expected, the obtained bearing capacities are greater than those obtained for the plane strain condition and tends to decrease with an increase of foundation length.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Comparison of bearing capacity estimations for b = 1.0 m, γ 18 kN/m³, φ = 35° and two thicknesses of the sand layer: t = 0.5 m and t = 1.0 m.

https://doi.org/10.1371/journal.pone.0231992.t001

3. Numerical procedure

The numerical procedure for a one Monte Carlo simulation is provided below (see also Fig 5):

1. Step I introduces the values of the initial parameters. Set i = 1.
2. Step II starts Monte Carlo loop. If i ≤ N follow the steps below. If i > N go to step VII.
3. Step III generates 16 independent values of undrained shear strengths via a pseudo-random number generator ). All values are generated from the underlying normal distribution ( and ) for the initial stationary lognormal random field ( and ).
4. Step IV finds the optimal bearing capacity based on the soil parameters from step III transformed to the lognormal distribution by using the optimization procedure described in section 2.3.
5. Step V determines the correlated undrained shear strengths by transforming the independent values of undrained shear strength obtained in step III to correlated values using the covariance matrix (Eq (13)) determined for the optimal failure geometry from step IV. For this purpose, the algorithm shown in Appendix A is used. This algorithm is based on the Cholesky decomposition (e.g., [25]) of the covariance matrix C_X. Generally, the correlated parameters are obtained by calculating the product of a vector of uncorrelated parameters obtained in step III that was standardized and a triangular matrix that is the result of the Cholesky decomposition of C_X. Finally, a set of correlated values for the undrained shear strength is obtained (). These values are correlated by the covariance matrix C_X.
6. Step VI calculates the final bearing capacity by using the correlated undrained shear strengths () obtained in Step V. For this purpose, the optimization procedure from section 2.3 is used again. The resulting failure geometry and the corresponding bearing capacity are returned as the final results of the algorithm. Set i = i + 1, and go to step II.
7. Step VII provides the set of N bearing capacity estimates.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Flow chart of the numerical procedure.

A detailed description is provided in the text.

https://doi.org/10.1371/journal.pone.0231992.g005

In this study, for all performed analyses, the number of Monte Carlo simulations N is equal to 1000. The proposed approach is computationally efficient; namely, one three-dimensional simulation takes approximately 1 minute for one core of a standard notebook.

The numerical procedure detailed above allows consideration of three-dimensional soil spatial variability in the bearing capacity estimations for two-layered soil. Unfortunately, there are no published results that can be compared with the approach proposed in this study. However, this is possible in the case of plane strain conditions. In this particular case, the random finite limit analysis can be used for a comparison by using OptumG2 software. The comparison is made for and for two thicknesses of the top layer: t = 0.5 m and t = 1.0 m. The results for θ_v = 1.0 m plotted as a function of θ_h/θ_v are shown in Fig 6. The observed differences between both approaches are acceptable. Moreover, to some extent, the differences can be explained by the different correlation function assumed in the analyses performed by OptumG2 software (the Markovian one).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Comparison of the results obtained with the algorithm developed in this study and those obtained with OptumG2 software.

Two scenarios were considered for t = 0.5 m and t = 1.0 m. (a) mean values of bearing capacity; (b) standard deviations of bearing capacity.

https://doi.org/10.1371/journal.pone.0231992.g006

4. Considered scenarios

The main objective of this study is to assess the impact of a homogenous sand layer (top layer) on the overall random bearing capacity estimation in the case of a two-layered soil. A spatially variable, purely cohesive soil was assumed as the bottom layer. Moreover, the considered scenarios were selected to ensure creation of the failure mechanism shown in Fig 2. Therefore, relatively small thicknesses of the sand layer were examined. These thicknesses are suitable for working platform problems. The impact of vertical and horizontal fluctuation scales on random bearing capacity was investigated. Three vertical fluctuation scales were selected for numerical analyses, namely, θ_v = 0.25 m, θ_v = 0.50 m and θ_v = 1.00 m. These values were recently indicated by many authors, e.g., [26]. In the case of horizontal fluctuation scales, the values of θ_h = θ_v, θ_h = 5θ_v, θ_h = 10θ_v and θ_h = 20θ_v were examined. The two-dimensional failure mechanism was extended to three dimensions (Fig 2) to enable consideration of three-dimensional spatial averaging. As a result, the length of the averaging domain a was also examined. The following lengths were considered in the numerical analyses: a = 1.0 m, a = 3.0 m, a = 6.0 m and a = 10.0 m. The analyses that considered the impact of fluctuation scales were performed for two values of the undrained shear strength coefficient of variation: and . The mean value of undrained shear strength was assumed to be . Numerical analyses of all of the abovementioned parameter combinations were performed for two thicknesses of the top layer, i.e., t = 0.5 m and t = 1.0 m. The information for all scenarios that considered the impact of fluctuation scales is juxtaposed in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Combinations of fluctuation scales and averaging domain length a used in the numerical analyses.

Note that the shown scenarios are repeated for two thicknesses of the sand layer and two values of , as detailed in the text. Therefore, 192 scenarios were analysed.

https://doi.org/10.1371/journal.pone.0231992.t002

Moreover, for comparison purposes, scenarios with assumed infinite horizontal and vertical fluctuation scales were examined. By assuming θ_v = θ_h = ∞, a simple situation where the undrained shear strength is described by a single random variable is obtained. In such a case, there is no spatial averaging; therefore, the obtained standard deviations or coefficients of variation of bearing capacity are conservative estimates of the true values. These scenarios were intended to be analysed to separate the influence of dumping the spatial variability impact on the final bearing capacity coefficient of variation for the deterministic top layer of sand.

5. Results

Based on the Monte Carlo analyses, the mean values of the resulting bearing capacity μ_p and corresponding bearing capacity coefficients of variation were obtained. In Fig 7, the bearing capacity mean values obtained for t = 0.5 m are shown. The dashed lines represent , and the solid lines represent . Lower μ_p values are observed for greater values. This effect can be a result of using a lognormal probability distribution of c_u; a similar effect was described in [27]. In Fig 7a, 7b and 7c, the results for θ_v = 0.25 m, θ_v = 0.5 m and θ_v = 1.0 m are provided, respectively. Some interesting observations can be made based on Fig 7. First, μ_p, shown as a function of the horizontal fluctuation scale, reaches a minimum value for θ_v = 0.5 m and θ_v = 1.0 m. In the case of θ_v = 0.5 m, the minimum value in most of the considered cases is observed for θ_v/θ_h = 10. However, for θ_v = 1.0 m, the minimum value is observed for θ_v/θ_h = 5. By multiplying the vertical fluctuation scale by the ratio for which the minimum is observed in both cases, the value of the horizontal fluctuation scale, θ_v = 5.0 m, is obtained. This observation indicates that the observed worst-case correlation length (e.g., [28]) for the mean value is independent of the vertical fluctuation scale. Note that Fig 7a does not refute this observation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Bearing capacity mean values μ_p for t = 0.5 m as a function of the horizontal fluctuation scale and length of the averaging domain a.

Three vertical fluctuation scales are considered. The dashed lines represent , and the solid lines represent .

https://doi.org/10.1371/journal.pone.0231992.g007

Fig 7 clearly shows the differences in the bearing capacity mean values for different values of a. In each case, a longer averaging domain provides a higher μ_p. This can also be explained by the abovementioned effect of obtaining a lower μ_p for a greater coefficient of variation c_u. Thus, a longer averaging domain a provides greater averaging; therefore, the overall coefficient of variation c_u decreases as the variances defined in Eq (5) are more significantly reduced. It is worth mentioning that all obtained bearing capacity mean values are below the values of the deterministic scenario that is denoted by the horizontal dashed line. Note that the relative differences in the obtained μ_p values are not very significant (see values on the vertical axis). Analogous observations can be made for t = 1.0 m (Fig 8). However, the relative differences in the obtained μ_p values are smaller than those for t = 0.5 m. This demonstrates the effect of damping the impact of soil spatial variability on bearing capacity with increasing thickness of the homogenous layer t. This effect is described later in more detail.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Bearing capacity mean values μ_p for t = 1.0 m as a function of the horizontal fluctuation scale and length of the averaging domain a.

Three vertical fluctuation scales are considered.

https://doi.org/10.1371/journal.pone.0231992.g008

Very significant differences among the considered scenarios are observed in the coefficient of variation of bearing capacity v_p. The results for are shown in Fig 9. The horizontal dashed and dotted lines shown in Fig 6 indicate v_p for special scenarios. The dashed line indicates the coefficient of variation of bearing capacity for t = 0 and infinite values of the vertical and horizontal fluctuation scales. In this case, the considered issue is degenerated into a one-layered purely cohesive soil with no spatial variability and . Therefore, the resulting coefficient of variation of bearing capacity (as shown in Fig 9). However, the dotted lines in Fig 9a and 9b indicate the resulting coefficient of variation of bearing capacity for the value of t specified in each figure and infinite values of the vertical and horizontal fluctuation scales. For t = 0.5 m, the corresponding v_p,t = 0.211; thus, the coefficient of variation of bearing capacity decreases in comparison with the scenario in which t = 0 m. A further decrease is observed for t = 1.0 m, for which v_p,t 0.157. Thus, it is easy to understand why the v_p values observed for scenarios with finite fluctuation scales (see Table 1) are lower than v_p,t. As shown in Fig 9, a shorter vertical fluctuation scale provides a smaller value of v_p. The dotted lines represent the limit values of v_p when the fluctuation scales reach infinity. The information is intentionally presented in the form shown in Fig 9 to separate the impact of the homogenous sand layer from the impact of fluctuation scales on the resulting bearing capacity coefficients of variation. As expected, for a longer averaging domain a, a lower v_p is obtained. The observed differences between the extreme values, i.e., a = 1 m and a = 10 m, are significant, and they indicate the necessity of considering the real dimension of the structure if soil spatial variability is taken into account. In such situations, the two-dimensional simplification of three-dimensional issues provides more conservative bearing capacity estimates. This conclusion has recently been noted [1,21,29].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Bearing capacity coefficient of variation v_p for as a function of the horizontal fluctuation scale and length of the averaging domain a.

Three vertical fluctuation scales and two thicknesses of the sand layer, i.e., t = 0.5 m (a) and t = 1.0 m (b), are considered. Detailed descriptions are provided in the text.

https://doi.org/10.1371/journal.pone.0231992.g009

Results analogous to those shown in Fig 9 for are shown in Fig 10. The horizontal dashed and dotted lines have the same meaning as those in Fig 9.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Bearing capacity coefficient of variation v_p for as a function of the horizontal fluctuation scale and length of the averaging domain a.

Three vertical fluctuation scales and two thicknesses of the sand layer, i.e., t = 0.5 m and t = 1.0 m, are considered. Detailed descriptions are provided in the text.

https://doi.org/10.1371/journal.pone.0231992.g010

Both Figs 9 and 10 indicate that calculating v_p,t can be a very easy method to conservatively estimate the random bearing capacity in the case of a lack of information about fluctuation scales. The reduction in the bearing capacity coefficient of variation is a result of considering only the homogenous soil layer. Therefore, it can easily be estimated. To illustrate that possibility, Figs 11 and 12 show the results obtained by analysing scenarios with θ_v = θ_h = ∞. First, in Fig 11, the obtained bearing capacity mean values are plotted as a function of the homogenous sand layer thickness t. Fig 11 presents the results for and . For each undrained shear strength coefficient of variation, three values of internal angle are considered, i.e., φ = 30°, φ = 35° and φ = 40°. Generally, the bearing capacity mean values increase with an increase in the sand layer thickness. Moreover, greater μ_p values are observed for greater angles of internal friction that characterize the sand layer. The differences between the scenarios with and are practically negligible.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Bearing capacity mean values μ_p for and as a function of the sand layer thickness t.

Three values of the angle of internal friction are considered.

https://doi.org/10.1371/journal.pone.0231992.g011

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Bearing capacity coefficient of variation v_p for and as a function of the sand layer thickness t.

Three values of the angle of internal friction are considered.

https://doi.org/10.1371/journal.pone.0231992.g012

Fig 12 shows the bearing capacity coefficients of variation that correspond to the mean values shown in Fig 11. The obtained results indicate that v_p decreases with an increase in the sand layer thickness. Moreover, the observed differences in v_p between different friction angles are relatively narrow. A simplified method for which the results are shown in Figs 11 and 12 can be used to estimate the reduction factor in the bearing capacity coefficient of variation that results from the thickness of the homogenous sand layer.

6. Summary and conclusions

An original adaptation of the two-dimensional failure mechanism for two-layered soil proposed by Michałowski and Shi [3] to random analyses is shown in the present study. The adaptation was necessary for numerical analyses that are based on the upper bound approach and are dedicated to the scenario of a spatially variable bottom layer. A method of Vanmarcke’s spatial averaging in conjunction with a kinematic failure mechanism [15] was used in this study. The use of this method requires its adaptation to the considered scenario. As a result, an original approach for random analysis of the bearing capacity of two-layered soil and three-dimensional spatial averaging was proposed in this study. The created algorithm is capable of considering all important parameters, such as geometrical parameters, thickness of the top layer, its unit density and fluctuation scales, that describe the variability of undrained shear strength along the vertical and horizontal directions. The proposed approach is very efficient (one three-dimensional simulation takes approximately 1 minute for one core of a standard notebook). The method is dedicated to a situation in which the top layer is homogenous sand (generally cohesionless soil) and the bottom layer is purely cohesive soil that is spatially variable. This arrangement of the soil layers has practical meaning because it reflects the conditions in which working platforms operate, where the top layer is man-made, and the assumption that it is homogenous is justified. The proposed approach is in accordance with the upper bound theorem, and thus, the obtained bearing capacity must be treated as the upper limit of the true limit load. The proposed approach is limited to bearing capacity estimations and is not dedicated to serviceability limit states, e.g., Johari et al. [30]. In the performed analyses, seismic loadings are not examined, and thus, the obtained results are reliable for a static scenario (this can influence the obtained results, e.g., [31]). However, this type of analyses will be explored in future research.

Based on the results presented in section 5 and the description of the proposed algorithm given in section 2 and section 3, the following conclusions can be drawn:

1. The algorithm proposed in this study demonstrates that the method of joining Vanmarcke’s spatial averaging with kinematic failure mechanisms can deliver an efficient approach for analysing two-layered soil problems if the bottom layer is spatially variable. Further efficiency improvement is possible if a constant covariance matrix is applied [32].
2. The obtained results indicate that the observed reduction in the bearing capacity coefficient of variation (see Figs 6 and 7) is caused mostly by three factors: thickness of the homogenous sand layer, length of the averaging domain a, and fluctuation scales that characterize the spatial variability of the undrained shear strength within the bottom layer. The impact of surcharge load is not considered in the performed analyses. The soft layer is assumed to be weightless (see arguments given in section 2.1); however, in future studies, it would not be very difficult to extend the analyses by soil weight consideration.
3. Three-dimensional issues are of special importance when soil spatial variability is taken into account in the numerical analyses. Two-dimensional simplifications provide a more conservative approach that can be significant in some scenarios, as shown in Figs 6 and 7.
4. As shown in the results section, the impact of dumping the feature of the thickness of the homogenous sand layer can be separated from the other impacts detailed in the second conclusion. An example of this separation is shown in Fig 9, where the reduction in the bearing capacity coefficient of variation is caused only by the homogenous sand layer. As expected, the reduction is greater for a greater sand thickness. This simple approach can be utilized in practical applications for creating graphs that allow nearly instant determination of the reduction in the bearing capacity coefficient of variation.

Appendix A

The algorithm for computing average undrained shear strengths in the case of 11 rigid block failure mechanisms.

1. Step 1. Find the mean value and the variance of the underlying normal distribution based on and from step I (see section 3). (A.1) (A.2)
2. Step 2. Take 16 values of undrained shear strength generated in step I in section 3 as a vector of stochastically independent components .
3. Step 3. Transform vector to the lognormal distribution , where i = 1, …, 16.
4. Step 4. Determine the covariance matrix [C_X] (see Eq (13)) based on the optimal failure geometry established in step II in section 3.
5. Step 5. Compute the correlation coefficients for the covariance matrix determined in step 4. (A.3)
   where i, j = 1, …, 16.
6. Step. 6. Transform [r_X] to the corresponding normal underlying distribution correlation matrix [r_Y], (see [7]): (A.4)
   where (A.5)
7. Step 7. Transform [r_Y] to the corresponding normal underlying distribution covariance matrix [C_Y], (A.6)
   where (A.7)
   where i, j = 1, …, 16.
8. Step. 8. Calculate the Cholesky decomposition of [C_Y]: (A.8)
   where [L] is the lower triangular matrix. This matrix has real and positive diagonal entries (see [25]).
9. Step. 9. Standardize vector obtained in step I in section 3. (A.9)
   where i, j = 1, …, 16.
10. Step 10. Compute vector of correlated entries by using the following theorem (see [17]): if [C_Y] is a positively definite n × n matrix, is an n-dimensional vector whose components are independent Gaussian standard random variables, and [L] is a lower triangular matrix such that Eq (A.8) is true, then, the random vector defined as is the Gaussian vector with the covariance matrix [C_Y]. As a result of using this theorem, entries of the vector have expected values equal to zero and the covariance matrix [C_Y]. To shift the components of to proper mean values, they are added to , where i = 1, …, 16.
11. Step 11. Transform to lognormal distribution: (A.10)
    where i = 1, …, 16.
    The obtained components of are used in step IV in the algorithm given in section 3.