Bearing capacity of eccentrically loaded strip footing on spatially variable cohesive soil

In recent decades, following the widespread growth of interest in risk assessment and management, the concept of reliability-based design has developed in civil engineering. The idea of designing a structure in such a way that the calculated probability of failure does not exceed certain specified values (depending on the function of the object, the type of failure or its consequences) has spread and, although rarely used in engineering practice, is nowadays part of European and International standards for construction design (EN 1990, 2002, ISO 2394, 2015). The approach associated with this trend has been developing for several decades also in geotechnics, where spatial variability of soil parameters is considered the most important source of uncertainty (e.g., Cami et al., 2020; Pieczyńska-Kozłowska et al., 2021). It is worth mentioning that variability of soil parameters has often a much greater effect on the structure reliability than variability of other structure or material parameters which are usually (with some exceptions, e.g., Wyjadłowski et al., 2018) easier to control. Furthermore, the problem of spatial variability applies to all structures founded on the soil (i.e., the vast majority of the total number of structures).

Several works in the last two decades have considered modelling the spatial variability of soil parameters using so-called random fields (Vanmarcke, 2010). The application of random fields in geotechnical design is an approach explicitly recommended by Annex D of the latest edition of ISO 2394 (2015). In such an approach, the analysis is performed within the Monte Carlo framework by solving the given problem an appropriate number of times, each time with different ‘realisation’ of random fields (modelling the properties considered) and then analysing the probabilistic characteristics of the results obtained. The numerical variant consists of generating realisations of random fields discretised over a finite element (FE) or a finite difference (FD) mesh, so that in each mesh element, the parameters are constant and usually represent an average of the field values over the area of the element. Despite the facts that the concept of this random finite element method (RFEM) in application to foundation bearing capacity was formulated two decades ago (Griffiths and Fenton, 2001; Fenton and Griffiths, 2003) and that the random fields can also be used in some other semi-analytical approaches (cf. Chwała, 2020a, 2020b), the mentioned numerical method remains very popular (probably due to its universality). So far (together with the analogous random finite difference method), it has been used to analyse a number of geotechnical problems, in particular the random bearing capacity of strip footings (e.g., Fenton and Griffiths, 2003; Vessia et al., 2009; Pieczyńska-Kozłowska et al., 2015), rectangular footings (e.g., Kawa and Puła, 2020; Li et al., 2021), slopes (e.g., Griffiths and Fenton, 2004; Jha and Ching, 2013a; Huang et al., 2019) and different types of retaining walls (e.g., Fenton and Griffiths, 2005; Sert et al., 2018; Kawa et al., 2019; Kawa et al., 2021).

As indicated in the work by Kawa and Puła (2020), little attention has been paid so far to probabilistic modelling of eccentrically loaded strip or rectangular footings. This problem was considered, for example, in the work by Soubra (2009). However, the influence of soil spatial variability was neglected (in the cited work, soil parameters were assumed as random variables, which means that they are constant over the considered space). To the best of authors’ knowledge, the only work that considers the spatial variability of the medium parameters in the analysis of the bearing capacity of an eccentrically loaded strip footing is the conference paper by Ali et al. (2016). In the cited study, the boundary problem is solved using the random finite element limit analysis (RFELA) technique. The effect of the length of the fluctuation scale (correlation radius) on the failure surface is analysed. This surface is presented in the vertical force–bending moment space (the two quantities were assumed as loads for the strip footing). An important simplification adopted in the paper is that the random cohesion field considered is assumed to be isotropic (the fluctuation scales in the horizontal and vertical directions are identical), which is not consistent with the results of other studies (e.g., Lloret-Cabot et al., 2014) and may lead to neoconservative results (cf. Vessia et al., 2009; Pieczyńska-Kozłowska et al., 2015).

The present paper complements and expands the study mentioned above. Similar to the work of Ali (2016), it was assumed that the soil medium is purely cohesive and that cohesion is the only parameter modelled by a random field. However, in the present work, the field of cohesion was assumed to be an anisotropic random field with the horizontal fluctuation scale larger than the vertical one. For easier scaling of the problem solution, it was assumed that the soil is weightless. The analysis was carried out using a typical RFEM; the individual random field realisations were generated using the Fourier series method (FSM; Jha and Ching, 2013b) and the considered boundary value problems were solved in ABAQUS software. For eight combinations of typical values of the fluctuation scale (representing both weak and strong anisotropy), the effect of the eccentricity of the vertical foundation load on the mean, standard deviation (SD) and coefficient of variation (CoV) of the bearing capacity of the soil was investigated. In the final part of the paper, an example reliability analysis was also performed showing the influence of eccentricity on the reliability index.

This paper is organised as follows. The next section briefly describes the basic assumptions for the study. In the following section, the numerical model, the preliminary deterministic tests and the probabilistic analysis performed are described. The results of the analysis are then presented. A brief section of conclusions closes the article.

As mentioned in the ‘Introduction’ section, the RFEM in geotechnical applications involves modelling soil parameters using a random field discretised over a FE mesh. The phrase ‘random field’ refers to a spatial generalisation of a stochastic process (see, e.g., Vanmarcke, 2010) and is associated with a whole group of different field classes. The class considered in this paper is restricted to ergodic and weakly stationary random fields, which is a typical assumption in RFEM applications. For this class of fields (Doob, 1953), the probability distribution does not depend on the position of a point, where it is tested, and similarly, the autocorrelation function defining the correlation between values of field in two different points does not depend on the position of these points, but only on their relative positions, that is, the direction and length of the vector between them (e.g., Kawa et al., 2021). Thus, the point probability distribution and the stationary (relative position dependent) autocorrelation function with specified parameters, that is, the so-called scales of fluctuations (SOFs), form a complete set of information needed to define such a field. A satisfactory description of probability distribution can usually be obtained by choosing the distribution type (the normal or log-normal distributions are most commonly accepted) and assuming or determining its first two moments, that is, mean and SD. It is much more difficult to credibly estimate the values of SOF. Identification of these parameters, which can be understood as the average dimensions of a cluster within which the field values are strongly correlated, is a problem widely described in the literature (Lloret-Cabot et al., 2006; Cami et al., 2020; Bagińska et al., 2018; Ching et al., 2018). It is generally known that the value of SOF in the horizontal direction is larger than the value of SOF in the vertical direction. Due to the way in which soil investigations are carried out, horizontal SOF is considered more difficult to identify and is often the subject of parametric studies.

If bearing capacity of footings is considered, only the strength parameters of the soil medium are usually modelled by random fields (cf. Puła and Zaskórski, 2015). In the present paper, it is assumed that the soil under consideration is a purely cohesive medium, and thus, the only parameter modelled by the random field is its cohesion c. Furthermore, it is assumed that the soil considered is weightless. Although such a model (due to the lack of both internal friction and unit weight of the soil) may seem simplistic, it should well describe the bearing capacity of the strip foundation in undrained conditions. Furthermore, both of these assumptions enable the use of a well-known limit state solution for the validation of the problem and allow easy scaling of the obtained numerical solution for other values of the width of foundation B (or rather B* as will be explained in a further section) than assumed in the calculations: for such a model, the bearing capacity (in terms of stress) does not depend on the width of the footing and SOFs, but only on the width/SOFs ratio.

It was assumed that the random field of cohesion considered is described by a normal probability distribution with mean μ_c = 20 kPa and SD σ_c = 2 kPa. For positive mechanical parameters (such as cohesion), the normal distribution, due to its ability to take negative values, can only be assumed for relatively small values of the CoV of the modelled quantity δ_c = σ_c/μ_c (cf. Kawa et al., 2019), say not greater than 10%–20%. The value of 10% was assumed here for this reason. This value, although mentioned in the literature, is associated with relatively low variability of cohesion in the soil (according to the work by Phoon and Kulhawy, 1999, the typical values of CoV for cohesion are in the range of 10%–50%). However, the assumed normal distribution for cohesion and the linear dependence of the bearing capacity on the cohesion value in the considered case (for random field, this linear dependence applies to the field-constant multiplier of the cohesion value in each element) should allow the problem to be scaled in terms of SD and CoV (within the δ_c range allowed for the normal distribution). Therefore, the results obtained in this study for δ_c = 10% can be scaled to the other values of δ_c (e.g., δ_c = 20%); to do this, the values of the bearing capacity CoV obtained here should be multiplied by the ratio between new chosen value of δ_c and 10% (in the case of δ_c = 20% by two).

The autocorrelation function describing the internal correlation structure of the field was assumed to be a two-dimensional Markov function, that is, (1) ρ(τv,τh)=exp(−2|τv|θv+−2|τh|θh) \rho \left( {{\tau _v},{\tau _h}} \right) = \exp \left( {{{ - 2\left| {{\tau _v}} \right|} \over {{\theta _v}}} + {{ - 2\left| {{\tau _h}} \right|} \over {{\theta _h}}}} \right) where τ_v and τ_h are distance lags and θ_v and θ_h are SOF values in the vertical and horizontal directions, respectively. In each series of calculations, the values of θ_v and θ_h were assumed differently (the horizontal SOF always being larger), which corresponds to the anisotropic random field structure. The description of the numerical model and the SOF values considered will be presented in the following section.

The boundary value problem of the bearing capacity of the strip foundation was defined and solved using the Abaqus FE code under the assumption of a plane strain. The FE model used is presented in Fig. 1. The model consists of two parts: a 10.8 m × 3.6 m soil domain and a 2 m × 0.5 m concrete strip footing. These dimensions were verified in multivariate numerical simulations (including numerous realisations of the cohesion random field) as sufficient to determine the bearing capacity for the considered case (plastic zones for none of the simulations performed were near the edge of the domain). The contact between the footing and the soil has been assumed to be perfectly smooth. To discretise the continuum model, the soil domain has been divided into 54 ×18 squared FEs and the footing into 10 × 3 FEs, resulting in an element of size h = 0.2 m. The second-order (eight nodes) elements were used. This relatively coarse mesh was adopted (based on several tested meshes) as a compromise between accuracy and numerical efficiency. The boundary conditions are as follows: the bottom of the soil domain is fixed in both directions, while on the sides of the domain, only horizontal displacements are fixed, allowing the nodes for a vertical displacement. It has been assumed that the footing is elastic with Young's modulus E = 32 GPa and Poisson's ratio ν = 0. The soil model was assumed to be elastic–perfectly plastic, with the Mohr–Coulomb plasticity function. Both internal friction ϕ and dilation angle ψ were assumed to be 0º (Tresca). Young's modulus and Poisson's ratio values for the soil were assumed to be E = 200 MPa and ν = 0.33, respectively. Please note that these values should not affect the obtained bearing capacity (cf. Puła and Zaskórski, 2015). As mentioned above, the mean cohesion value was assumed to be 20 kPa and its SD was 2 kPa. All these values are summarised in Table 1.

Figure 1

During the computations, six positions of the vertical force P acting on the foundation were considered, one at the central point of the foundation and the other five eccentric to the foundation with eccentricity values from 0.1 to 0.5 m with an interval 0.1 m (see Fig. 1). Bearing capacity calculations were performed as displacement controlled; the boundary condition for vertical displacement was applied to the bottom node of the foundation in the position that corresponded to the location of the load, and the displacement of the node (towards the soil) increased from 0 to 5 cm. During the displacement of the foundation, the reaction in the displaced node increased until a critical value was reached, which was considered the bearing capacity of the foundation. Due to some differences in results between models using small- or large-strain frameworks, and since the assumed displacement was relatively large and the computational time for both cases was similar, a large-strain framework was utilised in calculations. For the same reason, the possibility of large sliding was assumed for the foundation–soil contact.

Before running the Monte Carlo simulations (MCSs), some preliminary deterministic computations were performed for all the values of eccentricity considered. In these calculations, a constant cohesion value equal to the mean value μ_c = 20 kPa was assumed for all soil modelling elements. The obtained results in the form of function of reaction value in the displaced node with respect to the displacement value are shown in Fig. 2. The obtained bearing capacities Q (maximum values of P) are also shown in the figure.

Figure 2

As seen, the reaction after reaching the maximum value exhibits some small decreases, which can be an effect of a large sliding assumption. The peak decreases obtained for larger eccentricities (which also appear after reaching the maximum value) are related to the passing of the foundation and the soil nodes (during larger foundation slides).

To verify the correctness of the numerical computations performed, the model was validated based on Prandtl's solution of the strip-bearing capacity problem (for cohesive and weightless soil), that is, (2) Q=(2+π)c B' Q = \left( {2 + \pi } \right)c\;B' and Meyerhof theory (1953) assuming a linear decrease of the effective foundation width B’ with eccentricity e, that is, (3) B'=B−2e B' = B - 2e The bearing capacities calculated from the above equations analytically Q_an for B = 2 m and c = 20 kPa were compared with the values determined numerically Q_FEM (presented in Fig. 2) in Table 2. The table also shows the relative differences between the numerical method and the analytical solution. As seen, the respective values are close to each other, although for larger eccentricity values, the differences between the values obtained are significant.

It is worth mentioning that the bearing capacity obtained numerically may correspond to a footing width B* larger than the modelled width B. This is because the displacement distribution under the foundation does not end with the footing boundary, but extends to the first elements beyond the compressed zone on each side. This effect, shown in Fig. 3, although rarely taken into account, can (especially when using larger elements) lead to relative errors of several percentages. For the case without eccentricity and linear distributions of displacement in FEs (or FD zones), the equivalent strip footing width B* corresponding to the results obtained can be estimated as larger than the modelled one by half of the width of element adjusted to the last displaced node on both sides (Itasca, 2011). For second-order elements, the increase in equivalent width B* should be smaller, about one-fourth of the element width h on both sides (since the middle node is not displaced). However, due to the nonlinear interpolation function, the exact value is more difficult to estimate. One way of its estimation would be to compare Eq. 2 to the value obtained numerically for zero eccentricity (assuming that the width of the strip footing is unknown). The width B*, obtained in this way, that is, 2.113 m, is quite close to B + 2×1/4h = 2.10 m, and thus seems probable. However, it is possible that some of the differences between the numerical and theoretical models are not due to the B* value, and the real B* value should be different. The results and relative differences obtained for the mentioned value of B* are also summarised in Table 2. The small values of the differences between the analytical and numerical models seem to confirm the correctness of the assumptions made.

Figure 3

After the validation of the model, the actual MCSs were carried out. The bearing capacity for each of the six-load eccentricity values considered was analysed for eight pairs of SOF values assumed to generate a cohesion random field. These pairs were combinations of the two considered vertical SOF values θ_v, that is, 0.5 and 1.0 m, and the four considered horizontal SOF values θ_h, that is, 2, 5, 10 and 30 m. For both SOFs, these are typical values from the literature covering both weak and strong anisotropy of the subsoil (cf. Lloret-Cabot et al., 2014; Bagińska et al., 2018; Kawa et al., 2019; Kawa et al., 2020). The random field was generated using the FSM algorithm (Jha and Ching, 2013b), which (in the case of rectangular element meshes) enables generation of local averages of the field values over the area of the elements. The ratio between the element dimension h, that is, 0.2 m, and the smallest modelled SOF, that is, 0.5 m, corresponds to the values adopted in the literature (cf. Set et al., 2018; Kawa et al., 2021), and thus should enable a reliable reconstruction of the random field. For each pair of SOF and eccentricity value, N = 1000 MCSs were performed. The cohesion values obtained in the exemplary realisation of the random field for θ_v = 1 m and θ_h = 10 m are presented in grey scale in Fig. 4 (the white colour represents the high values of cohesion). The cumulative equivalent plastic strain values obtained for the same realisation for two different eccentricity values (0.0 and 0.5 m) are shown in Fig. 5. The probabilistic analysis of the obtained results will be presented in the next section.

Figure 4

Figure 5

As mentioned in the previous section, the analysis included eight pairs of SOF values and six eccentricity values, that is, 48 problems in total. For each of these problems, 1000 MCSs with different realisations of random field modelling cohesion were carried out. The time needed to solve a single series of 1000 MCSs on a modern PC is about 24 hours.

The bearing capacity values obtained for each problem were collected and analysed probabilistically. First, the assumed normal distribution of the results was estimated and tested with the Kolmogorov–Smirnov goodness-of-fit test. The test results obtained were in the range of 30%–70%. Exemplary probability density functions (PDF) for the estimated normal distribution for θ_v = 1 m and θ_h = 10 m and six eccentricity values are shown in Fig. 6.

Figure 6

The mean values μ_Q, SDs σ_Q and CoVs δ_Q of the bearing capacity for all considered cases are shown in Figs 7–9, respectively. The mean and CoV values are also summarised in Table 3. As can be observed, the mean values μ_Q for a given value of the eccentricity of the load and different values of the horizontal SOF are almost identical and similar to those obtained in the deterministic calculations. On the contrary, SD σ_Q and CoV δ_Q significantly increase with increasing horizontal (and comparing the values for θ_v = 0.5 m and θ_v = 1 m also with vertical) value of SOF. However, both of these quantities stabilise for higher values of horizontal SOF, say θ_h greater than 10 m. This type of behaviour is often reported in the literature for probabilistic analysis of bearing capacity problems involving the use of random fields (e.g., Pieczyńska-Kozłowska et al., 2015; Kawa and Puła, 2020).

Figure 7

Figure 8

Figure 9

As seen in the figures, for small eccentricity values, the stabilisation of the SD of the bearing capacity is slower than in the case of larger eccentricity values. The graphs also show that as the eccentricity increases, not only the mean value but also the SD of the bearing capacity decreases. However, the resulting CoV increases with eccentricity. All these effects are due to a reduction in the effective width of the foundation. This reduction is associated with the smaller size of the failure mechanism, and thus a smaller averaging area of the cohesion values (cf. Chwała and Kawa, 2021), which, according to local averaging theory (Vanmarcke, 2010), means a smaller reduction in variance. For larger eccentricities, this results in both faster stabilisation of SD (when the mechanism size is smaller than SOF, the variance reaches its maximum value) and relatively larger variation in bearing capacity (smaller reduction of variance). Although easy to explain, the increase in CoV of the bearing capacity δ_Q is important, as it can mean lower reliability for eccentrically loaded foundations.

An interesting behaviour regarding small eccentricity values can be observed in Fig. 8. For both the considered vertical SOF values (Fig. 8a, b), the graph for e = 0.0 m intersects the graph for e = 0.1 m. However, this behaviour is not reflected in the CoV graph. The specific behaviour of these characteristics for small eccentricities can be observed in more detail in Figs 10 and 11, where the effect of eccentricity on the values of SD σ_Q and CoV δ_Q, respectively, is shown. As seen in the figures, for small values of eccentricity, both graphs behave non-linearly, which is more visible in the SD graph. This effect is interesting, but needs further investigation.

Figure 10

Figure 11

In the last step of the analysis, the effect of eccentricity on the value of reliability index β and the probability of failure p_f was investigated. These calculations were performed for only one pair of SOFs, that is, θ_v = 1 m and θ_h = 10 m, assuming that the estimated distributions presented in Fig. 6 are exact (since only 1000 realisations were performed, this cannot be the case; cf. Kawa and Puła, 2020; Kawa et al., 2021). The load P value for the reliability analysis was assumed to be such that according to EN 1990 (2002), it exactly meets the ultimate limit state reliability condition for a typical structure (β = 3.8) in the case without eccentricity. For the SOFs considered, this value was calculated based on data presented in Table 3 as 159.95 kN/m. In addition to the eccentricity values considered earlier, the reliability analysis also included smaller eccentricities, that is, e = 0.01 m and e = 0.02 m. Since the numerical analysis for these eccentricities, due to the location of the displaced node, would have to be performed for a finer mesh, which is associated with a significant increase in the computation time, the mean μ_Q and CoV δ_Q for these eccentricities were interpolated (linearly) based on the μ_Q and δ_Q obtained for e = 0.0 m and e = 0.1 m. This interpolation may raise some doubts because, as mentioned above, for small eccentricities, the graphs for SD and CoV of bearing capacity are non-linear. However, for CoV, the impact of this non-linearity is small, especially for some pairs of fluctuation scales, including the one adopted for the analysis, namely θ_v = 1 m and θ_h = 10 m. The interpolated values for this pair of SOFs are also shown in Table 3. The values of the reliability index β and the probability of failure p_f determined for the assumed load for all eccentricities considered are summarised in Table 4.

As seen in the table, an increase in eccentricity causes a significant decrease in the reliability index β. The small values of this index or even its negative values (p_f greater than 50%) obtained for larger eccentricities are not surprising; this is the effect of the mentioned strong linear decrease of μ_Q. However, it is worth noting that an increase in eccentricity of even 0.01 m results in a significant decrease in the index value. For an eccentricity of 0.02 m, the probability of failure p_f is already three times greater than the allowable value. Eccentricities with such values are typical in all types of constructions and are often assumed even for the model without the bending moments. The analysis shows that even small eccentricity is an important factor significantly affecting structural reliability, and at least small random eccentricities should be taken into account when performing probabilistic analysis of foundation bearing capacity. This is easy to apply using the RFEM approach, but has not been practiced so far. Adding this and perhaps some other elements can be necessary to enable the use of probabilistic analysis in engineering practice.

In this paper, a probabilistic analysis of the bearing capacity of an eccentrically loaded shallow foundation on a purely cohesive soil was carried out. The soil was assumed to be weightless. The cohesion was modelled by an anisotropic random field, for which eight pairs of different vertical and horizontal SOFs were considered. For each of these cases, six different values of load eccentricity were assumed, resulting in a total of 48 problems. For each of these problems, N = 1000 MCSs were performed. The paper analyses the effect of the SOFs and eccentricity on the probabilistic characteristics of the bearing capacity results obtained. In the final part of the paper, the effect of eccentricity on the reliability index value (and failure probability) of the considered structure for the selected case of SOFs is shown. The following conclusions can be drawn from the paper: i)

Risk management for geotechnical structures requires modelling that includes the effect of spatial variability of soil parameters. A frequently used approach in this case is the RFEM, which allows the modelling of various structures located in soil, including strip footings (which were among the first to be modelled using this method). Although in typical structures, the loading of foundations almost always occurs with a certain eccentricity, this has, so far, rarely been taken into account in probabilistic studies. To the authors’ knowledge, this paper is the first study to analyse the bearing capacity of eccentrically loaded strip footings on a soil modelled by an anisotropic random field.