Bearing capacity of embedded foundations using quasi-kinematic limit analysis

Computers and Geotechnics 117 (2020) 103275

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Bearing capacity of embedded foundations using quasi-kinematic limit analysis Debasis Mohapatra, Jyant Kumar

T

⁎

Department of Civil Engineering, Indian Institute of Science, Bengaluru 560012, India

ARTICLE INFO

ABSTRACT

Keywords: Limit analysis Failure Foundations Plasticity Semidefinite programming

Bearing capacity of rectangular foundations, with varying aspect ratios (1 ≤ L/B ≤ 5) and embedded at shallow to medium depths (D/B ≤ 5), has been evaluated by using the three-dimensional cell-based smoothed finite elements quasi-kinematic limit analysis. The results have also been obtained for circular and strip foundations by performing exclusively the axisymmetric and plane strain analyses. The associated optimization problem has been solved by employing the semidefinite programming (SDP). Shape factors (sc, sq and sγ) and depth factors (dc, dq and dγ) are given as a function of soil internal friction angle (ϕ). The shape factors have been expressed as a function of L/B, on the other hand, the depth factors have been given as a function of D/B and L/B. A thorough comparison of the results has been made with the different solutions available in literature. The variations of (i) power dissipation function, (ii) maximum plastic shear strain rates, and (iii) the nodal velocities patterns, have also been examined to interpret the failure mechanism.

1. Introduction The ultimate bearing pressure ( pu ) for a footing embedded at a depth D below ground surface is often evaluated by means of the following equation:

pu = P u / A = cNc sc dc + qNq sq dq +

1 BN s d 2

(1)

where (i) Pu is the ultimate vertical load for the foundation having an area A; (ii) Nc , Nq and N refer to the bearing capacity factors for a rough strip footing placed at ground surface (D/B = 0); (iii) sc , sq and s define the shape factors; (iv) dc , dq and d correspond to the depth factors; (v) q is the uniform external surcharge pressure exerted at ground surface (not simply the soil overburden pressure at the footing base level as usually defined in literature [1–3]); and (vi) and c refer to unit weight and cohesion of soil mass, respectively. The effect of soil overburden pressure and the shear strength due to overlying soil mass above the footing level has been duly accounted in the depth factor d as defined by Eq. (1). The exact values of the bearing capacity factors Nc and Nq can be obtained analytically from the plane strain analysis by using the theory of plasticity [4,5]. Although the exact analytical solutions are not feasible for the bearing capacity factor N , very accurate numerical solutions for the factor N of a strip footing are available [6–10]. By means of experimental and simplified theoretical analysis, shape and depth factors have been established by several researchers. ⁎

The formulae for the shape factors proposed by different researchers have been summarized in Table 2 of [11]. The formulae for depth factors proposed by the different researchers are given in Table 1. The shape and depth factors proposed by Meyerhof [1], Hansen [2] and Vesic [3] are based on small scale model as well as prototype tests’ results in combination with the analytical method. The numerical solutions for the shape factors are available on the basis of (i) the upper bound limit analysis approach based on the multi-blocks rupture mechanism [12], (ii) the displacements based elasto-plastic finite element method [13], (iii) the finite difference method by using FLAC [14], (iv) lower and upper bound based finite element limit analysis approach [15–18]. The authors have recently carried out a three-dimensional (3D) kinematic limit analysis, based on the radial point interpolation method (RPIM), to compute shape factors of circular and rectangular footings on cohesive-frictional soil mass [11]. Depth factors for strip and circular footings for cohesive soils under undrained condition ( = 0 ) have been established numerically by different researchers [15,19–21]. A comprehensive numerical analysis for determining the bearing capacity of a rectangular footing embedded in sand at different depths has been performed by Lyamin et al. [17] with the usage of the rigorous lower and upper bound 3D finite element limit analysis in conjunction with the nonlinear optimization. However, the differences between the lower and upper bound results from this numerical approach were found to be quite substantial, and in order to predict the bearing

Corresponding author. E-mail addresses: [email protected] (D. Mohapatra), [email protected] (J. Kumar).

https://doi.org/10.1016/j.compgeo.2019.103275 Received 22 February 2019; Received in revised form 20 September 2019; Accepted 20 September 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

2. Problem definition

Table 1 Expressions for the depth factors proposed by the different researchers. Proposal

Depth Factors

Meyerhof [1]

dc = 1 + 0.2 kp

D B

any

dq = d = 1 + 0.1 kp Hansen [2]

dq = d = 1

=

dq = 1 + 2tan (1 Vesic [3]

> 10°

sin

dc = 1 + 0.4k any d = 1 any dc = 1 + 0.27

D B

k p = tan2 (45 + /2) ; k = D / B for D / B radians. Nq = k p e tan ; Nc = (Nq 1)cot .

Nq k Nc

> 0°

sin ) 2 k any = 0°

dq = 1 + (0.00440 + 0.356) d = 1 any

sin )2

any

)2 k

d = 1 any

dq = 1 + 2tan (1 Salgado et al. [15] and Lyamin et al. [17]

D B

0° = 0°dc = 1 + 2(1

dc = 1 + 0.4k

The present study deals with the determination of the bearing capacity of a rigid rough rectangular footing of width B and length L which is subjected to vertical load without any eccentricity. The footing is embedded in a homogeneous cohesive-frictional soil medium at a depth D below the ground surface. The yield in soil mass is governed by means of the Mohr-Coulomb criterion and the associated flow rule. The soil domain has volume V bounded by surface area A such that A = Au At and Au At = ; (i) along Au the velocities are prescribed, and (ii) along At the tractions are defined. The soil domain is acted upon by known body forces per unit volume (g) and known external surface tractions (t). In the present case, (i) the body forces exist on account of the unit weight of soil mass (γ), and (ii) the surface tractions are introduced by imposing uniform surcharge pressure (q) on ground surface. To avoid complications on account of the presence of the footing base and the associated vertical shaft (column), the footing is simply modelled as a zero-thickness rigid rough plate without any presence of the shaft. This modelling strategy is geometrically simple and produces a result that is close to the desired quantity (pure unit base resistance) with only a slight conservative bias when compared with other possible modelling options [17]. It should be mentioned that the soil overburden pressure, including its shear strength, above the footing level has been duly taken into consideration in the analysis. Due to symmetry, one-quarter of the total problem domain only needs to be modelled and corresponding governing boundary conditions are indicated in Fig. 1(b). The footing is subjected to uniform incremental vertical velocities (V0 ) at all the nodes located along the base of the rigid footing, while horizontal velocities are specified to be zero at these nodes. On account of the presence of the singularity around the corner of footing, a very fine mesh is generated around the footing edge as shown in Fig. 1(c). Along the chosen distant boundaries, away from the footing, both horizontal and vertical velocities were kept equal to zero. The analysis showed that these boundary conditions hardly affect the resulting collapse load as long as the boundaries remain far distant from the yielded zones. The distances of the chosen bottom horizontal and lateral vertical boundaries from the footing vary depending on different shape and depth of footing in order to ensure that the chosen boundaries did not interfere with the development of the collapse mechanism while maintaining adequate mesh density close to the footing. Different sizes of the problem domain used for different combinations of L/B and D/B are given in Table 2.

( )

1; k = tan

D B

0.28

any

1 (D / B )

for D / B > 1, k in

capacity, an average of the corresponding lower and upper bound solution was recommended by Lyamin et al. [17]. It is observed that the 3D bearing capacity solutions obtained by the different researchers do not match closely with each other; this will be revealed in the later part of this manuscript. There is a need to improve the existing bearing capacity solutions wherever the 3D analysis is involved especially for the cases where the foundation is embedded at some depths from ground surface. In the present research, a 3D upper bound limit analysis, using the cell-based smoothed finite elements method (CS-FEM), has been performed to compute the bearing capacity of rectangular foundations embedded at shallow to medium depths. Unlike the RPIM method, used by the authors in their recent study [11], where the strain rate is being evaluated only at the centroid of an element and is assumed to be satisfied throughout the element, in the CS-FEM, the strain rate is computed as the spatial average over the element. On account of its better computational stability as well as the efficiency in terms of its execution time, the CS-FEM has been employed in the present research rather than the RPIM [22,23]. In recent years, the CS-FEM has been implemented by a number of researchers to solve different plane strain stability problems in geomechanics using the kinematic limit analysis on the basis of the von-Mises and the MohrCoulomb yield criteria [23–25]. In this approach, the finite element is subdivided into different cells, and the strain rates are smoothed over each cell. The CS-FEM based upper bound limit analysis approach often provides very accurate, locking free and computationally inexpensive solutions [23–25]; a locked solution implies the forced restrained movement (s) in a particular direction at some nodes. On account of its overall versatility, the CS-FEM based kinematic limit analysis has been employed in the present research by enforcing the Mohr-Coulomb yield criterion. The optimization was performed with the help of the semidefinite programming (SDP) technique by means of an efficient primaldual interior point algorithm [26]. The bearing capacities of rectangular, circular and strip footings, placed at shallow to medium depths, have been computed. Axisymmetric and plane-strain analyses have also been carried out in addition to obtain the solutions for circular and strip foundations embedded at shallow to medium depths. The results are presented in terms of the shape and depth factors by computing the bearing capacities for rectangular footings with several combinations of L/B and D/B. The nodal velocity patterns, the distribution of power dissipation function and maximum plastic shear strain rates have also been plotted to interpret the collapse mechanisms in different cases

3. The CS-FEM based kinematic limit analysis The kinematic limit analysis theorem states that for a perfectly plastic material, following an associated flow rule with a kinematically admissible velocity field, the collapse load determined by equating the total internal plastic power dissipation to the rate of total external work done gives an upper bound of the actual collapse load. By using the expressions of the plastic power dissipation function and the associated flow rule constraints for the 3D Mohr-Coulomb yield criterion, as given in [27,28], the associated optimization problem for the determination of the magnitude of the collapse load (Pu) can be expressed as [11]:

Minimize

=

+

2c cos V 1 + sin

tr(E+)

dV

At

V0

=

u

in

V

tTu dAt

V

gTu dV (2a) (2b) (2c)

1 + sin tr(E+) = 1 sin tr(E )

u = u 0 on Au where (i) c and 2

in V

(2d) (2e)

refer to cohesion and internal friction angle of soil

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 1. (a) 8-noded hexahedral element divided into 4 smoothing cells [30]; (b) velocity boundary conditions; (c) discretization of domain for L/B = 1 and its zoomed view near footing.

mass, (ii) u and denote the velocity and plastic strain rate vector, T respectively; that is, u = [ux u y ]T and = [ xx yy zz xy yz zx ] , (iii) and u are related by a linear differential operator , (iv) u 0 is the prescribed velocity vector on Au , (v) To solve the problem by using the SDP, as earlier defined by Eq. (2b), is decomposed into (i) a dilatant part (vector),

+

=

+ xx

+ yy

+ zz

+ xy

+ yz

+ T , zx

The internal power dissipation function due to the development of the plastic strain (the first term in Eq. (2a)) is expressed as the product of the dilatant volumetric strain rate and the tensile strength; likewise, it can also be expressed as the product of the contractant volumetric strain rate and the compressive strength providing exactly the same answer [27]. It should be mentioned that + , and u form the basic variables in the present optimization problem framed by Eqs. (2a)–(2e). For an axisymmetric case, the previous work of the authors [29] can be referred. Eq. (2) is converted to a numerical optimization problem by using the CS-FEM based discretization technique. The key idea of the CS-FEM based kinematic limit analysis is to combine the standard finite element method discretization with a strain smoothing scheme [23,30]. The 3D problem domain ( ) is first discretized into a number of 8-noded hex1 2 3 nel ahedral finite elements such that and i j = i j . Elements are then subdivided into different smoothing cells (SC) (as shown in Fig. 1(a)) and the smoothing strain operations are performed along each SC. The application of a strain smoothing technique to the standard displacement finite element resolves the volumetric locking problem and also results in an efficient method which can provide accurate solutions with the minimal computational effort [23,25,30]. It should be noted that the cell-based smoothed finite element method, with just one smoothing cell per element (CS-FEM1), offers a very good combination of accuracy as well as computational efficiency [23,25,30], and the same was, therefore,

and (ii) a contractant T

= [ xx yy zz xy yz zx ] . Unlike the direct part (vector), and linear relationship between and u using the operator , there are no direct expressions relating either (i) + and u or (ii) and u . Further, E is a symmetric plastic strainratetensor, which is similarlydecomposed + into a dilatant part (E ) and a contractant part (E ), that is, E = E+ E . The components of E , E+ and E are given as

E=

xx

xy

xz

yx

yy

yz

zx

zy

zz

+

E =

E =

+ xx + yx + zx

; + xy + yy + zy

+ xz + yz + zz

xx

xy

yx

yy

yz

zx

zy

zz

;

and

xz

;

3

Computers and Geotechnics 117 (2020) 103275

17.27B × 12.28B 17.27B × 12.28B 17.27B × 12.28B 17.27B × 12.28B 19.38B × 12.67B 22.32B × 14.50B 25.25B×14.50B 28.99B × 14.50B 115,596 115,596 121,524 147,264 166,144 206,682 216,384 280,623

11,186 11,186 11,186 11,186 12,126 14,061 15,151 16,241

implemented in the present analysis. One can employ greater number of smoothing cells per element, however, the computational time to obtain the solution will obviously further increase. The velocity field over the eight noded brick element is varied by using the standard shape functions. A constant smoothed strain rate field is generated by using gradient smoothing technique all over the ~ domain of SCs. The smoothed strain rate ( ) over the smoothing cell Ce is defined as follows: ~ ~ = (x c ) = = u (x) (x)d e (x) (x)d e (3) c c where (i) x c is any point in the smoothing cell Ce , and (ii) distribution (smoothing) function which is defined as:

115,596 115,596 121,524 147,264 166,144 206,682 216,384 260,680 82,485 82,485 86,715 121,524 149,742 188,160 188,160 207,360 12.41B × 12.41B × 10B 12.41B × 12.41B × 10B 12.10B × 12.41B × 10.78B 16.07B × 15.92B × 10.93B 19.78B × 18.36B × 14.61B 20.87B × 20.62B × 13.10B 20.87B × 20.62B × 13.10B 21.71B × 21.46B × 19.35B 12.41B × 12.41B × 10B 12.41B × 12.41B × 10B 12.41B × 12.41B × 10.78B 16.07B × 15.92B × 10.93B 19.78B × 18.36B × 14.61B 20.87B × 20.62B × 13.10B 20.87B × 20.62B × 13.10B 21.71B × 21.46B × 19.35B 82,485 82,485 86,715 88,830 100,674 156,114 156,114 168,948

(x)d

e c

(x) =

{

=1

(4a) e c e c

1/ Vc x 0 x

(x) is a

(4b)

Vc = e d : the volume of the smoothing domain ec . By using the c smoothing function given in Eq. (4) and the Green’s divergence the~ orem, the smoothed strain rate field ( ) can be calculated using the boundary integral as given below: ~ 1 1 (x c ) = u (x)d = ncT (x) u (x)dS c c Vc Vc (5) where Sc is the closed boundary of the smoothing domain c and n c (x) is the matrix of the normals associated with the boundary Sc and is defined as:

nx 0 0 n c (x) = ny 0 nz

0 ny 0 nx nz 0

0 0 nz 0 ny nx

(6)

Substituting Eq. (6) in Eq. (5), the smoothed strain rate field can be expressed as: ~ ~ (x c ) = B C d (7) where d = [u1 v1 w1 vector for each element.

u8 v8 w8 ]T refers to the nodal velocity

~ B

~ B

C

~ ~ ~ = BC 1 BC 2 B ~ Ni, x (x c) 0

~ BCi =

=

1 Vc

C3

0 ~ Ni, y (x c)

(8)

C8

0 0 ~ Ni, z (x c)

0 0 ~ ~ Ni, y (x c)/2 Ni, x (x c)/2 0 ~ ~ 0 Ni, z (x c)/2 Ni, y (x c)/2 ~ ~ Ni, z (x c)/2 0 Ni, x (x c)/2

~ 1 Ni, m (x c) = Vc

sc

Ni (x) nm (x)dS

(m = x , y , z )

(9) (10)

nb

Ni (x Gb) nmb Ab

(11) ~ where (i) Ni, m is the derivative of the smoothed shape function, (ii) nb refers to the number of boundary surfaces associated with the c , and (iii) x Gb is the Gauss point of the boundary surface Scb which has area Ab and the outward normal nmb . From Eq. (11), it is evident that the smoothed strain rate is computed based on the integration of the

0.0 0.2 0.4 0.6 1.0 2.0 3.0 5.0

12.41B × 12.41B × 10B 12.41B × 12.41B × 10B 12.41B × 12.41B × 10.78B 12.41B × 12.41B × 10.93B 14.05B × 14.05B × 10.65B 18.16B × 18.16B × 13.55B 18.16B × 18.16B × 13.55B 18.16B × 18.16B × 15.25B

82,485 82,485 86,715 121,524 149,742 188,160 188,160 207,360

16.07B × 15.92B × 10B 16.07B × 15.92B × 10B 16.07B × 15.92B × 10.77B 19.44B × 19.75B × 10.20B 19.44B × 19.75B × 12.09B 20.90B × 20.31B × 15.41B 24.29B × 24.29B × 14.75B 27.36B × 26.22B × 14.74B

16.07B × 15.92B × 10B 16.07B × 15.92B × 10B 16.07B × 15.92B × 10.77B 19.44B × 19.75B × 10.20B 19.44B × 19.75B × 12.09B 20.90B × 20.31B × 15.41B 24.29B × 24.29B × 14.75B 28.43B × 28.42B × 14.75B

Domain Nele Domain Nele Domain Nele Domain Nele Domain Nele Domain

2.0 1.5 L/B = 1 D/B

Table 2 The sizes of the problem domain and the number of elements for different combinations of L/B and D/B.

3.0

5.0

∞ (Strip Footing)

Nele

D. Mohapatra and J. Kumar

4

b=1

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

boundary segment Scb Sc and no shape function derivatives are required. The values of the shape function at the Gauss quadrature points of the boundary of the smoothed domain are only required for doing the associated numerical integration. As the shape functions remain linear in nature over the boundaries, one-point Gauss quadrature is sufficient for doing the numerical integration over the boundary segment Scb . The final optimization problem in Eq. (2) can be numerically expressed as: ~+ nc nt nc V k tr(E (xi)) A t (xi)TNi di V g (xi )T Ni di i = 1 ci i i = 1 ti i = 1 ci Minimize V0

Table 3 A comparison of the bearing capacity factor Nc for ϕ = 30° by using different limit analysis approaches.

384 1536 9600 21,600 38,400

(12a)

i

=

(m)iT

~+

~

i

~+ i

i

~ = BCi di

+ (n)iT

~ i

=0

i = 1,

nc

(12c)

i = 1,

nc

(12d) (12e)

u = u 0 on Au

where (i) nc is total number of smoothing cells, (ii) Vci implies the volume of the smoothing cell i, (iii) Ati is area of the face i on the traction boundary, and (iv) N is the matrix of shape functions. The parameters m, n, b and k are given by:

mT = [ b

b

b];

nT = [1 1 1]; (1 sin ) b= ; and (1 + sin ) 2c cos k= 1 + sin

FELA-T6

FELA-Q4

Analytical

Nc

CPU (s)

Nc

CPU (s)

Nc

CPU (s)

32.675 31.632 30.786 30.581 30.475

0.4 0.9 5.6 15.7 30.7

42.805 35.785 32.508 31.633 31.292

0.5 1.4 9.1 29 52

32.675 31.632 30.786 30.581 30.475

0.5 0.9 5.9 16.2 35.9

30.140

elements, (iii) FELA-Q4 - the conventional finite elements limit analysis by employing standard four-noded quadrilateral element using reduced integration (1 Gauss point per element). A domain of 8B × 4B was discretized uniformly by using 4 noded quadrilateral elements and 6 noded triangular elements. The analytical solution of the bearing capacity factor Nc for this problem is 30.140 [4]. For different number of chosen elements, the obtained numerical solutions by using these three different approaches are reported in Table 3. From this comparison of the results, it is observed that in all cases the CS-FEM1 provides more accurate solutions than the FELA-T6. On the other hand, the results by using the CS-FEM1 and FELA-Q4 are found to be exactly the same in all the cases. In terms of its computational efficiency, the CS-FEM1 requires minimum computational time as compared to the other two discretization techniques. The difference between the computational times associated with the CS-FEM1 and the FELA-Q4 is found to be very marginal. However, when elements with the distorted shapes are used, as concluded by [23], the FELA-Q4 fails to compute the solution in many cases due to the problems associated with in finding the inverse of the Jacobian matrix due to its badly ill conditioned nature. On the other hand, the CS-FEM1 works well in all the cases. In the present analysis, square shaped quadrilateral elements were used, therefore, the CS-FEM1 and the FELA-Q4 solutions are found to match exactly with each other. Unlike the FELA-Q4, the CS-FEM1 does not require finding either the derivative of the shape function or the inverse of the Jacobian matrix. In the FELA-Q4, unlike the CS-FEM1, the flow rule is not guaranteed to be satisfied everywhere in the domain since the strain rates in the FELA-Q4 are not constant [23].

(12b)

~

CSFEM1

Nele

(13)

Eq. (12) comprises of a linear objective function since V0 is known constant footing velocity subjected to dual SDP and linear constraints which can be solved by using the primal-dual interior-point method solver MOSEK [26]. For the sake of convenience and efficiency, the input to MOSEK is specified in its dual form as explained in [11]. The dual variables of the solution comprise of the basic unknown plastic strain rates and the nodal velocity vectors. All the simulations were performed on a Windows 10 based desktop system having Intel core I77700 CPU @ 4.20 GHZ and 16 GB RAM. The CPU time required for the interior point optimization for the strip footing is found to be approximately 1 min and that for a square footing, with D/B = 0, is found to be 10 mins. The maximum CPU time required to solve the problem of the largest problem domain, with L/B = 5 and D/B = 5, is found to be approximately 7 hrs. The CPU time required for doing the analysis is found to increase continuously with an increase in the number of the elements. It should be mentioned here that the enforcements of the boundary conditions in the CS-FEM and the conventional FELA are identical. In the RPIM, the enforcement of the boundary condition is, however, not fully guaranteed due an involvement of the chosen basis functions [11]. In the CS-FEM, by using the smoothed strain rates, which will be constant over a cell, the flow rule is guaranteed to be satisfied everywhere in the domain as indicated by Le et al. [23]. The approximate strain rates, using the strain smoothing defined in Eq. (7), relax the compatibility condition. Hence, it will not be possible to guarantee that the solution obtained by using Eq. (12) will be a strict upper bound. In order to make a comparison of the results obtained from the CSFEM with the other discretization approaches, in terms of the accuracy and the computational efficiency, a strip footing problem for a weightless cohesive-frictional medium with =30°, was analyzed by using the kinematic limit analysis. Three different discretization approaches, namely, (i) CS-FEM1 - four noded quadrilateral element with one smoothing cell per element, (ii) FELA-T6 – the conventional finite elements limit analysis using six noded linear strain triangular

4. Results The average ultimate bearing pressure ( pu ) for a footing base area (A) is expressed as:

pu =

Pu = cNc' + qNq' + 0.5 BN ' (L × B )

(14a) (14b)

= cNc sc dc + qNq sq dq + 0.5 BN s d

and N for a rectangular footing The bearing capacity factors become a function of the aspect ratio L/B, embedment ratio D/B and soil internal friction angle ; for circular and strip foundation Nc' , Nq' and N ' depend on D/B and . The equivalent bearing capacity ex-

Nc' ,

Nq'

'

pression associated with the surcharge (q'=q + D) at the footing base will be given by:

pu = cNc sc dc + q'Nq sq dq' + 0.5 BN s d

(14c)

where dq' = dq/(1 + D / q) . For a special case, when D / q=0, dq' will become simply equal to dq . These factors have been computed for different values of L/B ranging from 1 to 5, D/B varying from 0 to 5 and, the value of ranges from 0 to 45°. The results in terms of Nc' , Nq' and N ' for different combinations of L/B, D/B and have been presented in Tables 4–6. In these Tables, the results exclusively from the plane strain (L/B = ∞) and axisymmetric analyses (circular footing) by using four noded quadrilateral elements, CS-FEM and SDP have also been 5

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Table 4 The values of Nc' for different combinations of L/B, D/B and . = 0°

= 10°

= 20°

= 30°

= 40°

= 45°

D/B 0.0

Circular Footinga 6.167

0.2 0.4 0.6 1.0 2.0 3.0 5.0

8.393 10.023 11.392 13.334 13.336 13.336 13.336

D/B 0.0

Circular Footing 11.265

0.2 0.4 0.6 1.0 2.0 3.0 5.0

15.064 17.970 20.319 23.728 28.079 32.136 39.654

D/B 0.0

Circular Footing 24.015

0.2 0.4 0.6 1.0 2.0 3.0 5.0

31.348 37.348 42.302 49.887 65.650 81.298 112.888

D/B 0.0

Circular Footing 63.581

0.2 0.4 0.6 1.0 2.0 3.0 5.0

80.213 94.725 106.885 127.021 177.596 229.024 336.162

D/B 0.0

Circular Footing 231.678

0.2 0.4 0.6 1.0 2.0 3.0 5.0

281.635 324.293 361.804 431.657 603.857 784.504 1165.430

D/B 0.0

Circular Footing 528.386

0.2 0.4 0.6 1.0 2.0 3.0 5.0

628.348 710.933 787.701 930.005 1282.090 1650.719 2427.706

b

L/B = 1.0 5.908 (6.124) 8.159 9.827 11.092 12.685 12.74 12.771 12.788

1.5 5.763 (5.892) 7.633 9.050 10.292 12.169 12.521 12.537 12.555

2.0 5.653 (5.772) 7.302 8.587 9.688 11.573 12.331 12.342 12.369

3.0 5.524 (5.625) 6.941 8.064 9.035 10.744 12.11 12.127 12.132

5.0 5.413 (5.496) 6.639 7.621 8.483 10.018 11.911 11.933 11.937

5.170 (5.156) 6.123 6.890 7.578 8.858 11.541 11.543 11.543

L/B = 1.0 10.672 (11.004) 14.561 17.497 19.789 22.784 26.718 30.484 37.100

1.5 10.204 (10.398) 13.349 15.825 18.008 21.438 25.068 28.242 33.910

2.0 9.859 (10.039) 12.569 14.743 16.634 19.988 23.899 26.732 31.868

3.0 9.445 (9.595) 11.701 13.53 15.132 17.998 22.166 24.738 29.111

5.0 9.083 (9.201) 10.953 12.491 13.846 16.286 20.284 22.243 25.938

8.419 (8.377) 9.674 10.737 11.681 13.427 16.858 17.95 19.107

L/B = 1.0 22.656 (23.240) 30.065 35.959 40.922 47.820 62.151 76.542 104.202

1.5 21.046 (21.381) 27.015 31.999 36.528 43.899 56.218 68.115 90.753

2.0 19.862 (20.212) 24.867 29.105 32.927 39.881 51.899 62.374 82.451

3.0 18.464 (18.734) 22.475 25.891 28.967 34.605 45.228 53.436 71.424

5.0 17.209 (17.406) 20.381 23.102 25.553 30.024 38.458 45.104 58.256

14.964 (14.891) 16.742 18.291 19.68 22.245 26.926 29.595 32.886

L/B = 1.0 60.013 (61.264) 76.803 91.143 103.229 122.352 167.459 214.480 308.178

1.5 54.032 (54.468) 67.585 79.409 90.450 109.063 147.842 186.441 261.776

2.0 49.403 (50.216) 60.401 70.259 79.423 96.423 132.558 166.890 232.457

3.0 44.077 (44.649) 52.486 60.064 67.059 80.121 109.451 138.177 193.264

5.0 39.247 (39.666) 45.561 51.208 56.411 66.017 86.877 106.751 146.669

30.489 (30.288) 33.172 35.589 37.806 41.892 49.250 54.783 63.132

L/B = 1.0 220.182 (222.331) 270.669 314.221 350.881 413.049 570.666 739.155 1069.167

1.5 192.565 (192.283) 233.595 269.889 303.897 362.934 496.667 635.068 907.243

2.0 169.576 (171.849) 202.039 231.990 260.344 313.227 436.045 559.058 801.113

3.0 143.612 (145.089) 167.240 189.059 209.488 248.204 344.255 446.080 635.834

5.0 120.164 (121.295) 136.737 151.875 166.126 190.229 255.884 322.115 456.698

76.457 (75.779) 81.069 85.249 89.240 96.370 110.115 121.709 141.298

L/B = 1.0 504.679 (505.077) 606.907 695.082 765.875 897.227 1217.632 1582.796 2239.827

1.5 436.619 (440.007) 520.221 593.399 661.761 781.333 1058.225 1333.337 1886.675

2.0 376.818 (380.914) 442.585 504.044 561.116 667.576 919.075 1191.803 1646.739

3.0 309.765 (312.599) 356.362 399.625 440.382 517.447 712.303 916.969 1310.569

5.0 249.488 (251.686) 281.045 309.987 337.137 386.761 512.958 646.754 918.555

136.198 (135.060) 142.702 148.519 154.500 164.280 184.434 202.041 232.803

Note: The values within parenthesis corresponds to UB-LA based on RPIM results reported in Mohapatra and Kumar [11]. a By using present axisymmetric analysis. b By using present plane strain analysis. c By using UB-FELA using linear strain elements with continuous quadratic velocity field, SOCP and adaptive meshing.

6

Strip Footingc 5.167 6.126 6.904 7.546 8.801 11.513 11.524 11.569 Strip Footing 8.402 9.597 10.687 11.627 13.323 16.876 17.982 19.486 Strip Footing 14.987 16.682 18.226 19.726 22.23 27.172 30.044 33.351 Strip Footing 30.697 33.401 35.951 38.697 42.062 50.075 57.388 64.635 Strip Footing 78.947 82.936 86.714 90.862 101.469 115.201 125.711 146.598 Strip Footing 140.421 149.814 153.628 165.096 176.332 196.541 214.624 249.610

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Table 5 The values of Nq' for different combinations of L/B, D/B and . = 10°

= 20°

= 30°

= 40°

= 45°

D/B 0.0

Circular Footinga 2.984

0.2 0.4 0.6 1.0 2.0 3.0 5.0

3.656 4.168 4.582 5.182 5.948 6.663 7.977

D/B 0.0

Circular Footing 9.730

0.2 0.4 0.6 1.0 2.0 3.0 5.0

12.409 14.593 16.395 19.150 24.891 30.565 42.068

D/B 0.0

Circular Footing 37.662

0.2 0.4 0.6 1.0 2.0 3.0 5.0

47.381 55.694 62.695 74.781 103.523 133.224 195.083

D/B 0.0

Circular Footing 195.188

0.2 0.4 0.6 1.0 2.0 3.0 5.0

237.315 272.112 304.555 363.177 507.636 659.276 978.907

D/B 0.0

Circular Footing 528.894

0.2 0.4 0.6 1.0 2.0 3.0 5.0

629.346 711.867 788.755 930.694 1283.085 1651.669 2429.102

b

L/B = 1.0 2.882 (2.960) 3.567 4.085 4.486 5.018 5.812 6.375 7.534

1.5 2.799 (2.837) 3.354 3.791 4.176 4.780 5.420 5.971 6.98

2.0 2.738 (2.775) 3.216 3.600 3.933 4.525 5.214 5.714 6.614

3.0 2.666 (2.695) 3.065 3.386 3.668 4.179 4.911 5.357 6.143

5.0 2.602 (2.624) 2.923 3.203 3.441 3.871 4.577 4.922 5.573

2.484 (2.477) 2.706 2.893 3.059 3.367 3.972 4.165 4.369

L/B = 1.0 9.247 (9.517) 11.944 14.107 15.867 18.442 23.622 28.859 38.881

1.5 8.662 (8.799) 10.833 12.646 14.297 16.979 21.461 25.792 34.065

2.0 8.232 (8.375) 10.051 11.595 12.987 15.516 19.884 23.736 31.010

3.0 7.723 (7.831) 9.182 10.426 11.539 13.596 17.467 20.816 26.988

5.0 7.267 (7.342) 8.418 9.408 10.302 11.845 14.999 17.417 22.204

6.446 (6.420) 7.094 7.657 8.163 9.096 10.800 11.772 12.969

L/B = 1.0 35.616 (36.553) 45.341 53.628 60.426 71.503 97.685 124.832 178.931

1.5 32.199 (32.519) 40.026 46.849 53.224 63.982 86.358 108.643 152.137

2.0 29.529 (30.061) 35.874 41.566 46.858 56.673 77.513 97.354 135.064

3.0 26.448 (26.819) 31.312 35.682 39.702 47.240 64.216 81.053 112.583

5.0 23.659 (23.925) 27.305 30.568 33.569 39.116 51.160 62.635 85.683

18.603 (18.487) 20.152 21.547 22.827 25.186 29.435 32.629 37.449

L/B = 1.0 186.031 (188.258) 228.157 264.653 294.694 348.820 479.857 621.232 899.103

1.5 162.617 (162.661) 197.038 227.488 256.011 305.521 418.284 533.891 762.278

2.0 143.329 (145.483) 170.541 195.663 219.467 263.845 366.892 470.121 673.223

3.0 121.533 (122.914) 141.336 159.645 176.789 209.226 289.907 373.591 534.537

5.0 101.823 (102.876) 115.736 128.436 140.395 162.229 215.717 271.291 384.219

65.154 (64.586) 69.024 72.533 75.882 81.864 93.397 103.126 119.563

L/B = 1.0 505.930 (507.615) 607.854 696.071 766.807 898.167 1218.631 1583.791 2240.371

1.5 437.647 (441.805) 521.237 594.365 662.773 782.337 1059.226 1334.545 1887.520

2.0 377.843 (382.566) 443.605 505.109 562.093 668.324 920.077 1192.813 1647.663

3.0 310.836 (313.979) 357.479 400.689 441.483 518.544 713.315 917.973 1311.562

5.0 250.453 (252.906) 282.044 310.993 338.143 387.759 513.953 647.755 919.708

137.198 (136.060) 143.072 149.520 155.500 165.278 185.430 203.041 233.801

Strip Footingc 2.480 2.687 2.884 3.051 3.351 3.978 4.174 4.435 Strip Footing 6.455 7.071 7.633 8.180 9.091 10.889 11.935 13.138 Strip Footing 18.722 20.284 21.756 23.341 25.284 29.910 34.132 38.317 Strip Footing 67.244 70.591 73.761 77.242 86.142 97.665 106.484 124.010 Strip Footing 141.421 150.814 154.628 166.096 177.332 197.541 215.624 250.610

Note: The values within parenthesis corresponds to UB-LA based on RPIM results reported in Mohapatra and Kumar [11]. a By using present axisymmetric analysis. b By using present plane strain analysis. c By using UB-FELA using linear strain elements with continuous quadratic velocity field, SOCP and adaptive meshing.

presented. Furthermore, for a strip footing, the results were also separately obtained by using six noded triangular elements [31], adaptive meshes [32] and second order cone programming (SOCP). The results for rectangular surface footing reported by Mohapatra and Kumar [11] by using the RPIM based upper bound limit analysis also have been presented in these Tables. It can be noted that in most of the cases, the values of the bearing capacity factors obtained, from the CS-FEM with the usage of quadrilateral elements are found to very marginally smaller (reflecting indirectly a better upper bound solution) than that obtained with the employment of six noded triangular elements using

the adaptive meshes. The present results are found to be in good agreement with that reported by the authors [11] by using the RPIM for D/B = 0. Two sets of the results for strip and circular footings are nevertheless found to be very close to each other confirming the validity of the solutions established. Note that the bearing capacity factors for a circular footing on the basis of the axi-symmetric analysis are found to be generally slightly greater as compared to that established from the present 3D upper bound analysis for a square footing. This is in agreement with the observations drawn by Salgado et al. [15] and Lyamin et al. [17] The

7

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Table 6 The values of N ' for different combinations of L/B, D/B and . = 10°

= 20°

= 30°

= 40°

= 45°

D/B 0.0

Circular Footinga 0.363

0.2 0.4 0.6 1.0 2.0 3.0 5.0

1.721 3.087 4.686 8.818 18.197 29.348 56.122

D/B 0.0

Circular Footing 2.601

0.2 0.4 0.6 1.0 2.0 3.0 5.0

8.457 14.514 21.377 36.685 82.576 139.994 291.468

D/B 0.0

Circular Footing 16.430

0.2 0.4 0.6 1.0 2.0 3.0 5.0

40.711 65.833 93.587 154.812 351.782 608.121 1313.376

D/B 0.0

Circular Footing 127.299

0.2 0.4 0.6 1.0 2.0 3.0 5.0

258.097 390.337 531.263 844.194 1832.079 3132.875 6731.934

D/B 0.0

Circular Footing 426.344

0.2 0.4 0.6 1.0 2.0 3.0 5.0

788.919 1152.526 1531.326 2347.783 4900.273 8218.246 17305.850

L/B = 1.0

1.5

2.0

3.0

5.0

0.507 (0.492) 1.678 3.045 4.606 8.062 17.509 28.249 53.525

0.539 (0.524) 1.653 2.920 4.329 7.516 16.717 26.766 49.886

0.577 (0.542) 1.650 2.842 4.166 7.097 15.946 25.597 47.400

0.589 (0.559) 1.708 2.804 3.975 6.649 14.792 23.721 43.879

0.595 (0.574) 1.645 2.703 3.843 6.266 13.665 22.005 39.787

L/B = 1.0

1.5

2.0

3.0

5.0

3.092 (2.980) 8.183 14.269 21.053 35.987 80.119 135.016 276.387

3.163 (3.106) 7.967 13.307 19.363 32.952 74.613 124.616 248.296

3.184 (3.173) 7.803 12.787 18.177 30.334 68.607 115.472 229.016

3.239 (3.221) 7.653 12.132 16.957 27.557 60.484 101.023 200.449

3.299 (3.268) 7.493 11.588 15.901 25.165 53.578 87.636 167.959

L/B = 1.0

1.5

2.0

3.0

5.0

16.787 (17.212) 39.535 65.092 92.323 153.185 341.226 587.562 1242.399

17.079 (17.205) 37.274 59.139 83.415 137.897 309.606 525.907 1093.321

16.899 (17.129) 35.675 55.188 76.205 123.336 275.152 476.303 985.677

16.752 (16.870) 33.636 50.695 68.507 107.839 231.345 393.771 817.085

16.672 (16.669) 31.893 46.674 61.913 94.726 194.595 320.109 638.830

L/B = 1.0

1.5

2.0

3.0

5.0

126.953 (128.997) 254.697 389.379 530.155 835.394 1785.796 3025.891 6380.958

124.600 (121.943) 232.539 345.677 469.154 740.974 1595.282 2720.855 5564.783

118.320 (118.063) 213.631 311.057 413.696 641.568 1376.246 2359.300 4921.739

111.961 (110.516) 192.122 272.249 355.148 535.966 1102.017 1858.191 3902.047

103.898 (103.944) 172.797 238.449 305.423 446.077 875.641 1436.383 2864.481

L/B = 1.0

1.5

2.0

3.0

5.0

424.416 (424.964) 788.702 1159.999 1534.326 2343.970 4813.969 8111.321 16510.060

405.459 (392.569) 707.912 1014.025 1345.118 2059.362 4273.316 7214.366 14709.290

377.741 (372.481) 635.599 894.288 1164.158 1756.393 3641.817 6145.971 13010.800

345.118 (338.668) 554.123 760.938 973.259 1430.821 2851.182 4741.865 9799.381

311.638 (307.592) 482.035 644.903 809.741 1157.429 2193.794 3525.712 6988.569

b

0.461 (0.486) 1.591 2.551 3.519 5.607 11.681 18.844 32.969

Strip Footingc 0.457 1.578 2.524 3.487 5.546 11.679 18.916 33.389 Strip Footing

2.947 (2.965) 7.069 10.497 13.887 21.158 41.684 65.007 115.896

2.943 7.001 10.393 13.858 21.005 41.842 65.517 119.877 Strip Footing

15.361 (15.117) 28.359 39.381 49.892 72.244 132.572 199.893 351.030

15.597 28.277 38.969 49.778 72.590 137.098 205.801 363.495 Strip Footing

87.729 (86.839) 136.524 177.552 215.089 293.762 494.767 714.082 1203.498

88.391 138.569 185.390 227.003 311.987 520.720 770.726 1291.834 Strip Footing

239.989 (237.102) 346.758 435.378 516.003 681.674 1093.610 1534.179 2505.229

243.760 359.627 498.455 540.764 719.589 1207.579 1641.828 2763.583

Note: The values within parenthesis corresponds to UB-LA based on RPIM results reported in Mohapatra and Kumar [11]. a By using present axisymmetric analysis. b By using present plane strain analysis. c By using UB-FELA using linear strain elements with continuous quadratic velocity field, SOCP and adaptive meshing.

bearing capacity factors are found to increase with (i) decreases in value of L/B, (ii) increase in value of D/B, and (iii) increase in value of . By keeping D/B = 0, the shape factors, sc, sq and sγ due to the components of cohesion, surcharge and unit weight were determined by the dividing the bearing capacity factors Nc' , Nq' and N ' for a given L/B with the corresponding values of the bearing capacity factor (Nc, Nq and

Nγ) for a strip footing. On the other hand, for a given L/B, the values of the depth factors dc, dq and dγ were determined by dividing the values Nc' , Nq' and N ' for a given D/B with the corresponding value of the respective bearing capacity factor (Nc' , Nq' and N ' ) for D/B = 0. Note that for a given ϕ, the shape factors will remain only a function of L/B whereas the depth factors become functions of both D/B and L/B. The variation of the shape factors sc, sq and sγ with L/B for different

8

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

only up to D/B = 1–2 beyond which the value of dc becomes constant. The maximum value of dc increases very marginally with an increase in L/B. 3. For ϕ ≥ 10°, the values of all the depth factors (dc, dq and dγ) increase continuously with an increase in D/B. The values of the depth factors generally tend to become greater for smaller values of L/B. This is also in agreement with the findings from Lyamin et al. [17]. As compared to dc and dq, the depth factor dγ increases much more extensively with an increase in D/B. Note that for a given D/B, the difference in the values of the depth factor between square and strip footings increases continuously with an increase in ϕ. 5. Comparison of Results 5.1. Strip and circular foundations with D/B = 0 By considering the plastic equilibrium of continuum with the usage of the Mohr-Coulomb yield criterion, by satisfying the governing equations along two families of slip lines (stress characteristics), one can obtain very accurate solutions for determining the bearing capacity factors for strip and circular footings placed on ground surface (D/ B = 0). For a strip footing placed on ground surface, based on method of stress characteristics [5], a closed form solution is available for evaluating the bearing capacity factors Nc and Nq. On the other hand, even for D/B = 0, by using the slip line method, one needs to use a numerical procedure for obtaining Nc' and Nq' for a circular footing. Likewise, by using the slip line method, for determining the bearing capacity factor N ' for strip and circular footings, a numerical technique becomes always essential. Note that for the circular footing, while using the slip lines method, an assumption is always essential in order to define the hoop stress (σθ). For D/B = 0, Table 7 provides a comparison of the obtained values of the bearing capacity factors Nc' , Nq' and N ' for strip and circular footings with the available results of Prandtl [4] based on analytical method and Martin [9] based on the method of stress characteristics. It can be observed that the present analysis provides a very close match with the method of characteristics. In all the cases, due to an upper bound nature of the current approach, the bearing capacity factors from the present approach are noted to be very marginally greater than the corresponding values using the method of stress characteristics. This close comparison, therefore, confirms the accuracy of the present approach for plane strain and axi-symmetric cases. Fig. 2. The variation of sc, sq and sγ with L/B and ϕ.

5.2. Nc' for a square foundation for ϕ = 0 with D/B = 0 The obtained values of Nc' for a square footing, with ϕ = 0 and D/ B = 0, were compared with (i) the upper bound analysis of Michalowski [12] using rigid multi blocks mechanism, (ii) the lower and upper bound finite element analyses of Salgado et al. [15] using nonlinear programming (NLP), (iii) the elasto-plastic finite element analyses of Gourvencec et al. [33] and Nguyen and Merifield [34] using ABAQUS, and (iv) the upper bound limit analysis of Mohapatra and Kumar [11] using the RPIM. The comparison of these results is provided in Table 8. Note that the results from the present three-dimensional analysis compares very closely with the finite elements solutions of Gourvencec et al. [33] and Nguyen and Merifield [34]. The results from the present analysis lie in between the corresponding lower and upper bound results of Salgado et al. [15]. Note that the upper bound values of Nc' of Michalowski [12], Salgado et al. [15] and Mohapatra and Kumar [11] remain greater than the present results.

values of ϕ has been shown in Fig. 2. Likewise, the variation of depth factors dc, dq and dγ with changes in L/B and D/B for different values of ϕ has been presented in Figs. 3–5. The variation of depth factor dq' with different L/B, D/B and ϕ for γD/q = 10 and 20 has also been plotted in Fig. 6. Note that the factor dq' decreases continuously with an increase in the value of D / q . Following observations have been drawn with respect to the variation of the shape and depth factors: 1. The shape factors sc and sq become continuously smaller with an increase in L/B and decrease in ϕ. On the other hand, the trend of the variation of sγ with L/B and ϕ does not follow a very unique pattern. The similar observation has also been reported earlier in the references [11,13,14,17,18]. 2. For ϕ = 0, the depth factor dc increases with the embedment ratio

9

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 3. The variation of dc with L/B, D/B and ϕ.

5.3. Shape factors for rectangular foundation (D/B = 0)

upper and average solutions of Lyamin et al. [17] based on the threedimensional finite element limit analysis in combination with the nonlinear optimization for the variation of sγ, and (vi) The upper bound limit analysis results of Mohapatra and Kumar [11] obtained by using the RPIM. The comparison of the shape factors from the different analyses as a function of L/B for two different values of ϕ, namely, 30° and 40°, has been presented in Fig. 7(a)–(f). Following observations have been drawn from these comparisons: The shape factors sc and sq reported by Meyerhof [1] and Hansen [2] are found to be very conservative as compared to the present solution. The values of sc and sq reported by Michalowski [12], based on the 3D upper bound failure mechanism, are found be significantly greater. On

The values of the shape factor sc, sq and sγ computed from the present analysis were compared with the different existing results provided by (i) Meyerhof [1] and Hansen [2] by using the limit equilibrium method or slip line method in combinations with the experimental findings based on model tests, (ii) the upper bound analysis of Michalowski [12] using rigid blocks collapse mechanism for the variation of sc and sq, (iii) the elasto-plastic finite element analysis of Zhu and Michalowski [13] for the variation of sc and sq, (iv) the upper bound finite element limit analysis, using the continuous quadratic velocity field of Antão et al. [18] for the variation of sq and sγ, (v) the lower,

10

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 4. The variation of dq with L/B, D/B and ϕ.

the other hand, the values of sc and sq from the present analysis match well the solution of Zhu and Michalowski [13] based on the three-dimensional elasto-plastic finite element method. The values of the shape factors reported by Zhu and Michalowski [13] are found to be only marginally lower. The values of sq and sγ reported by Antão et al. [18], based on the upper bound finite element limit analysis, have been found to be very close to the present solution confirming the accuracy of the present analysis. The values of sγ from the present analysis lie in between the lower and upper curves of Lyamin et al. [17] based on the three-dimensional finite element limit analysis and the nonlinear optimization. Note that the gap between the lower and upper bound curves of Lyamin et al. [17] remains quite extensive and the present analysis compares reasonably well with the average of the corresponding lower

and upper bound results of Lyamin et al. [17]. The values of shape factors reported by Mohapatra and Kumar [11] based on the upper bound limit analysis using RPIM, have been found to match very closely the results from with the present analysis. 5.4. Depth factors for strip and square foundations The comparison for the variation of the depth factors were made for two extreme values of L/B, that is, square and strip foundations with ϕ = 30° and 40°. This comparison has been presented in Figs. 8 and 9 associated with the factors, dc and dγ, respectively. The values of the depth factor (dc) obtained from the present analysis were compared with the correlations given by Meyerhof [1],

11

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 5. The variation of dγ with L/B, D/B and ϕ.

Hansen [2] and Vesic [3]. Note from Fig. 8 that for a strip foundation, the current values of dc lie between the results of Meyerhof [1] and Hansen [2]. On the other hand, for a square foundation, the recommendations of Meyerhof [1], Hansen [2] and Vesic [3] are found to be quite conservative as compared to the present solution. Amongst Meyerhof [1], Hansen [2] and Vesic [3], the recommendation of Hansen [2] is found to be the most conservative. The values of the depth factor dγ from the present analysis, for strip and square foundations, were compared with the solutions given by (i) Meyerhof [1] and Hansen [2], and (ii) Lyamin et al. [17] using the three-dimensional finite element limit analysis. In the references [1] and [2], the effect of the unit weight of soil mass on the bearing

capacity was calculated as a superposition of the soil overburden pressure along the footing base and self-weight of soil mass below the footing While making the comparison, the depth factor given in the references [1] and [2] was modified accordingly as per the definition adopted in the present study. The comparison of dγ for ϕ = 30° and 40° is presented in Fig. 9. It can be noted for both strip and square foundations, the results from the present analysis compare very closely with the average bound solution of Lyamin et al. [17]. For a square foundation, the present solution lies in between the lower and upper bound curves of Lyamin et al. [17]; for a strip foundation, Lyamin et al. [17], however, did not report the lower and upper bound values of dγ. For a strip foundation, the difference between the different solutions has been

12

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 6. The variation of dq with L/B, D/B and ϕ for γD/q = 10 and 20.

Table 7 A comparison of the bearing capacity factors Nc' , Nq' and N ' for strip and circular footings based on the slip line method. ϕ (°)

0 10 20 30 40 45 a b c

Nq'

Nc'

N'

Circular Footinga

Circular Footingc

Strip Footinga

Strip Footingb

Circular Footinga

Circular Footingc

Strip Footinga

Strip Footingb

Circular Footinga

Circular Footingc

Strip Footinga

Strip Footingc

6.167 11.265 24.015 63.581 231.678 528.386

6.05 11.09 23.67 62.70 228.62 520.30

5.170 8.419 14.964 30.489 76.457 136.198

5.141 8.036 14.839 30.139 75.313 133.874

– 2.984 9.730 37.662 195.188 528.894

1.0 2.94 9.45 36.50 189.19 502.74

– 2.484 6.446 18.603 65.154 137.198

1.0 2.417 6.401 18.401 64.195 134.874

0.00 0.363 2.601 16.430 127.299 426.344

0.00 0.32 2.41 15.54 124.10 419.44

0.00 0.461 2.947 15.361 87.729 239.989

0.00 0.433 2.839 14.754 85.566 234.213

Bearing capacity factors of circular and strip footings obtained by using the present analysis. Bearing capacity factors of strip footings obtained by the analytical method [4]. Bearing capacity factors of circular and strip footings obtained by using the method of stress characteristics [9].

13

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Table 8 A comparison of the factor Nc' for a square footing with ϕ = 0 and D/B = 0. Study

Methods

Bearing capacity factors, Nc'

Michalowski [12]

Upper bound limit analysis

Salgado et al. [15]

FELA

Gourvencec et al. [33] Nguyen and Merifield [34] Mohapatra and Kumar [11] Present Analysis

FEA-Abaqus FEA-Abaqus RPIM based upper bound limit analysis CS-FEM based upper bound limit analysis

Continuous mechanism, UB: 6.830 Multi-block mechanism: 4 symmetry planes, UB: 6.823 2 symmetry planes, UB: 6.561 LB: 5.523 UB: 6.221 5.91 5.950 6.124 5.908

Fig. 7. A comparison of the shape factors, sc, sq and sγ for ϕ = 30° and 40°.

14

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 8. A comparison of depth factors dc for ϕ = 30° and 40°.

found to be only marginal. On the other hand, for a square foundation, the differences between the different solutions become relatively much more extensive. For a square foundation, the results of Meyerhof [1] for dγ become the most conservative. The expressions for the shape and depth factors obtained by Meyerhof [1] have been determined empirically based on the experimental data of Meyerhof [35,36] and some other experimental results. The shape and depth factor formulas of Hansen [2] and Vesic [3] are based on experimental observations in combination with the application of the slip line method. These expressions generally tend to provide conservative estimate as compared to the different numerical results [11,13–14,16–18]. Unlike the conventional theories, in the present analysis the soil mass above the footing level is not simply modelled as an equivalent surcharge and the shear strength of the soil overburden above the footing level has been duly taken into account while determining the bearing capacity factor N ' . However, the variation of the soil internal friction angle with the stress level has not been taken into account; it understood that at very low stress level, the value of the ϕ tends to become greater [37]. Further, the difference between the values of ϕ under plane strain condition and axi-symmetric conditions has also not been incorporated.

of the foundation, for different combinations of L/B, D/B and ϕ, the contour plots were drawn for the variation of (i) the internal plastic power dissipation per unit volume (dp) normalized with respect to c, (ii) the maximum plastic shear strain rate max normalized with respect to V0/B, and (iii) the nodal velocities normalized with respect to V0. The variation of the internal power dissipation function has been depicted in Fig. 10(a)–(f) while computing Nc' for three different values of L/B, namely, 1, 5 and infinity, and two different values of ϕ, namely, 0° and 40° and with D/B = 5. It is observed that the maximum power dissipation occurs usually in a region below and around the footing edge. For ϕ = 0° and D/B = 5, irrespective of the value of L/B, the failure zone, that is, the region involving the zone of the significant plastic power dissipation, develops only locally around the footing. This observation justifies the observed trend of a constant value of Nc' beyond D/B = 1–2. Note that for ϕ = 40°, even for D/B = 5, the failure zones extend almost up to the ground surface justifying a continuous increase of dc with depth. Since the footing base is assumed to be rough, in a very small region immediately below the footing, no power dissipation occurs indicating the existence of a non-plastic zone below the rough base. Due to development of the plastic strain, no power dissipation can occur in a cohesionless (c = 0) soil medium. In order to assess the failure patterns in such cases, one can examine the distribution of the absolute value of the maximum shear strain rate ( max ) ; from the known strain rate vector at a point, one can determine the principal values of strain rates ( 1, 2 and 3 ), and accordingly, the absolute maximum value

6. Power dissipation function, maximum shear strain rate distribution and nodal velocity pattern plots To understand the collapse mechanism during ultimate shear failure

15

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 9. A comparison of the depth factor dγ for ϕ = 30° and 40°.

' of the shear strain rate is given by max = | 1 3|. While determiningN , the variation of max has been plotted in Fig. 10(g) and (h) for two values of L/B, 1 and 5, and ϕ = 40° with D/B = 5. Note the development of a curvilinear plastic shear zone below and around the footing. Once again, a small non-plastic region develops immediately below the footing base wherein the value of max becomes equal to zero. For greater values of D/B, the failure region spreads continuously outward and eventually reaching ground surface justifying a continuous increase of dγ with D/B. To examine the movement patterns, nodal velocities were shown in Fig. 11 while determining Nc' and N ' for ϕ = 30° using (i) two different values of L/B, namely, 1 and 5, and (ii) two different values of D/B, namely, 1 and 3. It should be mentioned that for drawing the nodal velocity patterns, a much coarser mesh was employed for better visualization. Note that for a small region just below as well as above the footing base, the movement takes place primarily in the vertical downward direction with the velocity (V0) same as that of the footing itself. Along the boundary of the curvilinear shear zone, the movement occurs in a direction away from the footing base and towards ground

surface. Along ground surface, heave can be clearly noted. 7. Remarks This study provides the estimation of the ultimate bearing resistance of the foundations placed at shallow to medium depths. The effect of soil overburden pressure and the shear strength due to overlying soil mass at the footing level has already been accounted in the analysis while determining the bearing capacity factor N ' . The footing has been modelled as a zero-thickness rigid plate embedded completely within soil mass. In reality, the footing and the associated vertical shaft (column) have some finite thickness. The error induced by this choice of foundation modelling is not very extensive. An artifact of this modeling approach is that, for the case of = 0° , a local failure mechanism around the footing develops, which does not extend up to the free surface as shown in Fig. 10(a), (b) and (e). This type of failure mechanism will, however, not form for the actual footings in the presence of shaft.

16

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 10. The variation of the plastic power dissipation (dp/c) and the absolute maximum value of the shear strain rate distibution close to footing base for difffernet material parameters and different values of L/B with D/B = 5.

17

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

Fig. 11. Nodal velocity patterns for different combination of L/B amd D/B and using different soil input parameters.

18

Computers and Geotechnics 117 (2020) 103275

D. Mohapatra and J. Kumar

8. Conclusions

[12] Michalowski RL. Upper bound load estimates on square and rectangular footings. Géotechnique 2001;51(9):787–98. [13] Zhu M, Michalowski RL. Shape factors for limit loads on square and rectangular footings. J Geotech Geoenviron Eng 2005;131(2):223–31. [14] Puzakov V, Drescher A, Michalowski RL. Shape factor sγ for shallow footings. Geomech Eng 2009;1(2):113–20. [15] Salgado R, Lyamin AV, Sloan SW, Yu HS. Two and three dimensional bearing capacity of foundations in clay. Géotechnique 2004;54(5):297–306. [16] Salgado R. The engineering of foundations. Boston, USA: McGraw-Hill; 2006. [17] Lyamin AV, Salgado R, Sloan SW, Prezzi M. Two and three-dimensional bearing capacity of footings in sand. Géotechnique 2007;57(8):647–62. [18] Antao AN, Silva MV, Guerra N, Delgado R. An upper bound-based solution for shape factors of bearing capacity of footings under drained conditions using a parallelized mixed finite element formulation with quadratic velocity fields. Comput Geotech 2012;41:23–35. [19] Gourvenec SM, Mana DSK. Undrained vertical bearing capacity factors for shallow foundations. Géotech Lett 2011;1(4):101–8. [20] Tapper L. Bearing Capacity of Perforated Offshore Foundations Under Combined Loading DPhil Thesis UK: Worcester College, University of Oxford; 2013 [21] Benmebarek S, Saifi I, Benmebarek N. Depth factors for undrained bearing capacity of circular footing by numerical approach. J Rock Mech Geotech Eng 2017;9(4):761–6. [22] Chen JS, Wux CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. Int J Numer Meth Eng 2001;50:435–66. [23] Le CV, Nguyen-Xuan H, Askes H, Bordas S, Rabczuk T, Nguyen-Vinh H. A cell-based smoothed finite element method for kinematic limit analysis. Int J Numer Methods Eng 2010;88(12):1651–74. [24] Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. Comput Mech 2007;39:859–77. [25] Le CV. Estimation of bearing capacity factors of cohesive-frictional soil using the cell-based smoothed finite element method. Comput Geotech 2017;83:178–83. [26] MOSEK, The MOSEK Optimization Software, 2016. (Available from: http://www. mosek.com). [27] Martin CM, Makrodimopoulos A. Finite element limit analysis of Mohr Coulomb material in 3D using semidefinite programming. J Eng Mech 2008;134(4):339–47. [28] Chen WF, Han DJ. Plasticity for Structural Engineering. New York: Springer-Verlag; 1988. [29] Mohapatra D, Kumar J. Upper bound finite elements limit analysis of axisymmetric problems for Mohr-Coulomb materials using semidefinite programming. J Eng Mech 2018;144(7). https://doi.org/10.1061/(ASCE)EM.1943-7889.0001472. [30] Xuan HN, Nguyen HV, Bordas S, Rabczuk T, Duflot M. A cell-based smoothed finite element method for three-dimensional solid structures. KSCE J Civ Eng 2012;16(7):1230–42. [31] Makrodimopoulos A, Martin CM. Upper bound limit analysis using simplex strain elements and second-order cone programming. Int J Numer Anal Methods Geomech 2007;31(6):835–65. [32] Martin CM. The use of adaptive finite-element limit analysis to reveal slip-line fields. Géotech Lett 2011;1(2):23–9. [33] Gourvenec SM, Randolph M, Kingsnorth O. Undrained bearing capacity of square and rectangular footings. Int J Geomech 2006;6(3):147–57. [34] Nguyen VQ, Merifield RS. Two- and three-dimensional undrained bearing capacity of embedded footings. Australian Geomechanics 2012;47:25–40. [35] Meyerhof GG. The ultimate bearing capacity of foundations. Géotechnique 1951;2(4):301–32. [36] Meyerhof GG. Some recent research on bearing capacity of foundations. Can Geotech J 1963;1(1):16–26. [37] Parry RHG. Mohr Circles, Stress Paths and Geotechnics. London: Toyler & Fancis group; 2005.

By making use of the cell based smoothed finite element method, the bearing capacity factors for rectangular, strip and circular rigid footings, embedded at shallow depths with D/B ≤ 5, have been determined on the basis of the quasi-kinematic limit analysis. The MohrCoulomb yield criterion, in its original form without any smoothing either in pi or meridional plane, has been employed. The optimization has performed by using the semidefinite programming (SDP) technique. Based on the three-dimensional computational analysis, the results in the form of non-dimensional shape and depth factors, for several combinations of L/B, D/B and ϕ, were obtained. Once the shape factors are simply defined as a function of L/B and ϕ, the depth factors for a given ϕ become a function of both L/B and D/B. The shape factors sc and sq decrease continuously with an increase in L/B. On the other hand, the factor s follows slightly a different trend. The depth factors increase continuously with depth except for ϕ = 0 for which case the depth factor dc becomes constant after a certain value of D/B (around 1–2); for ϕ = 0 at greater depths, the failure zone develops locally around the footing base, on the other hand, in all the other cases, the failure surface almost extends up to ground surface. As compared to the factors dc and dq , the factor d increases much more rapidly with D/B. A non-plastic zone was invariably observed immediately below the base in all the cases for the chosen rough footing. References [1] Meyerhof GG. Some recent research on the bearing capacity of foundation. Can Geotech J 1963;1(1):16–26. [2] Hansen JB. A revised and extended formula for bearing capacity. Geoteknisk Inst Bulletin 1970;28:5–11. [3] Vesic AS. Analysis of ultimate loads of shallow foundations. J Soil Mech Found Div 1973;99(1):45–73. [4] Prandtl L. Über die härte plastischer körper, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen. Mathematisch-Physikalische Klasse 1920;12:74–85. [5] Reissner H. Zum Erddruckproblem. In: Proc 1st Int Cong Applied Mechanics, Delft; 1924. pp. 295–11. [6] Kumar J, Rao VBKM. Seismic bearing capacity factors for spread foundations. Géotechnique 2002;52(2):79–88. [7] Kumar J, Rao VBKM. Seismic bearing capacity of foundations on slopes Géotechnique 2003;53(3):347–61. [8] Kumar J. N for rough strip footing using the method of characteristics. Can Geotech J 2003;40(3):669–74. [9] Martin CM. ABC: Analysis of bearing capacity, 2004. http://www.eng.ox.ac.uk/ civil/people/cmm/software. [10] Hjiaj M, Lyamin AV, Sloan SW. Numerical limit analysis solutions for bearing capacity factor Nγ. Int J Solids Struct 2005;42:1681–704. [11] Mohapatra D, Kumar J. Collapse loads for rectangular foundations by three dimensional upper bound limit analysis using radial point interpolation method. Int J Numer Anal Methods Geomech 2019;43(2):641–60.

19