Progressive failure analysis of shallow foundations on soils with strain-softening behaviour

This paper presents a finite element approach to analyse the response of shallow foundations on soils with strain-softening behaviour. In these soils, a progressive failure can occur owing to a reduction of strength with increasing the plastic strains induced by loading. The present approach allows this failure process to be properly simulated by using a non-local elasto-viscoplastic constitutive model in conjunction with a Mohr–Coulomb yield function in which the shear strength parameters are reduced with the accumulated deviatoric plastic strain. Another significant advantage of the method is that it requires few material parameters as input data, with most of these parameters that can be readily obtained from conventional geotechnical tests. To assess the reliability of the proposed approach, some comparisons with experimental results from physical model tests are shown. A fairly good agreement is found between simulated and observed results. Finally, the progressive failure process that occurs in a dense sand layer owing to loading is analysed in details, and the main aspects concerning the associated failure mechanism are highlighted.


Introduction
In the current applications, the bearing capacity of shallow foundations is usually evaluated using the well-known equation proposed by Terzaghi [1], i.e.
in which q lim is the bearing capacity pressure, q is the uniformly distributed surcharge replacing the overburden soil at the level of the foundation base, c 0 is the cohesion intercept of the soil, B is the foundation width, c is the soil unit weight, N q , N c and N c are the bearing capacity factors which depend on the soil shearing resistance angle, u 0 . Eq. (1) refers to strip footings resting on homogeneous and dry soil with centrally located vertical load and symmetrical failure pattern. To extend Terzaghi's solution to more general conditions than those above specified, a great number of theoretical studies were conducted in which numerical techniques were used to account for reliably the effects of important factors on the bearing capacity calculation, such as footing shape, roughness of base, inclination and eccentricity of loading, ground surface inclination, groundwater conditions and other factors . In most of these studies, it was assumed that failure occurs simultaneously along the slip surfaces that develop in the soil mass beneath the footing. However, this failure process is progressive in nature owing to the fact that the plastic strains induced in the soil by loading are markedly non-uniform. As a consequence, the soil shear strength is not simultaneously mobilised at all points of a slip surface. It was also assumed that the soil strength parameters remain unchanged even if large strains are induced by loading. This assumption is inadequate for soils that are characterized by a pronounced strain-softening behaviour, such as dense sands. In reality, in these materials it occurs that some portions of the soil first fail owing to loading, with the shear strain that is located in a zone of limited thickness (shear band). With increasing strain within this zone, soil strength reduces from peak towards the critical state. Owing to the consequent redistribution of stress, the shear band propagates in the soil and a slip surface progressively develops up to causing the collapse of the soil-foundation system. At failure, the average strength mobilised along the slip surface is generally less than the peak strength and greater than the strength at the critical state of the sand under consideration. The occurrence of a progressive failure in loading tests and centrifugal tests concerning footings on granular soils, was observed by several authors [27][28][29]. On the basis of these evidences, Perkins and Madson [29] proposed a semi-empirical procedure to estimate the bearing capacity of footings on dense sands. However, a progressive failure process can be successfully predicted using an approach that accounts for properly the strain-softening behaviour of the soil and it is able to simulate reliably the formation and development of shear zones within the soil. Finite element approaches with these characteristics were employed by Siddiquee et al. [30] and Banimahd and Woodward [31] to analyse the response of footings on granular soils to loading. Generalised elasto-plastic strain-softening models are incorporated in these approaches for modelling the behaviour of the soil involved. However, these models generally require that a significant number of specific constitutive parameters is determined.
In this study, a finite element approach is proposed to analyse the response of shallow foundations resting on soils with strain-softening behaviour. This approach utilises a non-local elasto-viscoplastic constitutive model in conjunction with a Mohr-Coulomb yield function in which the strength parameters are reduced with the accumulated deviatoric plastic strain. The proposed method requires few material parameters as input data. In addition, most of these parameters can be readily obtained from conventional geotechnical tests. The method is first applied to simulate the experimental results from some physical model tests concerning a rigid strip footing resting on a layer of dense sand in which a thin layer of weak material is located at different depths from the ground surface. Then, the progressive failure process that occurs in a homogeneous sand owing to loading is analysed, and the main aspects of this process are discussed.

Non-local elasto-viscoplastic model
Under the assumption of small strains, the total strain rate tensor for an elasto-viscoplastic material can be written as follows: where _ e ij is the total strain rate tensor, and _ hk is the effective stress rate tensor, and C ijhk is the elastic compliance tensor.
In the present study, the viscous component of the model is used with the primary objective of regularising the numerical solution, as it was suggested by Zienkiewicz and Cormeau [32]. Following Perzyna [33], the viscoplastic strain rate tensor is expressed by the following equation: in which U is the viscous nucleus that depends on the yield function F, and m ij is the gradient to the plastic potential function Q (i.e., m ij ¼ @Q =@r 0 ij ). The gradient of Q defines the direction of the viscoplastic strain rate tensor, and the yield function influences the modulus of this tensor by means of U. In this connection, di Prisco and Imposimato [34] proposed the following relationship for U: where c and a are constitutive parameters, and p 0 denotes the mean effective stress. A maximum value of 3 should be assumed for the product aF to prevent the exponent in Eq. (5) from becoming excessively large [36]. The parameter c influences the strain rate and consequently the rapidity with which strain occurs owing to a given stress increment. In particular, strain rate increases with increasing the value of c [36]. The Mohr-Coulomb failure criterion is adopted in the present study to describe the yield function F. Referring to cohesionless soils, the expression of F in terms of the principal effective stresses is where r 0 1 and r 0 3 are the major and minor principal effective stresses respectively (it is assumed that compressive stress is positive). Following several authors [18,[37][38][39][40][41], the strain-softening behaviour of the soil is simulated by reducing the soil shearing resistance angle u 0 from the peak value, u 0 p , to that at constant volume, u 0 cv , with increasing the accumulated deviatoric plastic strain. In the present study, this strain is expressed by the parameter k shear that is defined as follows [42]: t is time, and _ e p ij is the deviatoric plastic strain rate tensor the expression of which is in which d ij is the Kronecker tensor, and _ e p ij is the plastic deformation rate tensor. For simplicity, the relationship considered in the present study to relate the mobilised shearing resistance angle to k shear is schematised in Fig. 1a. As can be seen, this relationship is defined by two thresholds denoted as k p shear and k cv shear , respectively. The flow rule is of non-associated type with the plastic potential function which is expressed as in which w denotes the dilatancy angle of the soil. A relationship similar to that shown in Fig. 1a is also considered for w (Fig. 1b).
A significant advantage of this constitutive model is that it requires few soil parameters as input data, in comparison with other more sophisticated models. Specifically, these parameters are Young's modulus E 0 , Poisson's ratio v 0 , the parameters defining the shear strength of the soil (i.e., u 0 p and u 0 cv ), the dilatancy angle w P at peak, the strain thresholds, k p shear and k cv shear , and the viscous parameters c and a. Generally, the most part of these parameters are directly obtained from triaxial tests. In addition, considering the role principally attributed to the viscous component of the constitutive model (i.e., to regularise the numerical solution [32]), the viscous parameters c and a along with k p shear and k cv shear can be determined by matching the experimental results from triaxial tests with those obtained simulating numerically these tests, taken loading rate into account.
The described elasto-viscoplastic model could however lead to mesh-dependent results when materials with strain-softening behaviour are considered [34,35,43,44]. To overcome this drawback, di Prisco and Imposimato [34,35] proposed a non-local definition for U, according to which Eq. (4) is replaced by the following equation: where x i defines the position of a generic point with spatial co-ordinates (x, y, z), V is a spherical volume with centre in x oi and radius R, Eqs. (12) and (13) introduce an additional constitutive parameter that is the radius R of the spherical volume V. The choice of R depends on the problem under consideration. For example, a small value of R is usually selected when strain localisation is analysed at the scale of the soil specimens tested in the laboratory. In these cases, R may be empirically correlated to the grain size of the material. In this connection, Mühlhaus and Vardoulakis [45] suggested to assume a value of R of order of 20 times the particle diameter. The value of R also determines the thickness of the shear band. In particular, this thickness is about 2R for failure processes under two-dimensional conditions [46,47]. The described constitutive model was successfully used for performing stability analyses of slopes in soils with strain-softening behaviour [39,40]. This model is implemented in the finite element code Tochnog [42] which was used in the present study to perform the numerical analyses presented in the subsequent sections.

Analysis of physical model tests
Muscolino [48] carried out a series of physical model tests under 1 g plane strain conditions using a rigid footing placed on a dry sand layer in which a thin layer of weaker material was located at different depths from the surface. The main objective of Muscolino's study was to analyse the influence of this thin layer on the bearing capacity of shallow foundations. Further experimental results from these model tests can be found in Valore et al. [49]. A scheme of the tests is shown in Fig. 2 in which the dimensions of the soil layer and those of the footing are indicated. The sand was deposited in a box of plexiglas reinforced by a steel frame. The inner sides of this box were covered by a glass sheet smeared with oil to minimize the effect of friction at the interface between the soil and the box walls. The weak layer was 3 mm thick and consisted of wet bentonite or dry talcum. A uniformly distributed load p = 2.9 kPa was initially applied to level the layer surface. This load was then removed before placing the footing which was constituted by a block of aluminium with rectangular cross-section. Owing to the fact that a sheet of sand-paper was bonded at the lower surface of the footing (i.e., at the surface at contact with the soil), the foundation can be considered rough. The footing was placed at the centre of the box, and subsequently it was displaced downwards at a displacement rate of 0.6 mm/min. The vertical load and the associated settlement of the footing were measured during the tests. In addition, the displacement field and the slip surfaces which developed in the soil owing to loading were documented by digital photographic shots.
The soil used in the tests was a quartz sand with an average value of the particle diameter, d m = 1 mm and a unit weight of c = 16.8 kN/m 3 . The strength parameters were obtained by Muscolino [48] from the results of direct shear tests. These tests showed that the sand was characterized by an evident strain-softening behaviour [49]. On the basis of the available experimental results, a peak shear resistance angle u 0 p ¼ 50 with an effective cohesion equal to zero can be assumed for defining the soil strength at low effective stress (<60 kPa) [49]. In addition, it can be assumed that u 0 cv 33 . The dilatancy angle at peak, w p , was about 17°. Two loading tests in which the weak layer consisted of Luzenac talcum are considered in the present study. In these tests, this talcum layer was located at z = 0.5B and z = 1.2B (with B = footing width) from the free surface, respectively. The strength parameters of this material were u 0 , with effective cohesion and dilatancy angle nil [48]. A Prandtl-type failure mechanism was observed in both tests with a portion of the failure surface that developed in the talcum layer, as it is schematized in Fig. 3. The slope angle of the wedge-shaped zone, a, and the angle, b, with which the slip surface intersects the free surface of the sand layer are also indicated in this figure.
The mesh adopted in the analysis is shown in Fig. 4. Taking advantage of the symmetry of the soil-foundation system a half grid is considered. The base is fully fixed, and the lateral sides are constrained by vertical rollers. The mesh consists of triangular elements with three nodes and one Gauss point. The initial stress state of the the soil is calculated under K o -conditions. To account for the overconsolidation effect caused by the application and subsequent removal of the load p at the layer surface, the coefficient of earth pressure at rest is evaluated using the relationship K o ¼ ð1 À sin u 0 ÞOCR n in which OCR is the overconsolidation ratio and n is an empirical exponent [50]. Values of K o equal to 0.38 Fig. 2. Scheme of the model tests performed by Muscolino [48]. and 0.85 were assumed for the sand and talcum, respectively. To simulate a rigid footing, a uniform vertical displacement with the same rate used in the tests (0.6 mm/min) is imposed at the nodes of the soil-foundation interface. In contrast, the horizontal displacement at these nodes is prevented to reproduce a rough footing.
The non-local elasto-viscoplastic strain-softening model described in the previous section is used to model the behaviour of the sand which, as already said, experienced a strain-softening behaviour during the laboratory tests. The same constitutive model but with no strain-softening is adoped for the talcum. The soil parameters used in the analysis are indicated in Table 1. In addition, the values assumed for a and c are 61 and 2 Â 10 À6 s -1 , respectively. A part of these parameters (c, u 0 p , u 0 cv , w p ) was obtained experimentally by Muscolino [48] and Valore et al. [49]. The parameters which are not available (i.e., v 0 , E 0 , k p shear , k cv shear , c and a) are evaluated in this study by matching the experimental results from the loading tests with those obtained simulating these tests by the proposed numerical approach, on a trial and error basis. Lastly, owing to the fact that Tochnog allows an only value of R to be used, it is assumed that R = 5 mm for both sand and talcum. A comparison between experimental and theoretical results is presented in Figs. 5a and 6a in terms of the settlement curve which relates the foundation settlement normalized with respect to its width (s/B) to the associated footing pressure. There is a fairly good agreement between observation and simulation, although the calculated settlement at peak is slightly greater than that measured during the tests and some discrepancies between experimental and theoretical results can be observed in the post-peak relationship. However, the calculated value of the load at peak and that associated with the final settlement essentially coincide with those measured in the tests. In addition, Figs. 5b and6b show the accumulated plastic deviatoric strain (expressed by k shear ) calculated at the end of the simulation, for z = 0.5B and z = 1.2B respectively.
These results show that a failure mechanism of Prandtl type generally occurs, although this mechanism is influenced by the presence of the thin layer of talcum. Specifically, the slip surface starting from the edges of the footing defines a wedge-shaped zone below the loaded area. This surface develops in the underlying sand and in the talcum layer before involving the free surface. These results are consistent with what observed by Muscolino [48] and Valore et al. [49]. Moreover, the calculated values of the angles a andb ( Fig. 3) are in reasonable agreement with those measured during the tests, as it is documented in Table 2. These angles as well as the size of the failure mechanism depend on the location of the thin layer of talcum. In conclusion, the numerical simulation captures the main features of the foundation response observed in the physical model tests.

Analysis of progressive failure
Muscolino [48] also performed a test in which the footing was placed on a homogeneous layer of sand (i.e., without any weak inclusion). Unfortunately, an asymmetrical failure mechanism was observed in this test owing to some drawbacks of experimental nature that unexpectedly occurred (perhaps owing to a not perfectly central load or not completely uniform soil). This test is simulated in the present study using the proposed numerical approach with the primary objective of highlighting the main aspects of the failure process that develops in a strain-softening soil owing to loading. The constitutive parameters assumed for the sand are those validated by the comparisons with the experimental results from the model tests (Figs. 5 and 6 and Table 2). These parameters are specified in the previous section. Fig. 7 shows the settlement curve obtained using a different number of elements to show that the influence of the mesh on the results is not substantial. In Fig. 7, some loading stages (i.e., Fig. 4. Finite element mesh adopted in the present study to simulate the tests performed by Muscolino [48].

Table 1
Soil parameters used in the analyses. A-F) are also indicated to support the presentation of the results. As expected, the bearing pressure at peak and that at the final loading stage are higher than those obtained in the presence of the thin layer of talcum (Figs. 5a and 6a). For the sake of completeness, the value of the bearing capacity at peak is also compared with that predicted by the classical equation  Table 2 Comparison between calculated and observed values of the angles a and b, for z = 0.5B and z = 1.2B.   8. Accumulated plastic deviatoric strain (k shaer ) calculated at any loading stage indicated in Fig. 7.
where the bearing capacity factor N c is evaluated using the expression recently proposed by Salgado [51] and also suggested by Lyamin et al. [26], that is in which it is assumed that the soil shear resistance angle u 0 ¼ 50 The resulting value of q lim is 360 kPa. Considering that the mean effective stress along the slip surface should be of order of q lim /10 [52], this value of q lim is consistent with the selected value of u 0 at low stress (<60 kPa). In addition, the calculated value of q lim is greater than that obtained numerically using the present approach (335 kPa). This result should be ascribed to the occurrence of a progressive failure which causes a non-uniform strain state in the soil and consequently a different mobilization of the shear resistance angle. To document this failure process, Fig. 8 shows the accumulated plastic deviatoric strain (expressed by k shear ) calculated for any loading stage indicated in Fig. 7. At the loading stage A, a restricted area near the edge of the footing experiences plastic strains. With increasing the load, this area spreads downwards into the soil and the induced plastic strains increase (stages B-C). When the peak load is reached (stage C), a wedge-shaped zone develops in the soil beneath the footing along with a zone of radial shear that emanates from the edge of the footing. The soil located within the wedge-shaped zone remains however in a state of elastic equilibrium. Both the wedge-shaped zone and the radial zone are delimited by shear bands. Owing to the accumulated deviatoric plastic strain in these shear bands, the mobilized resistance angle of the soil reduces from the initial value of u 0 p in accordance with the relationship shown in Fig. 1a. In the post-peak stages (stages D-F), the area affected by plastic strains enlarges laterally involving the soil on the side of the footing. The deviatoric plastic strain accumulated in this portion of soil is however much lower than that in the former zones. Considering the values of k shear calculated at the final stage of loading (stage F) and the u 0 À k shear relationship assumed in the analysis (Fig. 1a), it can be concluded that different values of the shear resistance angle are mobilized in the above-mentioned zones. In this connection, Fig. 9 shows the portions of the soil where the angle u 0 cv is mobilized, those where the shearing resistance angle is at peak, and the zones in which the value of u 0 is intermediate between u 0 cv and u 0 p . For completeness, the mobilized values of w are also indicated in Fig. 9. Taking into account the results shown in   cv is the shear resistance angle mobilized in the shear zone that delimites the wedge below the footing, as documented in Fig. 9). On the other hand, the angle b (Fig. 3) is about 33°which is greater than 45 À u 0 P 2 (u 0 p is the resistance angle mobilized in the portion of soil close to the layer surface, as shown in Fig. 9). Following Terzaghi [1], this implies that the lateral compression of the soil in this latter zone is not great enough to produce a passive limit state. The non-uniform field of plastic strain and the associated non-simultaneous mobilization of soil strength account for the progressive failure that occurs in the soil.
The displaced soil mass associated with the loading stages indicated in Fig. 7 also documents the progressive nature of the failure mechanism (Fig. 10). Until peak load is reached, the displacements are concentrated in a small region of soil beneath the footing and are essentially downwards, whereas the soil adjacent to this region is essentially undisturbed (stages A-C). After peak, a different mechanism is activated with the zone affected by displacement that spreads laterally and involves the soil on the side of the footing (stages D-F). Owing to the friction at the soil-footing interface, the central zone beneath the footing behaves as if it is a part of the sinking footing. At the final stage of loading (stage F), the soil on the side of the foundation experiences displacements with a significant vertical component upwards. This is corroborated by the displacement vectors shown in Fig. 11. Considering the results presented by Potts and Zdravković [17] concerning the influence of the dilatancy angle on the failure mechanism of footings, this latter evidence should be also ascribed to the fact that reduction of w practically does not occur in the soil close to the layer surface on the side of the foundation, as it is documented in Fig. 9.

Concluding remarks
A finite element approach in which a non-local elasto-viscoplastic constitutive model in conjunction with the Mohr-Coulomb yield function is incorporated, has been proposed to predict the response of strip footings resting on soils with strain-softening behaviour. This behaviour is simulated by reducing the strength parameters with increasing the accumulated deviatoric plastic strain. Using the proposed approach, the formation and development of shear zones within the soil are reliably simulated and the post-peak response of the soil is reasonably captured. Another advantage of the method is that few soil parameters are required as input data. In addition, most of these parameters can be readily obtained from conventional geotechnical tests.
A fairly good agreement has been found between the theoretical results achieved using the present approach with the experimental results from some physical model tests concerning a rigid and rough strip footing on a dense sand layer in which a thin inclusion of weaker material is located at different depths from the layer surface. The case of a rough strip footing on a homogeneous sand layer has been also analysed. The associated settlement curve experiences a pronounced peak that drastically reduces to the ultimate value. In addition, the results have clearly shown the nature progressive of the failure process that occurs in the soil with a Prandtl-type mechanism. In particular, a wedge-shaped zone and two radial zones first develop in the soil beneath the loaded area. These zones are delimited by shear bands with very high values of the accumulated deviatoric plastic strain. Then, significant plastic strains are also induced in the outer zones of the soil on the sides of the footing. However, these strains are much lower than those generated in the former zones of the soil. Considering the strain-softening behaviour of the soil involved, this non-uniform field of shear strain and the consequent non-simultaneous mobilization of the strength parameters account for the occurrence of a progressive failure in the soil.