The use of the node-based smoothed finite element method to estimate static and seismic bearing capacities of shallow strip footings

H.C.Nguyen ,T.Vo-Minh

bDepartment of Civil and Environmental Engineering,Imperial College London,UK

b Department of Civil Engineering and Industrial Design,University of Liverpool,Liverpool,UK

c Faculty of Civil Engineering,Ho Chi Minh City University of Technology(HUTECH),Ho Chi Minh City,Viet Nam

Keywords:Limit analysis Node-based smoothed finite element method(NS-FEM)Second-order cone programming(SOCP)Seismic bearing capacity Strip footing

ABSTRACT The node-based smoothed finite element method(NS-FEM)is shortly presented for calculations of the static and seismic bearing capacities of shallow strip footings.A series of computations has been performed to assess variations in seismic bearing capacity factors with both horizontal and vertical seismic accelerations.Numerical results obtained agree very well with those using the slip-line method,revealing that the magnitude of the seismic bearing capacity is highly dependent upon the combinations of various directions of both components of the seismic acceleration.An upward vertical seismic acceleration reduces the seismic bearing capacity compared to the downward vertical seismic acceleration in calculations.In addition,particular emphasis is placed on a separate estimation of the effects of soil and superstructure inertia on each seismic bearing capacity component.While the effect of inertia forces arising in the soil on the seismic bearing capacity is non-trivial,and the superstructure inertia is the major contributor to reductions in the seismic bearing capacity.Both tables and charts are given for practical application to the seismic design of the foundations.

1.Introduction

The bearing capacity of strip footings under static conditions has been extensively studied by Prandtl(1920),Terzaghi(1943),Meyerhof(1951),Hansen(1970),Vesic(1973),and many others.In recent years,the effect of horizontal earthquake body forces on shallow foundation bearing capacity is an important topic of research in geotechnical engineering.The first studies were describing the seismic bearing capacity of shallow foundations concerning the works of Meyerhof(1953,1963),where the seismic forces were applied at the structure only as inclined pseudo-static loads.However,in these solutions,the inertia of the soil mass is not included.In the presence of seismic forces,the available theoretical researches are mainly based on(i)the limit equilibrium method;(ii)the upper and lower bound limit analysis;and(iii)the method of stress characteristics(slip-line method).

Using the limit equilibrium method,Sarma and Iossifelis(1990)determined the seismic bearing capacity of strip footings considering the failure mechanisms of three zones as the unsymmetrical active wedge,internal shear zone and passive wedge.Later on,Budhu and Al-Karni(1993)investigated an asymmetrical failure surface for seismic analysis of shallow foundations with similar three zones as that Vesic(1973)considered.Choudhury and Subba Rao(2005)considered two mechanisms,log-spiral curve,and planar active wedge and passive wedge,to calculate the seismic force on the structure and soil below.Recently,Saha and Ghosh(2015)investigated the seismic bearing capacity of strip footings using pseudo-dynamic coupled with the limit equilibrium method.

Based on the upper bound theorem,Richards et al.(1993)used two triangular planar wedges in the failure surface,each on the active and passive zones,to investigate the effect of soil and superstructure inertia on the seismic bearing capacity of strip footings.Dormieux and Pecker(1995)and Soubra(1999)investigated the seismic bearing capacity of strip foundation on the horizontal using pseudo-static method.Ghosh(2008)and Saha et al.(2018)used a pseudo-dynamic approach to explore the effect of shear and primary wave velocities and ampli fication factor on the seismic bearing capacity of strip footings using the upper bound limit analysis.More recently,Kalouzari et al.(2019)investigated the seismic bearing capacity of strip foundations in the vicinity of slopes using the lower bound finite element method(FEM)in conjunction with the linear programming technique.

Using the method of characteristics,Kumar and Mohan Rao(2003)and Kumar(2003)presented the variation of the seismic bearing capacity with changes in horizontal seismic coef ficients for different slope inclinations.Although seismic bearing factors reduce as the horizontal seismic acceleration increases,the evolution of characteristics line networks depends on soil and superstructure inertia in a complex manner.Recently,Cascone and Casablanca(2016)used the method of characteristics to compute the bearing capacity factors of strip footings and evolutions of failure patterns with the horizontal seismic acceleration for each seismic case.The effects of the soil inertia force and the superstructure inertia with variation in horizontal seismic coef ficients are distinguished in the analysis.

In recent decades,the FEMplays an important role in designing many practical engineering structural systems.Due to the simplicity,the 3-node linear triangular elements(FEM-T3)are popular.One of the notable drawbacks of FEM-T3 elements is the volumetric locking phenomenon when the mesh becomes highly distorted.In the standard FEM,the stiffness matrix is used to map basis function derivatives from the reference element to the natural element in the mesh.When distorted meshes are used in the FEM,the Jacobian becomes ill-conditioned,affecting the accuracy of the method.

To overcome this,Liu et al.(2007a,b,2009a,b)have combined the strain smoothing technique used in meshfree methods(Chen et al.,2001)into the FEM to formulate a series of smoothed FEMs used in the field of solid mechanics.The basic idea of strain smoothing based on the domain integration on the cells becomes line integration along the boundary of the cell.In the node-based smoothed FEM(NS-FEM),the strain smoothing technique over the cells associated with nodes is computed directly using only the shape functions themselves(not their derivatives),which eliminates the volumetric locking phenomenon.As a result,no coordinate transformation is required in the NS-FEM,and the problem domain can be discretized in highly distorted meshes.NS-FEMhas been successfully applied to several fields,including structural mechanics,solid mechanics,acoustic analysis,and electromagnetic problems.

Together with NS-FEM,other smoothed FEMhave been applied to computation of the collapse loads;and numerical results reveal that the smoothed FEMs are naturally“immune”from the volumetric locking.In particular,Le et al.(2010)presented an upper bound procedure using cell-based smoothed FEM,arriving at conclusions that the use of the strain smoothing technique can overcome the volumetric locking and provide accurate results of the collapse loads with minimal computational effort.In addition,Nguyen-Xuan et al.(2012)stated that combination of the nodebased FEM and primal-dual algorithm serves as a better means of calculating the limit and shakedown loads of structures when compared with the available solutions.Application of the smoothing technique over the edge of triangle elements to the kinematic limit analysis was presented by Le et al.(2013),in which the volume locking is removed for the plane stress and strain problems when using only one Gauss point for each smoothing domain.More recently,the cell-based smoothed FEM is applied to limit and shakedown analyses of structures(Ho et al.,2019)and the collapse of soil(Le,2017).

Recently,the application of the smoothed FEM to the upper bound limit analysis in geotechnical problems has been developed by Nguyen et al.(2011),Nguyen(2020,2021a,b),Vo-Minh et al.(2017,2018),Vo-Minh(2020),Vo-Minh and Nguyen-Son(2021),and Nguyen and Vo-Minh(2022).The Mohr-Coulomb yield criterion can be formed in a second-order cone programming(SOCP)using the NS-FEM in geomechanics problems.The problem can then be solved by the primal-dual interior-point method implemented in the MOSEK software package(MOSEK ApS,2009).This algorithm was proved to be a very effective optimization tool for limit analysis in geotechnical engineering.In order to assess the effects of the earthquake on the stability,the so-called pseudostatic approach has been widely used because of its simplicity to be implemented in the numerical procedure and its ability to produce acceptable solutions to plasticity problems in geotechnical engineering under seismic conditions(Ausilio et al.,2000;Sahoo and Kumar,2012,2014;Chakraborty and Kumar,2013;Krabbenhoft,2018).Seismic effects on the bearing capacity are represented by considering relevant inertia forces which are proportional to seismic accelerations in both the vertical and horizontal directions.Installing these forces in the simulations is straightforward as noted in Ausilio et al.(2000),Sahoo and Kumar(2012,2014),Chakraborty and Kumar(2013)and Krabbenhoft(2018),and is further discussed in the following section for each seismic bearing capacity problem.The dynamic problems turn now to solve the standard static problems which are augmented by inertial forces due to the soil weight,surcharge and the structure.

It is considered important that the effects of cyclic loads cannot be explored using limit analysis,thus the shakedown analysis(Nguyen-Xuan et al.,2012)or the dynamic analysis(Wang et al.,2021a,b)should be performed to further understand the bearing capacity of soils under cyclic loads.One of the main features of soil under cyclic loading is the accumulation of strain or deformation during a great number of cycles that result in signi ficant reduction of behaviors(Abadie,2015;Houlsby et al.,2017;Abadie et al.,2019).Under cyclic loads and seismic forces,due to rapid changes in shaking direction and amplitude during earthquake,the available soil shear strength under a foundation may repeatedly and momentarily be attained,inducing several instantaneous failures.In addition,due to the reduction of footing-soil contact area,the accumulated permanent rotation of the foundation reduces the bearing capacity of strip footings.The seismic loads can be represented by pseudo-static loads using corrective parameters,which are further discussed in the analysis by changing the direction of seismic acceleration.

In this paper,we present the upper bound limit analysis using the NS-FEMto estimate the seismic bearing capacity factors of strip footings based on the pseudo-static approach.In general,increasing the horizontal seismic acceleration of inertia force and superstructure inertia reduces the seismic bearing capacity factors of strip footings.The corrective coef ficients are presented to point out the reduction in the seismic bearing capacity factors due to the effects of soil and superstructure inertia.The numerical results obtained by the NS-FEM are compared with those from Cascone and Casablanca(2016)and previous studies to verify the proposed method’s accuracy and reliability.

The paper is arranged as follows.Section 2 describes a brief review of the NS-FEM for the upper bound limit analysis in the geomechanical problems.In Section 3,several numerical examples are performed and discussed to demonstrate the effectiveness of the NS-FEM.Finally,some concluding remarks are made in Section 4.

2.The NS-FEM for upper bound limit analysis of the geomechanical problems

2.1.A brief overview of the NS-FEM

where N(s)is the set containing nodes directly connected to nodek,dkis the nodal displacement vector,and the smoothed strain gradient matrixon the domaincan be determined as

2.2.The upper bound limit analysis for plane strain geotechnical problems using the NS-FEM

A two-dimensional(2D)problem domainΩbounded by a continuous boundaryS˙u∪St=SandS˙u∩St=?is considered.The rigid-perfectly plastic body is subjected to external traction load g onStand body force f on the boundarySuprescribed by the displacement velocity vector˙u.The strain rate can be expressed by

In the upper bound theorem,for a kinematically admissible displacement field˙u∈U,where U is a space of kinematically admissible velocity field,we have whereα+is the limit load multiplier of the external traction load g and body force f.

The external work can be determined as

The internal plastic dissipation of the 2D domainΩcan be written as

in which the space of kinematically admissible velocity field U is denoted by

De fining C={˙u∈U｜Wext(˙u)=1},the limit analysis problem is based on the kinematical theorem to determine the collapse multiplierα+yielding the following optimization problem:

For plane strain geotechnical problems,Makrodimopoulos and Martin(2007)proposed the internal plastic dissipation equation as follows:

wherecand?are the cohesion and friction angle of the soil,respectively.

Fig.1.The smoothing cells associated with the nodes in the NS-FEM.

whereμis a non-negative plastic multiplier,and the Mohr-Coulomb yield functionψ(σ)can be expressed in the form of stress components as

Using the NS-FEM,the problem is discretized byNetriangular elements and the total number of nodes isNn.The smoothed strain rate?˙εcan be calculated from Eq.(1).The upper bound limit analysis for plane strain geomechanics problems using the Mohr-Coulomb failure criterion can be written as

subject to

whereα+is the collapse load(i.e.in this study,α+is either the static or seismic bearing capacity factor),andAiis the area of nodei.The last constraint in Eq.(14b)is expressed in the conic form.As a result,the conic interior-point optimizer of the academic MOSEK package(MOSEK ApS,2009)is used for solving this problem.The upper bound using the NS-FEMhas been written using the Matlab language.The computations were performed on a Dell Optiplex 990(Intel Core?i5,1.6 GHz CPU,8 GB RAM)in a Windows XP environment.

3.Computation of static and seismic bearing capacity factors

3.1.Effect of soil inertia on seismic bearing capacity factors

Based on the expression for the bearing capacity proposed by Terzaghi(1943),the ultimate load of strip footings under seismic conditions that only consider the inertia forces arising in the soil is conveniently expressed as follows:

where the seismic bearing capacity factorsonly account for the effect of soil inertia.

Fig.2 illustrates the geometry and soil parameters to investigate the effect of soil inertia on the bearing capacity factors of strip footings.Due to the effects of an earthquake,the whole domain is considered,which should be large enough to eliminate the boundary effect on the solution,as shown in Fig.2a-c.The typical finite element mesh and displacement conditions for the upper bound limit analysis of shallow strip foundations using NS-FEMare illustrated in Fig.2d.The soil is described as a Mohr-Coulomb material with the cohesionc,unit weightγand surcharge loadq.The seismic acceleration coef ficient is denoted askhandkvin the horizontal and vertical directions,respectively,based on the pseudo-static approach.The horizontal displacements are free(u≠0)or fixed(u=0)along the ground surface to describe the smooth or rough interface condition between the soil and the strip footings.

Sweet-peas drooped3 over the boxes, and the rose-bushes shot forth4 long branches, which were trained round the windows and clustered together almost like a triumphal arch of leaves and flowers

Without the soil weight,the seismic bearing capacity factoris analogous to the static factorNcclosed-form solution reported by Prandtl(1920).Later on,Reissner(1924)established the analytical formulation for the static bearing capacityNqbased on the Prandtl failure mechanism.

Fig.2.The geometry and soil properties of the problem due to soil inertia.

Fig.3.Comparisons of plastic dissipation distribution with the characteristics line networks for the N s qE problem for?=30°,and k v=0:(a)k h=0,(b)k h=0.1,(c)k h=0.3,and(d)k h=0.5.

Fig.3 shows a comparison of plastic dissipation distribution using NS-FEM with that obtained by the characteristics line networks reported by Cascone and Casablanca(2016)in the case of ?=30°,kv=0,andkh=0-0.5.The power dissipations based on the NS-FEM are almost identical with those using the slip-line method from Cascone and Casablanca(2016).Under static conditions,the failure mechanism is symmetrical.Under seismic conditions,the failure pattern becomes asymmetrical.For increasingkh,the larger inertia force in the soil mass must be balanced,and the failure zone becomes longer.

To show the ef ficiency of the seismic bearing capacity factors using NS-FEM,we consider the computational cost based on the number of variables and CPU times for the terms(smooth interface)in the case of?=30°,kh=0.3,andkv=0.The reported CPU times refer to the time spent on the interior-point iterations for solving the seismic bearing capacity factors.The factorsthe number of variablesNvarand the CPU time obtained by the finite element analysis using triangular elements(FEM-T3),the edge-based smoothed FEM(ES-FEM-T3),and the NSFEM using triangular elements(NS-FEM-T3)are summarized in Tables 1 and 2.

The convergence rate archived by the present method(NS-FEM)is compared with those obtained by FEM-T3 and ES-FEM-T3,as shown in Fig.4.Although the coarse mesh is used,the seismic bearing capacity factorsusing the NS-FEM are more convergent than other existing methods such as FEM and ES-FEM using triangle elements.In addition,with the same number of elements,the total number of NS-FEMvariables is smaller than those of FEM-T3 and ES-FEM-T3.Regarding the maximum number of variables(Nvar=40,046),according to the CPU time for solving the optimization problem using NS-FEM based on the interior-point algorithm,very fast convergence is attained using the NS-FEM limit analysis(i.e.the CPU time of 5.8 s).These comparisons prove the effectiveness and rapid convergence of the NS-FEMwhen using the MOSEK optimizer to solve large sparse SOCP problems.

We conducted several simulations for the case of?=0°,in which the Mohr-Coulomb condition reduces to the unbonded Tresca yield criterion for plane strain problems,as noted in Tin-Loi and Ngo(2003),Vicente da Silva and Antao(2007)and Ciria et al.(2008).The numerical results are given in Table 3 for the cases ofkh=0-0.5andkv=0.In this case,the stable convergence of solutions and the accurate values of collapse loads as shown in Fig.5 prove that the numerical approach is naturally“immune”from the volumetric locking.Readers are referred to Meng et al.(2020)who stated that the application of NS-FEMand SOCP to elastoplastic analysis of static problems is free from the volumetric locking while using the loworder mixed element.In that study,the well-known“overly stiff”phenomenon was addressed by varying the Poisson’s ratio from 0.4 to 0.49999,while the solution still becomes stable and accurate.In addition,the special arrangement of linear elements(Fig.2),in which four elements forma rectangle domain with the central node in the intersection of the diagonals,can eliminate the locking phenomenon as noted in Vicente da Silva and Antao(2007).Readers can refer to Sloan and Kleeman(1995)who reported that the volume locking can be removed by increasing the ratio between the number of compatibility conditions and the number of degrees of freedom used in the discretization in the optimization problem.

Table 1Comparisons of seismic bearing capacity factor N s qE using the NS-FEM and other solutions(?=30°,k h=0.3 and k v=0,smooth interface).

Table 2Comparisons of seismic bearing capacity factorusing the NS-FEM and other solutions(?=30°,k h=0.3,k v=0,smooth interface).

Table 2Comparisons of seismic bearing capacity factorusing the NS-FEM and other solutions(?=30°,k h=0.3,k v=0,smooth interface).

?

Table 3The seismic bearing capacity factor NcE for the case of?=0°,k h=0-0.5,and k v=0.

Table 4Seismic bearing capacity factorof strip footing(k v=0).

Table 4Seismic bearing capacity factorof strip footing(k v=0).

?

Table 5Seismic bearing capacity factorof strip footing(k v=0).

Table 5Seismic bearing capacity factorof strip footing(k v=0).

?

Table 6Seismic bearing capacity factorsandof strip footing(k v=0).

Table 6Seismic bearing capacity factorsandof strip footing(k v=0).

?

Table 7Seismic bearing capacity factorof strip footing(k v=0).

Table 7Seismic bearing capacity factorof strip footing(k v=0).

?

Table 8Variations in seismic bearing capacity factors with various combinations of k h and k v forφ=30°.

Table 8(continued)

Table 9Comparisons of the numerical results of NqE with solutions available in literature(φ=30°and k v=0).

Table 10Comparisons of the numerical results of NγE with solutions available in literature(?=30°and k v=0).

Fig.4.The convergence rate of seismic bearing capacity factors in the case of?=30°,k h=0.3 and k v=0:(a)and(b)(smooth interface).

Fig.5.Calculation of(a)static and(b)seismic bearing capacity factors for?=0°,k h=0-0.5 and k v=0.

3.2.Effect of superstructure inertia on seismic bearing capacity factors

The ultimate load of strip footings under seismic conditions that only consider superstructure inertia is expressed as the following formulation:

Fig.6.Seismic bearing capacity factor N s qE:(a)Effect of horizontal acceleration;and(b)Effect of vertical acceleration for the case of?=30°.

Fig.9 illustrates the geometry and soil parameters to investigate the effect of superstructure inertia on the bearing capacity factors of strip footings.Due to the effects of seismic force,the rectangular domain is considered large enough to eliminate the boundary effect on the solution,as shown in Fig.9a-c.The typical finite element mesh and displacement conditions for the upper bound limit analysis of shallow strip foundations using NS-FEM are illustrated in Fig.9d.

Figs.10-12 show the power dissipation and displacement fields for strip footing problems using NS-FEM compared with the characteristics line networks reported by Cascone and Casablanca(2016)for?=30°,kv=0,andkh=0-0.5.Under static conditions(kh=0),the failure mechanism of shallow strip foundations is symmetrical,as shown in Figs.10a,11a and 12a.Under seismic conditions(kh＞0),the failure mechanisms become asymmetrical and identical to that of the characteristic line networks using the slip-line method from Cascone and Casablanca(2016).With an increasing value ofkh,the seismic beaering capacity factors reduces signi ficantly due to inclination of the applied load acting on the shallow foundation,leading to the observation that the length and the depth of failure zones become shallower than those for the case ofkh=0,as shown in Fig.10b-d,11b-d,and 12b-d.

Fig.7.Seismic bearing capacity factor N sγE:(a)Effect of horizontal acceleration for the smooth foundation;(b)Effect of horizontal acceleration for the rough foundation;(c)Effect of vertical acceleration for the smooth foundation(?=30°);and(d)Effect of vertical acceleration for the rough foundation(?=30°).

Fig.8.Effect of k v on the corrective coef ficients for the case of?=30°:(a)and(b)for smooth foundation;and(c)for rough foundation.

Fig.9.The geometry and soil properties of the problem due to superstructure inertia.

Fig.10.Comparisons of plastic dissipation distribution with the characteristics line networks for theproblem for?=30°,k v=0:(a)k h=0,(b)k h=0.1,(c)k h=0.3,and(d)k h=0.5.

Fig.11.Comparisons of displacement fields with the Hill’s failure mechanism for theproblem with smooth interface for?=30°and k v=0:(a)k h=0,(b)k h=0.1,(c)k h=0.3,and(d)k h=0.5.

Fig.12.Comparisons of displacement fields with the Prandtl’s failure mechanism for theproblem with smooth interface for?=30°and k v=0:(a)k h=0,(b)k h=0.1,(c)k h=0.3,and(d)k h=0.5.

3.3.The link between the static and seismic bearing capacity factors

For practical application,the seismic bearing capacity of strip footings can be evaluated by the combination of the effect of soil and superstructure inertia as follows:

Fig.13.Effect of superstructure horizontal acceleration k h on seismic bearing capacity factors:(a)for smooth foundation,and(d)for rough foundation.

Fig.14.Effect of superstructure horizontal acceleration k h on the corrective coef ficients:(a)for smooth footing,and(d)for rough footing.

Fig.15.The effect of soil and superstructure inertia on corrective coef ficients in the case of?=30°,k v=0,and k h=0-0.5:(a)eqE,(b)eγE for smooth footing,and(c)eγE for rough footing.1-Effect of soil inertia;2-Effect of superstructure inertia;3-Effect of soil and superstructure inertia.

Fig.16.Comparison of the present results with previous studies in the case of?=30°,k v=0,and k h=0-0.5:(a)NqE,(b)eqE,(c)NγE for rough footing,and(d)eγE for rough footing.

4.Conclusions

Using a pseudo-static approach,a NS-FEM was presented to estimate the static and seismic bearing capacities of shallow strip footings.A large number of simulations have been performed to study the bearing capacity of shallow strip footing under both static and seismic conditions,providing a new upper bound solution to all seismic bearing capacity components.Numerical results of both static and seismic bearing capacity factors are in good agreement with those using the slip-line method which correspond to the lower bound solutions to the bearing capacity.In addition,power dissipations obtained are identical to the slip-line networks given by Cascone and Casablanca(2016),providing knowledge of seismic effects on the failure mechanism of strip footing.Based on the numerical analyses,the following conclusions are made:

(1)In terms of the numerical procedure,combining the NS-FEM into the upper bound procedure can give accurate and reliable values of both static and seismic bearing capacity factors of shallow strip footing.The inclusion of horizontal and vertical seismic accelerations is presented in a straightforward and simple manner,providing an effective limit analysis to capture the ultimate load and the corresponding failure pattern under both static and seismic conditions.

(2)Separate calculations of all components of static and seismic bearing capacity factors to consider the effect of soil and superstructure inertia are performed and checked against the well-known Terzaghi’s formula,with very satisfactory results.

(3)The numerical results reveal that the seismic bearing capacity depends on the direction of seismic acceleration in a highly complex manner.An upward vertical seismic acceleration results in a decrease in the magnitude of seismic bearing capacity.On the other hand,the vertical seismic acceleration in the downward direction causes a rise in the bearing capacity

(4)The corrective coef ficients obtained show that superstructure inertia is a major contribution to reductions in the bearing capacity under seismic conditions compared with soil inertia.To facility practical applications,the superposition of both the soil inertia and the structure inertia has been presented to give safety values of seismic termsNcE,NqEandNγE.

(5)Variations of seismic bearing capacity with various combinations of seismic accelerationskhandkvare given in design tables and charts,providing useful knowledge of dependent seismic bearing capacity.These values can be used as references for engineers in the seismic design of a shallow foundation.

Declarationofcompetinginterest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to in fluence the work reported in this paper.

Acknowledgments

This is part of the TPS project.The first author is grateful to a Vied-Newton PhD scholarship and a Dixon scholarship from Imperial College London,UK,for supporting his studies at Imperial College London.He is also indebted to the Dean’s Fund from Imperial College London for financial support(2017-2020).