Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil using random finite element limit analysis

The assumed soil is of c - φ type. Both soil strength parameters were described with random fields. It is assumed that these fields are stationary and anisotropic. Since c and φ are usually reported as negatively correlated and such correlation typically results in a lower standard deviation of bearing capacity (cf. Cho & Park 2010) and consequently lower probability of failure, the fields were also assumed uncorrelated for conservatism. Both probability distributions were assumed to be lognormal. To generate the random fields, the FSM was utilised. This approach approximates the autocorrelation function by applying a finite Fourier series. The value of the field at a given point in 2D space with coordinates x and z is calculated as the following double sum: (1) WX=Wx,z=∑−mm∑−nnamn+ibmnexpi2πm⋅xLx+n⋅zLz W\left( {\bf{X}} \right) = W\left( {x,z} \right) = \sum\limits_{ - m}^m {\sum\limits_{ - n}^n {\left( {{a_{mn}} + i{b_{mn}}} \right)\exp \left( {i2\pi \left( {{{m \cdot x} \over {{L_x}}} + {{n \cdot z} \over {{L_z}}}} \right)} \right)} } where a_mn, b_mn are independent random variables with a normal distribution N(0,σ_mn) with σ_mn dependent on the chosen autocorrelation function and L_x and L_z are the dimensions of the domain over which the simulations are performed. The length of the series is determined by the numbers m and n which depend on the assumed precision and the assumed autocorrelation function. Please note that the description of the random field by a continuous function, that allows the value of the field realisation at any point to be determined independently of the element mesh, is a mandatory condition for RFELA analysis which implements the mesh adaptation procedure. The results obtained by the FELA method are strongly affected by the FE mesh fitting to the failure mechanism obtained in a given realisation. The implementation of an iterative procedure that allows the mesh shape to be adapted to the mechanism makes it possible to increase its accuracy without an excessive increase in the number of elements. A description of the mesh adaptation algorithm is provided in the later sections.

The FSM algorithm generates random fields with normal distribution. Since the distribution of soil parameters was assumed as lognormal, the fields with the parameters for the underlying normal distribution RF_Y(X) were first generated. The mean value (μ_g) and standard deviation (σ_g) of the underlying normal distribution can be calculated using the analytical formulas that link these moments to moments of the target lognormal distribution {μ,σ}, namely, (2) σg2=ln1+σ2μ2, μg=lnμ−12σg2 \sigma _g^2 = \ln \left( {1 + {{{\sigma ^2}} \over {{\mu ^2}}}} \right),\;\;\;{\mu _g} = \ln \left( \mu \right) - {1 \over 2}\sigma _g^2 After generation, the values of the underlying normal field were further transformed to lognormal field RF_X(X) according to formula (2). (3) RFXX=expRFYX, R{F_X}\left( {\rm{X}} \right) = \exp \left( {R{F_Y}\left( {\rm{X}} \right)} \right), where X={x,z} denotes the coordinate of the point that belongs to the random field.

Both random fields were assumed to be characterised by a Gaussian autocorrelation function according to the following formula: (4) ρτx,τz=exp−πτxθx2+τzθz2 \rho \left( {{\tau _x},{\tau _z}} \right) = \exp \left\{ { - \pi \left[ {{{\left( {{{{\tau _x}} \over {{\theta _x}}}} \right)}^2} + {{\left( {{{{\tau _z}} \over {{\theta _z}}}} \right)}^2}} \right]} \right\} where τ_x, τ_z are, respectively, the distance between two points on the horizontal (x) and vertical (z) axes, and θ_x, θ_z are the horizontal and vertical scales of fluctuation. The Gaussian autocorrelation function was chosen because of the low length of the associated Fourier series (compared with, e.g., exponential function) which affects the time of random filed values calculation.

As mentioned, the RFELA method utilises FELA to solve individual realisations of the given problem. The latter is one of the possible numerical methods capable of solving the limit state problem, e.g., of the load-bearing capacity of a shallow foundation. The method requires several assumptions to be made which significantly simplify the model while omitting the knowledge of the deformation state. On the other hand, it enables a considerable reduction in calculation time compared with, e.g., the FE or the finite difference methods. Most importantly, FELA provides rigorous lower and upper limits of bearing capacity.

The LB and UB solutions are determined, respectively, by searching for the maximum load corresponding to the statically admissible stress field and the minimum load for which the kinematically admissible displacement velocity field is developed. Both issues represent a mathematical optimisation problem. This article uses the approaches described by Lyamin and Sloan (2002a) in the context of LB estimation and Makrodimopoulos and Martin (2008) regarding UB estimation. A perfectly plastic material model with an associated flow rule is assumed. The FE mesh consists of triangular elements for which a linear shape function is used for the stresses and a quadratic shape function is used for the displacement velocities when estimating the LB and UB, respectively. As previously mentioned, the Coulomb–Mohr material model with strength parameters c and ϕ is assumed. This model can easily be written in the form of a conic quadratic constraint, which enables finding the solution of both UB and LB using second-order cone programming. The letter appears to be a very efficient method for solving optimisation problems.

The LB and UB problems are weakly dual, which means that their solutions differ by a certain gap due to discretisation error. On the other hand, strongly dual problems with the same solutions can be formulated to each of the task (Ciria Suárez, 2004). The dual problem to the kinematic approach involves searching for an admissible stress distribution (slightly different from the one in static approach), while the dual problem to the static approach involves searching for an admissible failure mechanism (slightly different from the one in kinematic approach). In this article, both problems are formulated by searching for the stress distribution (i.e., the original static problem and the dual to kinematic problem are solved), because, as noted by Krabbenhoft et al. (2005), a dual problem to the kinematic approach could be easier to solve.

Both mentioned problems can be described by the following equations provided by Podlich (2018): (5) maximise:αsubject to:BTσ=αp+p0fσe≤K0 ∀ e=1,2,…,ne \matrix{{{\rm{maximise:}}} \cr \alpha \cr {{\rm{subject}}\;{\rm{to:}}} \cr {{B^T}\sigma = \alpha p + {p_0}} \cr {f\left( {{\sigma ^e}} \right){ \le _K}0\;\forall \;e = 1,2, \ldots ,{n_e}} \cr } where α is the searched load value multiplier and σ is the global stress vector different for LB and UB analysis. Additionally, B^T contains the coefficients of the static equilibrium equations significantly different for the LB and UB limit analysis. Vectors p and p₀ contain the optimised and permanent load. In both cases, the Coulomb–Mohr plasticity condition is defined for all elements in the form of SOCP (Second Order Cone Programming).

On the other hand, in both approaches, the mesh adaptation algorithm utilises parameters of plastic multipliers, which are obtained by solving duals to the mentioned problems. To solve the problems, the commercial software MOSEK run from the Matlab environment was used. According to the knowledge of the authors of this study, this is the most effective software currently available for solving many different types of optimisation problems (including SOC) (Podlich et al., 2014). The software uses the primal–dual method, which allows the automatic solution of a primal and dual problem, providing a full range of information on stresses, plastic multipliers and displacement velocities for both LB and UB estimates.

The degree to which the mesh is adapted to the shape of the failure mechanism has a significant impact on the results obtained from the optimisation problems. Obviously, increasing the mesh density over the entire domain will greatly reduce the difference between the LB and UB estimates, however, the number of nodes resulting from such mesh is computationally prohibited. For more efficient computations, elements should be densified in a specific manner only in the zone of interest.

For mesh generation, a tool called MESH2D (Engwirda, 2014) for the Matlab environment, which allows the generation of triangular meshes suitable for use in the FELA method, is utilised. This generator is based on the Delaunay method. It includes options to modify the meshes both by subdividing the indicated elements and requesting a density function in the form of a dimension expectations map. The mesh adaptation procedure utilised in the calculations described below assumes that the yielding elements should be subdivided. These elements are indicated based on the plasticity multipliers obtained from both the UB and LB analyses.

In the first step, the same low-density mesh is generated regardless of the realisation of the random field. Smaller elements occur only at the corners of the foundation due to the use of so-called ‘element fans’, which are necessary to reduce the influence of the singularity located at the specified points (Lyamin & Sloan, 2003). Performing a FELA analysis with such a generated mesh is extremely time-efficient and gives a rough indication of the potential plastic zones. In the second step, a mesh density function, which assumes minimum dimensions for the elements that have yielded and maximum dimensions for all the others, is applied. This operation enables the refinement of the areas where the plasticity multiplier is present without an excessive increase in the number of elements. In the next steps, the elements with the highest plastic multiplier values are selected and subdivided using the MESH2D tool. The limit values of the plastic multiplier are chosen depending on the individual task. The algorithm must try to match the shape of the mesh with the results of both estimations. As mentioned, often the static and kinematic mechanisms may differ. The algorithm applied in the current work consists of four steps. A graphical representation of the discussed algorithm for a zoomed area around the foundation is illustrated in Figure 1.

Figure 1:

A graphical representation of the mesh adaptive procedure.

The study investigates the effect of the foundation depth and the value of vertical and horizontal scales of fluctuation on the bearing capacity of a shallow foundation placed on spatially variable soil. The first analysis results for a slightly different formulation of the FELA algorithm and a different mesh adaptive procedure were presented by Puła et al. (2022). The deliberation was limited to the shallow foundation with D=0.0 m only. A full analysis was carried out in the present work.

The foundation was assumed as a rigid body with a load uniformly distributed over the top face. It was allowed to rotate if such movement resulted from an optimal failure mechanism. Boundary conditions were adopted as indicated in Figure 2. The width and height of the model were assumed to be 20.0 m and 10.0 m, respectively, which allowed to ensure that no failure mechanism reached the area boundary. The width of the foundation was 1.0 m. Three foundation depth scenarios with D=0.0 (foundation placed on the ground surface), 0.5 and 1.0 m were considered. The surface model is a frequently adopted theoretical model. On the other hand, the typical direct foundation depth oscillates around 1.0 m. A value of half a metre was assumed as an intermediate value to present a potentially nonlinear influence of embedded depth. It was assumed that the foundation established at a depth greater than zero is a rigid body that extended above the ground surface (no backfill soil overlying the foundation was assumed). Please note that the foundation material was assumed as infinitely strong (conditions of the Mohr–Coulomb criterion were not checked within the foundation).

Figure 2:

Model geometry with boundary condition.

For each foundation level, a set of five horizontal fluctuation scales θ_x={1.0,2.0,4.0,8.0,16.0 m} and a set of three vertical scales θ_z={0.25,0.5,1.0 m} were considered. The analyses included pairs formed from all possible combinations of these values. The values of cohesion (μ_c,σ_c) and internal friction angle (μ_φ,σ_φ) were described by the random fields realisation generated for these values as described earlier. The parameters of the lognormal distribution assumed for generation for both fields were assumed as μ_c=29.0kPa, σ_c=7.0 kPa for cohesion and μ_φ=12.41°, σ_φ=1.15° for friction angle. The soil weight was assumed to be 20 kN/m³. The specified parameters were adapted after Zaskórski and Puła (2016), who utilised them in their investigations of clay soil layers. Probabilistic analysis was carried out using the Monte-Carlo method. For each problem (each assumed pair of scales of fluctuation), 10,000 MCS was solved. The LB and UB estimates of bearing capacity were obtained for each problem realisation. A summary of the results obtained for all considered vertical and horizontal scales of fluctuations and foundation levels is reported in the following section.

The results described in the previous chapter are based on 10,000 MCS realisations. By using RFELA for a single fluctuation scale, these amount of calculations were carried out in about 8–10 hours. It is worth noting that, although the lower estimation results always estimate the load conservatively, even such a large number of realisations do not allow a sufficiently accurate estimation of the probability of failure for the ultimate limit state of the foundation. According to EN-1990, for the standard 50-year reference period, the allowable reliability index should be greater than β=3.8 which corresponds to the failure probability lower than p_f=7.23∙10⁻⁵. According to Kolmogorov’s inequality, 4,194,304 realisations are required to obtain a credible (with 95% confidence and a relative error of less than 10%) Monte-Carlo estimator of such a low probability (Kawa, 2023). Such a number of simulations is not feasible in a reasonable time. This number could be probably significantly reduced by using one of the variance reduction techniques, e.g., Subset Simulation (Au & Beck, 2001), but using such a method is out of the scope of the current work.

As mentioned above, as the consequence of the rigorousness of FELA, the results described above are consistently characterised by a lower mean value and lower standard deviation for LB bearing capacity estimations and greater mean and greater standard deviation for UB bearing capacity estimations. While this relationship is less intuitive in the case of standard deviation, it was clearly observed in all considered cases. It is important to note that the lower mean value of bearing capacity distribution with fixed standard deviation always lead to greater probability of ultimate limit state failure and thus conservative estimation of allowable load. On the other hand, an increase in the standard deviation of structural response results in an increased probability of failure, which in the consequence also lead to conservative allowable load estimation. For this reason, adopting a distribution with a lower mean based on the results of the LB (static) approach and a higher standard deviation based on the UB (kinematic) approach appears to be a safe and conservative approximation of the probabilistic bearing capacity behaviour. Such a distribution with the following moments will be referred to as the mixed solution. (6) μ=minμLB,μUB=μLBσ=maxσLB,σUB=σUB \matrix{{\mu = \min \left\{ {{\mu _{LB}},{\mu _{UB}}} \right\} = {\mu _{LB}}} \cr {\sigma = \max \left\{ {{\sigma _{LB}},{\sigma _{UB}}} \right\} = {\sigma _{UB}}} \cr }

Based on the distributions with values of mean and standard deviation presented above including upper-, lower- and mixed-approach distributions, the values of the reliability-based design loads for the analysed shallow foundations were estimated. By following the mentioned EN-1990 guidelines, p_f=7.23∙10⁻⁵ (β=3.8) was assumed as an allowable probability condition. For a given distribution, the allowable reliability-based design load value can be obtained as: (7) qd=Φ−1pf {q_d} = {\Phi ^{ - 1}}\left( {{p_f}} \right) where Φ⁻¹ denotes the inverse cumulative distribution function of the considered distribution

Figure 11 illustrates the results in the form of diagrams of the allowable load q_d and the safety factor, understood as the ratio of the mean value to the value of the allowable load. Both variables were plotted as a function of the horizontal scale of fluctuation. Subsequent vertical scales of fluctuation have been marked with solid, dash-dotted and dashed lines, respectively. Once again, results for different foundation levels are presented in different columns. The red colour represents the upper estimate, the blue the lower estimate and the green the mixed approach. As the horizontal scale of fluctuation increases, a decrease in the allowable load value of the foundation and an increase in the safety factor can be observed. The safety factor determines the reduction in the ultimate bearing capacity that must occur for the assumption of an acceptable probability of structural failure to be satisfied. For the same values of the fluctuation scales, increasing the foundation level from D=0.0 to 1.0 m means a decrease in the factor by 10%–25%. Again, the general trends are similar to the values obtained in other studies (Pieczyńska-Kozłowska et al., 2015;, Kawa & Puła, 2020). The results from Figure 11 are also summarised in Tables 4, 5 and 6.

Figure 11:

Diagrams of allowable loads and safety factors as a function of the horizontal scales of fluctuations θ_x with a constant vertical scale of fluctuation for all considered foundation levels.

Please note that the approach utilising the distribution estimated based on the LB usually provides a lower (more conservative) estimate of the allowable load than the analogous approach based on the UB. However, since the rate of convergence differs between the LB, UB, mean and standard deviation, it may be possible that the LB results will estimate the mean value of the real solution better than the standard deviation. If such a standard deviation is too low, the allowable load distribution based on the LB may not be a conservative one. It should be noted that for the different scenarios, the mixed approach always provides a conservative and quite precise (as indicated by the small differences between the values provided by the different approaches) estimation. Moreover, this precision could be further improved (if the current precision level in the given problem is not satisfactory) by increasing the number of steps in the adaptive mesh refinement procedure. Since the moments of the upper and lower approaches can be accurately estimated even for lower numbers of realisations, the mixed approach provides a rigorous, conservative and effective approach for reliability analysis. With these features, the proposed approach seems to be applicable in engineering practice.

As has been mentioned in the introduction, the application of the FELA approach to the estimation of the bearing capacity of a shallow foundation has a recognisable history. However, until now, the rigorousness of this technique has not been fully exploited. The main motivation for writing this paper is to show a relatively simple application of this rigorousness for the conservative estimation of probability-based design values. As has been shown in the above section, such an application is possible.

It should be noted that the current study is subject to certain limitations. One of the most obvious is considering the structure in 2D as a plain strain problem. As a consequence, the value of the out-of-plane fluctuation scale was assumed to be infinite. The issue of modelling a problem of spatially variable soil, which is, in general, a 3D problem, in 2D space is well described in the literature for different issues (e.g., Liu et al. 2022). As shown in the paper (Chwała & Puła, 2020), omitting the fluctuation in out-of-plane direction can strongly influence the bearing capacity result by increasing the standard deviation of the structure response. This means that 2D analysis of strip foundation based on random field lead to crude but conservatism estimation of bearing capacity probability distribution. A full 3D analysis of spatially variable soil was beyond the scope of the current work. However, we believe that the proposed method for estimating failure probability or allowable load values based on the mixed approach remains applicable to 3D analyses.