Estimation of bearing capacity factors of cohesive-frictional soil using the cell-based smoothed finite element method

Computers and Geotechnics 83 (2017) 178–183

Contents lists available at ScienceDirect

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

Technical Communication

Estimation of bearing capacity factors of cohesive-frictional soil using the cell-based smoothed finite element method Canh V. Le Department of Civil Engineering, International University - VNU HCMC, Viet Nam

a r t i c l e

i n f o

Article history: Received 22 July 2016 Received in revised form 5 October 2016 Accepted 27 October 2016

Keywords: Limit analysis Strain smoothing Bearing capacity Cohesive-frictional SOCP

a b s t r a c t Numerical limit analysis provides a powerful tool for assessing the bearing capacity of geo-structures. Recently, a quasi-kinematic procedure, that uses the cell-based smoothed finite elements and second order cone programming, has been developed for limit analysis of purely cohesive soil. In this paper a more generally applicable formulation associated with Morh-Coulomb failure criterion, capable of determining bearing capacity of cohesive-frictional soil, is presented. Computed bearing capacity factors of geo-technical engineering problems such as footings and slope stability are presented and compared with selected published results. It is shown that the current procedure can provide solutions of accuracy being at least as good as those obtained using existing finite element limit analysis, while less computational effort is used. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction In general, analyzing the stability of geo-structures can be performed using one of the four main methods: slip-line methods [1], limit equilibrium [2,3], limit analysis [4,5], and the incremental displacement-based finite element method [6,7]. However, numerical limit analysis procedures based on finite elements and mathematical programming have been found to be the most effective for predicting the collapse load of a geo-structure at limit state [8]. Indeed, the method combines the advantages of numerical discretization with the power of bound theorems, and hence enables the solutions of stability problems with arbitrary geometries, heterogeneous soil profiles, complicated loadings and anisotropic material properties. For elastic-perfectly plastic materials that obey the associative normality rule, the actual collapse load can be bracketed from above and below using upper and lower bound theorems, respectively [9]. In the finite element lower bound formulation, the assumed stress fields are often expressed in terms of nodal stress values and required to satisfy a priori equilibrium conditions within elements and at their interfaces [10,11]. Due to these equilibrium conditions, the construction of such fields is often difficult. Compared with the equilibrium approach, the displacement formulation is more popular due to the facts that compatibility conditions can be satisfied straightaway in the assembly scheme, and

E-mail address: [email protected] http://dx.doi.org/10.1016/j.compgeo.2016.10.023 0266-352X/Ó 2016 Elsevier Ltd. All rights reserved.

that essential (kinematic) boundary conditions can be enforced directly [12–18]. Note that when low-order elements are used, the flow rule can be guaranteed to be satisfied everywhere in the problem domain, and hence strict upper bound on the true bearing capacity can be obtained. However, in order to avoid volumetric locking appeared in the upper bound limit analysis of purely cohesive soil discontinuous velocities must be introduced [13,18]. In recent years, the development and application of strain smoothing finite element methods have attracted much attention due to their properties, e.g. providing accurate solutions, locking free and computational inexpensive [19]. In these methods, smoothed strains are determined by spatially averaging of the standard/compatible strain fields using the divergence theorem. Finite element mesh is divided into smoothing cells, and then by applying the divergence theorem, interior integration on the smoothing cells is transformed into boundary integration. The strain smoothing finite element methods have been developed for limit analysis of structures associated with kinematic theorem [20,21]. Note that the cell-based and edge-based smoothed finite elements have been applied to quasi-kinematic limit analysis of footing problem associated with von Mises yield criterion [20,21]. It has been shown that both methods can prevent volumetric locking phenomena, but the cell-based method with one smoothing cell can provide better approximation of bearing capacity factors than the edge-based method. It is, therefore, relevant to investigate the performance of the cell-based strain smoothing finite elements for limit analysis of cohesive-frictional soil materials.

179

C.V. Le / Computers and Geotechnics 83 (2017) 178–183

The objective of this paper is to extend the recently developed cell-based smoothed finite element (CS-FEM) limit analysis procedure for von Mises yield criterion of [20] to the Morh-Coulomb failure criterion, allowing the calculation of bearing capacity factors for cohesive-frictional soil. Taking advantages of the primal-dual interior-point optimization algorithm, the CS-FEM discrete formulation is also cast in the form of a second-order conic programming to enable solutions to be obtained rapidly. The CS-FEM based numerical procedure is applied to assess the bearing capacity factors of geo-structures such as footings and stability analysis of slopes. The results of the analyses show that the accuracy of the present CS-FEM solution appears to be at least as good as those realizable using finite element limit analysis [14,15], while the number of elements used in the present numerical procedure is very much smaller that used in [14,15]. 2. CS-FEM based kinematic formulation The kinematic theorem is based on the principle of a kinematically admissible velocity field that satisfies the velocity boundary conditions and the plastic flow rule. Consider a rigid-perfectly plastic body of area X 2 R2 and subjected to surface tractions q. Let u be plastic velocity or flow fields that belong to a space U of kinematically admissible velocity fields. The upper bound on the actual collapse multiplier kþ can be determined by solving the following optimization problems [20]

kþ ¼ min DðuÞ;

ð1Þ

u2C

where C ¼ fu 2 U j FðuÞ ¼ 1g; DðuÞ ¼ maxr2B aðr; uÞ is the plastic dissipation rate, B is the yield condition, and FðuÞ; aðr; uÞ are the external and internal virtual work respectively given by

Z

FðuÞ ¼

qT u dX

X

aðr; uÞ ¼

Z X

ð2Þ

rT  dX ¼

Z X

rT ru dX

ð3Þ

where r is the stress vector,  is the strain rate vector given by

 ¼ ½xx ; yy ; cxy T ¼ ru, in which r is the differential operator.

Making use of Morh-Coulomb failure criterion and associated flow rule, the power of plastic dissipation can be determined by [14]

DðÞ ¼

Z

X

tc cos u dX

ð4Þ

where kqk 6 t, with q is the auxiliary variable vector defined as

q ¼ ½q1 ; q2  ¼ ½ðxx  yy Þ; cxy  , and k  k denotes the Euclidean T

T

norm, i.e. kv k ¼ ðv T v Þ . Hence, the upper-bound limit analysis problem for plane strain can be expressed explicitly as 1=2

Z kþ ¼ min tc cos u dX  F 0 ðuÞ X 8 kqk 6 t > > > < xx þ yy ¼ t sin u s:t: > u¼0 on Cu > > : FðuÞ ¼ 1

~ ~ h ¼ Bd

ð6Þ

where d ¼ ½u1 ; v 1 ; . . . ; un ; v n T is nodal velocity vector in each ele~ is the strain-displacement matrix. ment and B Using the smoothed strain rates, the optimization problem (5) can be rewritten as

kþ ¼ min

Nc X

Aj t j c cos u  F 0 ðdÞ

j¼1

8 d¼0 on Cu > > > > < FðdÞ ¼ 1 s:t: ~ j j ~ yy > Bxx d þ B d ¼ t j sin u j ¼ 1; 2; . . . ; N c > > > : kq k 6 t j ¼ 1; 2; . . . ; N c j j

ð7Þ

where N c is the number of smoothing cells over the problem domain; the total number of variables of problem (7) is equal to 2  n þ 3  N c , and

"

2

3

# j j ~ xx ~ yy q1j B dB d 5; qj ¼ 2 ¼ 4 j qj ~ xy d 2B

ð8Þ

The problem (7) will be set up using MATLAB, and solved using the Mosek optimization solver. 3. Results and discussion In this section, the performance of the proposed adaptive scheme is assessed via a number of benchmark problems for which upper- and/or lower-bound solutions have previously been reported in the literature. It should be noted that the cell-based smoothed finite element method with only one smoothing cell per element (CS-FEM1) offers a good combination of accuracy and computational efficiency for plastic limit analysis problems [20], and hence it will be used to solve all problems examined in this section. 3.1. Strip footing bearing capacity The first example comprises the classical plane strain problem originally investigated by Prandtl [22]. On a weightless cohesivefrictional soil (c P 0; u P 0; c ¼ 0) in the absence of surcharge, the Prandtl collapse pressure for a surface footing is given as qf ¼ c N c , where N c is a dimensionless bearing capacity multiplier depending on u. The analytical value of N c can be expressed as [22]

h p u i Nc ¼ ep tan u tan2 þ  1 cot u 4 2

ð9Þ

Taking advantage of symmetry, only half of the foundation is considered. The rectangular region of L and H should be chosen sufficiently large to ensure that rigid elements show up along the

ð5Þ

where c is the cohesive parameter, u is the internal friction angle of soil, F 0 ðuÞ is the external work of loads not subjected to kþ . Now, problem (5) will be discretized using the cell-based finite element method. In the settings of the cell-based strain smoothing finite element method, the smoothed version of the strain rates can be expressed as [20]

Fig. 1. Finite element mesh of 2560 Q4 elements for CS-FEM model and displacement boundary conditions.

180

C.V. Le / Computers and Geotechnics 83 (2017) 178–183

limit analysis procedures using three node triangular element (FEM-T3) and the edge-based smoothed finite element (ES-FEM) limit analysis method presented in [21] are also performed. Collapse multipliers and associated errors for various discretization meshes and methods are shown in Table 1. The convergence rate is also illustrated by Fig. 3. It can be observed that the CS-FEM method provides more accurate solutions and converges faster than FEM-T3 and ES-FEM methods. The evaluation of the bearing capacity factor N c for cohesivefrictional soil is next considered. Note that as the internal friction angle increasing, the width of the rectangular modeling regions, L, should be chosen so that far away from the footing the soil remains rigid and the velocity is equal to zero. Computed collapse factor for various internal friction angles ranging from 5° to 35° are shown in Table 2 and Fig. 4. It is evident that for the case when u ¼ 35° solution obtained using the ES-FEM model with 5760 T3 elements is lower than the result of 50.85 reported in [14] using up to 18,719 T3 elements with discontinuities, and the CS-FEM approach using 5500 quadrilateral elements can provide more accurate solution than the 6-node triangular element method presented in [14] employing a mesh of 18,719 elements, 46.143 com-

Fig. 2. Finite element mesh of 2560 T3 elements for ES-FEM model and displacement boundary conditions.

entire boundary. The punch is represented by a uniform vertical load and appropriate boundary conditions were applied as shown in Figs. 1 and 2. The performance of the proposed method for purely cohesive soil is studied first. The exact collapse multiplier for purely cohesive soil (u ¼ 0) is N c ¼ 2 þ p. For comparison purpose, kinematic

Table 1 The strip footing problem: bearing capacity factor for purely cohesive soil. N

FEM-T3

360 1000 2560

ES-FEM Error (%)

Nc

Error (%)

Nc

Error (%)

5.305 5.235 5.198

3.18 1.81 1.10

5.255 5.207 5.181

1.61 1.27 0.77

5.150 1.146 5.144

0.17 0.09 0.05

5.32

5.28 5.26 5.24 5.22 5.2

FEM−T3 ES−FEM CS−FEM

0.5

0

−0.5

−1

10

5.18

1 log (Relative error in collapse multiplier)

FEM−T3 ES−FEM CS−FEM

5.3

5.16 5.14 0

CS-FEM

Nc

500

1000 1500 2000 number of nodes

2500

−1.5 2.5

3000

3 log (number of nodes)

3.5

10

Fig. 3. The strip footing problem: convergence analysis.

Table 2 The strip footing problem: bearing capacity factors for cohesive-frictional soil. Internal friction angles (u)

Models 5

10

15

20

25

30

35

ES-FEM Error (%)

6.659 2.62

8.588 2.91

11.327 3.19

15.348 3.47

21.490 3.71

31.380 4.12

47.978 (N ¼ 5760) 4.03

CS-FEM Error (%)

6.491 0.04

8.345 0.001

10.977 0.005

14.834 0.000

20.724 0.017

30.143 0.011

46.143 (N ¼ 5500) 0.05

Exact

6.4888

8.3449

10.9765

14.8340

20.7205

30.1396

46.12

181

C.V. Le / Computers and Geotechnics 83 (2017) 178–183

50 45 40

CS−FEM ES−FEM Exact solution

35 30 25 20 15 10 5 5

10

15

20

25

30

35

Fig. 4. The strip footing problem: bearing capacity factors of cohesive-frictional soil.

pared with 46.37. Furthermore, although it has been pointed out that the strain smoothing finite element based procedure cannot be guaranteed to provide strict upper bound solutions, it is clear that all the solutions obtained are higher than the exact values. The efficacy of the proposed CS-FEM method is also studied. Table 3 compares computational effort between the CS-FEM method and the finite element limit analysis using 6-node triangular elements [14]. It is evident that the present CS-FEM based numerical limit analysis can produce accurate solutions with minimal computational effort. Now, the performance of the proposed strain smoothing finite elements is tested by a more challenging problem, the calculation of N c for weighty cohesionless soil without surcharge (c ¼ 0; u > 0; c > 0). The self-weight bearing capacity problem was investigated analytically using the limit equilibrium method [23–25] and the method of stress characteristics [26,27], and finite element limit analysis [15,14,29]. Computed values of self-weight bearing capacity N c using the strain smoothing finite elements for both smooth and rough soil-footing interfaces and various fric-

Table 3 On the performance of the proposed method: u ¼ 35°. Present CS-FEM solutions

a

Results in [14]

Elements

CPU time (sa)

Nc

Elements

CPU time (s)

Nc

768(Q4) 1980 5500

0.3 0.7 3.2

47.44 46.19 46.14

774(T6) 6308 18719

1.4 22.5 36.5

49.25 46.99 46.37

Solution time taken by MOSEK solver on a 2.8 GHz Pentium 4 PC.

Table 4 The strip footing problem: bearing capacity factors for weighty cohesionless soil.

u°

Smooth interface

5 10 15 20 25 30 35 40 45

N

Rough interface

ES-FEM

CS-FEM

ES-FEM

CS-FEM

0.0974(15.3) 0.3134(11.6) 0.7803(11.6) 1.7525(11.0) 3.8120(10.1) 8.4469(10.4) 19.509(11.0) 48.261(11.7) 133.52(13.5)

0.0934(9.6) 0.3021(6.1) 0.7388(5.5) 1.6645(4.8) 3.6198(3.2) 8.0031(3.5) 18.131(3.1) 44.239(2.3) 118.73(1.1)

0.1666(46.9) 0.5654(30.5) 1.4764(25.0) 3.4069(20.0) 7.6336(17.6) 17.215(16.7) 40.032(16.1) 102.30(19.6) 303.60(29.6)

0.1420(25.2) 0.5001(15.4) 1.3135(11.2) 3.0512(7.7) 6.9321(6.8) 15.411(4.5) 35.899(4.1) 88.375(3.3) 241.24(3.0)

10240 10240 16464 16464 18816 19360 21060 21060 21060

ðÞ indicates the computed errors based on exact solutions reported in [28].

140

350

120

300

CS−FEM ES−FEM Hjiaj (2005) Makrodimopoulos et al. T6 Makrodimopoulos et al. T3 Exact solution

100 80

250 200

60

150

40

100

20

50

0 0

10

20

30

CS−FEM ES−FEM Hjiaj (2005) Makrodimopoulos et al. T6 Makrodimopoulos et al. T3 Exact solution

40

50

0 0

10

20

30

Fig. 5. The strip footing problem: bearing capacity factors for cohesionless-frictional soil (N c ).

40

50

182

C.V. Le / Computers and Geotechnics 83 (2017) 178–183

tion angles ranging from 5° to 45° are presented in Table 4. In Fig. 5, the present self-weight bearing capacity N c are also compared with upper bound solutions reported in [15] using threenode triangular elements with discontinuities and with those obtained in [14], where 6-node triangular elements were used. In general, the ES-FEM results are in good agreement with those of [14] employing discontinuous three-node triangular elements, while our CS-FEM solutions agree with those reported in [15,14] (with the computation using 6-node triangular elements). It is observed in the extreme case when u ¼ 45° that the present CSFEM solutions are lower than upper bounds of [14] (for both smooth and rough footings) and of [15] (for smooth footing) while less computational cost is used, e.g. 21,060 quadrilateral elements compared with 31,481 6-node elements in [14] and with 20,815 3node elements and 31,119 discontinuities in [15]. 3.2. Slope stability The second example is the classical stability of a homogeneous slope in a cohesive-frictional soil having different inclinations from 50° to 90° and height of H. In this typical slope stability analysis the only load considered is the soil weight, c. The quantity of interest is

80 70

β= 50 β= 60 β= 70 β= 80 β= 90

60 50 40 30 20 10 0 0

10

20

30

40

Fig. 7. Slope analysis: factors of safety for various slope angles and internal friction angles.

14

cH

the factor of safety or the dimensionless critical height, N s ¼ c . The factors of safety for this problem have been investigated analytically by Chen [4] and numerically by Lyamin and Sloan [17], Krabbenhoft et al. [18] and Makrodimopoulos and Martin [14]. Note that solutions presented in [18] are exactly identical to those reported in [17], but the applied formulations are different. In this example, AUTOMESH-2D [30] is used to generate quadrilateral meshes for slopes with different inclined angles, see Fig. 6 for the finite element mesh of 70° slope. The results of analyses

Present CS−FEM solutions Results by Krabbenhoft et al. (2005) Results by Chen et al. (1990)

13 12 11 10 9 8 7 6 5 50

2

1.5

70

80

90

Fig. 8. Slope analysis: the present CS-FEM factors of safety in comparison with other results.

1

0.5

0

0

60

0.5

1

1.5

2

2.5

3

Fig. 6. Slope of 70° inclined angle: finite element mesh of 3708 quadrilateral elements and 3793 nodes.

for wide ranges of slope angles and internal friction angles are shown in Table 5 and Fig. 7. In the case when internal friction angle u ¼ 20°, the present CSFEM solutions are compared with results reported by Krabbenhoft et al. [18] and by Chen [4], as shown in Fig. 8. It is observed that for all slope angles the present CS-FEM factors of safety improve on the solutions of [18], and slightly higher than analytical results of [4]. For b ¼ 70° and u ¼ 20°, 35° our solutions are also lower than upper bounds of [14] using unstructured meshes of 3-node and 6-node elements, despite the fact that the number of elements used in our computation is very much smaller than that of [14],

Table 5 Slope stability: factors of safety for various slope angles and internal friction angles using CS-FEM. Internal friction angle u°

b

50 60 70 80 90

0

5

10

15

20

25

30

35

40

5.353 5.163 4.763 4.305 3.821

6.874 6.141 5.457 4.819 4.193

8.529 7.268 6.257 5.374 4.594

10.704 8.666 7.200 6.042 5.039

13.752 10.451 8.337 6.796 5.534

18.302 12.819 9.753 7.678 6.104

25.746 16.190 11.570 8.747 6.774

39.954 21.232 14.014 10.103 7.584

73.523 29.722 17.539 11.854 8.570

C.V. Le / Computers and Geotechnics 83 (2017) 178–183

183

Fig. 9. Analysis of 70° slope: plastic dissipation distribution for various friction angles.

only 3708 quadrilateral elements compared with 28,864 triangular elements. Plastic dissipation distribution for 70° slope with various friction angles are plotted in Fig. 9. It can be seen that the size of plastic zones decreases as the friction angle increases. 4. Conclusions The cell-based smoothed finite element limit analysis described in [20] has been extended to allow the computation of bearing capacity factors of cohesive-frictional soil governed by the MorhCoulomb failure criterion and obeyed the associated flow rule. Solutions for the ultimate bearing capacity of footings and slope stability with a wide range of internal friction angles are obtained. It has been demonstrated that the CS-FEM method can result in very accurate estimation of bearing capacity factors with a minimal computational effort, even for cases of high-friction angles, e.g. u P 35°. Obtaining solutions of comparable accuracy using finite element limit analysis procedures would necessitate the use of a very much finer mesh, and hence longer run times. Acknowledgement The authors acknowledges the support of Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant reference 107.01-2015.41. References [1] Sokolovskii VV. Statics of granular media. Oxford, UK: Pergamon Press; 1965. [2] Bishop AW. The use of slip circles in the stability analysis of earth slopes. Geotechnique 1955;5:7–17. [3] Narita K, Yamaguchi H. Bearing capacity analysis of foundations on slopes by use of log-spiral sliding surfaces. Soils Found 1990;30:144–52. [4] Chen WF. Limit analysis and soil plasticity. Amsterdam: Elsevier; 1975. [5] Yu HS. Plasticity and geotechnics. Newyork: Springer; 2006. [6] Sloan SW, Randolph MF. Numerical prediction of collapse loads using finite element methods. Int J Numer Anal Methods Geomech 1982;6:47–76. [7] de Borst R, Vermeer PA. Possibilities and limitations of finite elements for collapse. Geotechnique 1984;34:199–210. [8] Sloan SW. Geotechnical stability analysis. Geotechnique 2013;63:531–72.

[9] Drucker DC, Greenberg W, Prager W. Extended limit design theorems for continuous media. Quart Appl Math 1952;9:381–9. [10] Sloan SW. Lower bound limit analysis using finite elements and linear programming. Int J Numer Anal Methods Geomech 1988;12:61–77. [11] Krabbenhoft K, Damkilde L. Lower bound limit analysis of slabs with nonlinear yield criteria. Comput Struct 2002;80:2043–57. [12] Vicente da Silva M, Antao AN. A non-linear programming method approach for upper bound limit analysis. Int J Numer Methods Eng 2007;72:1192–218. [13] Sloan SW, Kleeman PW. Upper bound limit analysis using discontinuous velocity fields. Comput Methods Appl Mech Eng 1995;127:293–314. [14] 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:835–65. [15] Hjiaj M, Lyamin AV, Sloan SW. Numerical limit analysis solutions for the bearing capacity factor N c . Int J Solids Struct 2005;42:1681–704. [16] Bottero A, Negre R, Pastor J, Turgeman S. Finite element method and limit analysis theory for soil mechanics problems. Comput Methods Appl Mech Eng 1980;22:131–49. [17] Lyamin AV, Sloan SW. Upper bound limit analysis using linear finite elements and non-linear programming. Int J Numer Anal Methods Geomech 2002;26:181–216. [18] Krabbenhoft K, Lyamin AV, Hjiaj M, Sloan SW. A new discontinuous upper bound limit analysis formulation. Int J Numer Methods Eng 2005;63:1069–88. [19] Liu GR, Nguyen-Thoi T. Smoothed finite element methods. NewYork: CRC Press, Taylor and Francis Group; 2010. [20] Le CV, Nguyen-Xuan H, Askes H, Bordas S, Rabczuk T, Nguyen-Vinh H. A cellbased smoothed finite element method for kinematic limit analysis. Int J Numer Methods Eng 2010;83(12):1651–74. [21] Le CV, Nguyen-Xuan H, Askes H, Rabczuk T, Nguyen-Thoi T. Computation of limit load using edge-based smoothed finite element method and secondorder cone programming. Int J Comput Meth 2013;10(01):1340004. [22] Prandtl L. Ueber die haerte plastischer koerper. Nach Akad Wissen Gottingen. II. Math-Phys Klasse II 1920;12:74–85. [23] Meyerhof GG. The ultimate bearing capacity of foundations. Geotechnique 1951;2(4):301–32. [24] Meyerhof GG. Some recent research on the bearing capacity of foundations. Can Geotech J 1963;1(1):16–26. [25] Kumbhojkar AS. Numerical evaluation of Terzaghi’s N c . Int J Geomech ASCE 1993;119(3):598–607. [26] Bolton MD, Lau CK. Vertical bearing capacity factors for circular and strip footings on Mohr-Coulomb soil. Can Geotech J 1993;30(6):1024–33. [27] Kumar J. N c for rough strip footing using the method of characteristics. Can Geotech J 2003;40(3):669–74. [28] Martin CM. Exact bearing capacity calculations using the method of characteristics. Proceedings of the 11th international conference of IACMAG, Turin, vol. 4. p. 441–50. [29] Ukritchon B, Whittle AJ, Klangvijit C. Calculations of bearing capacity factor N c using numerical limit analyses. J Geotech Geoenviron Eng ASCE 2003;129 (5):468–74. [30] AUTOMESH-2D: a robust quadrilateral mesh generator by Dr. Xinwu Ma. Available from .