SciELO - Scientific Electronic Library Online

vol.14 número13Viscoelastic Substrates Effects on the Elimination or Reduction of the Sandwich Structures Oscillations Based on the Kelvin-Voigt ModelExperiment and Calculation Method of the Dynamic Response of Deep Water Bridge in Earthquake índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados




Links relacionados


Latin American Journal of Solids and Structures

versão impressa ISSN 1679-7817versão On-line ISSN 1679-7825

Lat. Am. j. solids struct. vol.14 no.13 Rio de Janeiro  2017 


Static Limit Analysis of Reinforced Soil Structures by a Simple Finite Element and Second-Order Cone Programming

Eric Luis Barroso Cavalcantea  * 

Eliseu Lucena Netob 

Denilson José Ribeiro Sodréc 

a Tribunal de Contas da União, 70042-900 Brasília, DF, Brazil.

b Instituto Tecnológico de Aeronáutica, 12228-900 São José dos Campos, SP, Brazil.

c Universidade Federal do Pará, 66075-110 Belém, PA, Brazil.


To discretize reinforced soil structures in plane strain and predict their collapse load, a simple three-node triangular finite element is formulated based on the static theorem of the limit analysis. The element satisfies the equilibrium equations and the mechanical boundary conditions in a weak sense. A modified Mohr-Coulomb yield surface is adopted to describe the reinforced soil behavior from a macromechanics point of view. It is also taken into account the possibility of tension failure of the reinforcement and failure of the reinforcement interface. The stated nonlinear convex optimization problem is cast as second-order cone programming. Numerical examples illustrate the predictive accuracy of the above scheme as well as the efficiency and speed of an interior-point method to reach optimal solutions.

Keywords: reinforced soil; static theorem; finite element; second-order cone programming


In the design of reinforced soil structures, the maximum load to be resisted at the impending collapse must be evaluated. The ability of the methods to accurately estimate ultimate limit states depends on the fulfillment of theoretical requirements derived from continuum mechanics concerning the equilibrium, strain-displacement relations, constitutive behavior and boundary conditions. Under certain restrictions, the theory of plasticity allows the prediction of the collapse load by means of the limit analysis based on the bound theorems (Drucker et al., 1952). Numerical solutions are most facilitated in this case coupling such theorems with the finite element method, as originally proposed by Lysmer (1970) and Bottero et al. (1980). This is a straightforward approach for computing lower and upper bounds on the load at incipient collapse which have gained increasing attention (Lyamin and Sloan, 2002; Krabbenhoft and Damkilde, 2003; Makrodimopoulos and Martin, 2006; Yu and Tin-Loi, 2006; Ciria et al., 2008; Munõz et al., 2009; Le et al., 2010; Bleyer and de Buhan, 2014; Nguyen-Thoi et al., 2015). This is partly due to the development of robust optimization methods and solvers on which they strongly rely (Drud, 1996; Sturm, 1999; Lyamin and Sloan, 2002; Murtagh and Saunders, 2003; Tütüncü et al., 2003; Andersen et al., 2003; Krabbenhoft and Damkilde, 2003).

The identification of failure modes through experiments has had an essential impact on the development of reinforced soil mechanics. It has been realized that, in the case of sufficiently good bonding between soil and reinforcement, reinforced soil behaves, in the macroscopic scale, like a homogeneous material. This important observation was confirmed by laboratory investigations reported in Long et al. (1972) and Chapuis (1980). Reinforced soil can then be treated as a composite material formed by the association of frictional soil and tension-resistant fiber elements, such as geosynthetics, with behavior presumed homogeneous and anisotropic from a macromechanics point of view. The material strength is estimated from the strength characteristics of its components and the interaction between them. Limit analysis procedures derived from the referred macromechanics approach have been successfully applied to predict the observed response of reinforced superficial foundations and earth walls (Sawicki, 1983; Kulczykowski, 1985; Sawicki and Lesniewska, 1987, 1989; de Buhan et al., 1989; de Buhan and Siad, 1989; Yu and Sloan, 1997).

In this paper, specific attention is focused on the static approach of the limit analysis coupled with the finite element method to provide limit load solutions for reinforced soil structures in plane strain. In the description of the soil behavior, the isotropic Mohr-Coulomb yield surface is modified to include the effect of anisotropy caused by the presence of reinforcement. It is also taken into account the possibility of tension failure of the reinforcement and failure of the soil-reinforcement interface (de Buhan and Siad, 1989). The finite element proposed for limit analysis, which has its roots in the classic paper by Anderheggen and Knöpfel (1972), satisfies the equilibrium equations and the mechanical boundary conditions on average. In this sense, a solution obtained with this weak equilibrium model is not a strict lower bound. The formulation of the static theorem is then written within the framework of a nonlinear convex optimization technique, known as second-order cone programming (Lobo et al., 1998).

There is no analytical approach for the solution of general convex optimization problems, but there are very effective methods for solving them (Boyd and Vandenberghe, 2009). The interior-point method employed in this paper works extremely well in practice. It was developed by Andersen et al. (2003) and implemented in the solver MOSEK (MOSEK ApS, 2016) to address nonlinear convex optimization problems, such as second-order cone programs. Predictions of the collapse load for some well-known stability problems of soil mechanics, namely the bearing capacity of footings and the stability of slopes, are computed after the soil is reinforced and the discrete problem is formulated as second-order cone programming (SOCP). For comparison, the examples are also treated as linear programming (LP) and the results with both schemes compared with exact or approximate solutions available in the literature.


The reinforced soil in this paper will be treated as a homogeneous composite material with anisotropic properties. The reinforcement is assumed to be unidirectional with thickness d very small compared to the space h between two reinforcements (d/h ≪ 1), as shown in Figure 1. Three stress tensors are defined at every point in the homogenized continuum: tensor ⎣σt σt τtn⎦ T for macrostresses, and tensors σs=σtsσnsτtnsT and σr=σtrσnrτtnrT for microstresses which act on the soil and reinforcement, respectively. The tensor components are related by

σt=σts+dhσtr=σts+σr σn=σns=σnrτtn=τtns=τtnr (1)

where σr=d/hσtr (Yu and Sloan, 1997).

Figure 1: Stresses on reinforced soil. 

The relations between the macrostress components in the orthogonal Cartesian coordinate systems xy e tn are

σx=σtcos2θ+σnsin2θ2τtnsinθcosθ σy=σtsin2θ+σncos2θ+2τtnsinθcosθτxy=(σtσn)sinθcosθ+τtn(cos2θsin2θ) (2)

where θ represents the angle between the horizontal axis x and the reinforcement direction t, measured counterclockwise. Substitution of (1) and relations

σts=σxscos2θ+σyssin2θ+2τxyssinθcosθ σns=σxssin2θ+σyscos2θ2τxyssinθcosθτtns=(σxsσys)sinθcosθ+τxys(cos2θsin2θ) (3)

into (2) yields

σxs=σxσrcos2θ σys=σyσrsin2θτxys=τxyσrsinθcosθ. (4)

The reinforcement inside the soil is supposed to act merely as tensile load carrying elements (de Buhan et al., 1989). The bounds

0σrσo (5)

are thereby imposed on σ r , where σo = (d/h)σyield is the tensile yield strength σyield of the reinforcement times its volumetric fraction.

The soil mass with cohesion c and angle of internal friction ϕ is initially assumed to obey the Mohr-Coulomb yield surface in plane strain (Muñoz et al., 2009):

F¯s=(σxsσys)2+(2τxys)2[2ccosϕ(σxs+σys)sinϕ] (6)

This expression is then modified to include the effect of anisotropy caused by the presence of reinforcement using relations (4):

Fs=(σxσyσrcos2θ)2+(2τxyσrsin2θ)2[2ccosϕ(σx+σyσr)sinϕ] (7)

The soil-reinforcement interface failure surface is assumed to be described by

Fi= |τtn| ci+σntanϕi (8)

where ci and ϕi denote the interface cohesion and angle of internal friction, respectively (Yu and Sloan, 1997). In view of (2), the failure surface (8) takes the form

Fi= 12|(σyσx)sin2θ+2τxycos2θ| ci+(σxsin2θ+σycos2θτxysin2θ)tanϕi=0 (9)

in the system xy.


The limit analysis relies on the assumption of an elastic-perfectly plastic material with a flow rule associated to a convex yield surface, and also on small displacement gradients so that the body does not undergo large deformation at collapse (problem readily stated on the undeformed body geometry). The static approach of the limit analysis requires that the assumed stress field must satisfy the equilibrium equations, the mechanical boundary conditions and the yield criterion everywhere. Under these idealized conditions, the computed limit load is a lower bound on the true collapse load (Chen, 1975).

3.1 Equilibrium Equations

Let Ω be the region occupied by the reinforced soil structure subjected to the body force b = ⌊bx by⌋ T expressed in the system xy. The macrostress field σ = ⌊σx σy σxy⌋ T must satisfy the equilibrium equations

Dσ+b=0 (10)

throughout the domain Ω, where the differential operator

D=[x0y0xx] (11)

3.2 Mechanical Boundary Conditions

Let Γ be the boundary of the domain Ω, with unit outward normal vector denoted by n = ⌊nx ny⌋ T . In addition to (10), the macrostress field must also satisfy the mechanical boundary conditions

t=Nσ=t¯ (12)

on the portion Γt of Γ on which the traction (stress vector) t = ⌊tx ty⌋ T is prescribed as being t¯ . Matrix

N=[nx0ny0nynx] (13)

contains the components of n.

3.3 Yield Criterion

To assure that the yield criterion is satisfied it is necessary to impose simultaneously

Fs0 Fi0 0σrσo (14)

As the function Fi defined in (9) has a term in modulus, the second of inequalities (14) splits into two other conditions given by

Fi=Ajσx+Bjσy+Cjτxyci0 j=1,2 (15)


A1=sin2θtanϕi12sin2θ B1=cos2θtanϕi+12sin2θC1=sin2θtanϕi+cos2θ A2=sin2θtanϕi+12sin2θB2=cos2θtanϕi12sin2θ C2=sin2θtanϕicos2θ (16)


Suppose that the reinforced soil is divided into a number of triangular elements and treated as an assembly of them. To apply the above equations to a typical element Ωє shown in Figure 2, the definition of the boundary should be extended to include the traction continuity on the interelement portion Γi:

(t)++(t)=0 (17)

where the superscripts “+” and “-” denote the two sides of Γi. Expressions (10), (12) and (17) may now be viewed as the equilibrium within the element and on its boundary Γt and Γi, respectively (Pian and Wu, 2006; Washizu, 1982).

Figure 2: Finite element Ωє with the unit normal vector n on its boundary Γє. 

Equations (10), (12) and (17) can be enforced to be satisfied on average over the element and on its boundary by means of

ΩewT(Dσ+b) dx dyΓtwT(tt¯) dsΓiwTt ds=0 (18)

where w = ⎣wx wy⎦ T is an arbitrary weight function that is continuous on the element domain and across its interfaces. The last integral when considered jointly with those of the neighborhood elements enforces (17), since it represents one of the terms in ΓiwT[(t)++(t)] ds=0 .

In view of the divergence theorem, one writes

ΩewT(Dσ) dx dy=ΓewTt dsΩe(DTw)Tσ dx dy (19)

where w is additionally assumed to be differentiable once with respect to x and y. Since Γe= Γu ⋃ Γt⋃ Γi, expression (18) is then simplified to

ΓuwTt ds+ΓtwTt¯ ds+ΩewTb dx dyΩe(DTw)Tσ dx dy=0 (20)

where the portion Γu of the element boundary Γe falls on the reinforced soil boundary with prescribed displacement.

The macrostress field is linearly approximated over the element by

σ=N1σ1+N2σ2+N3σ3 (21)


σi=σxiσyiτxyiT Ni=12A(αi+βix+γiy) i=1,2,3 (22)

are the macrostress nodal values and the shape functions, respectively, with A standing for the triangle area. To evaluate

αi=xjykxkyj βi=yjyk γi=xkxj (23)

from the node coordinates, the indices i, j, k should be permuted in a natural order (i ≠ j ≠ k). The weight function is also taken to vary linearly,

w=N1w1+N2w2+N3w3 (24)

with nodal values

wi=wxiwyiT i=1,2,3 (25)

Substitution of (21) and (24) into (20) yields

{w1w2w3}T(F[GGG]{σ1σ2σ3})=0 (26)


G=16[y2y30x3x20x3x2y2y3y3y10x1x30x1x3y3y1y1y20x2x10x2x1y1y2]F=Γt[N1N2N3]Tt¯ ds+Ωe[N1N2N3]Tb dx dy (27)


Ni=[Ni00Ni] (28)

Because (26) holds for any arbitrary weight function, it follows that

[GGG]{σ1σ2σ3}=F (29)

Elimination of the integral over Γu by choosing w null over there introduces zero components of w1, w2 or w3 into (26) which must be accounted for in (29) by removing the respective equations. Thus, the element enforces the equilibrium equations (10) and the mechanical boundary conditions (12) by the discrete equation (29) on average according to (20). Note that the traction continuity (17) is rigorously enforced because the assumed stress field (21) is continuous across the element interfaces. The element development has its roots in the classic paper by Anderheggen and Knöpfel (1972).


In the static limit analysis, the structure equilibrium and the yield criterion, expressed in terms of nodal stresses, are constraints of an optimization problem for the applied load maximization. The optimal solution

λ*={maxλ| lσ¯=λf1+f2 , g1(σ¯, σ¯r)0 , g2(σ¯)0 , 0σ¯rσo} (30)

identifies the collapse load, where the applied load has been split into two parts: ⋋f1 which is adjusted during the optimization by means of the load factor ⋋, and f2 which is kept constant. Vectors σ¯ and σ¯r collect the nodal values of σ and σ r respectively.

Under the continuity provided by the approximations (21) and (24), the equality constraint

lσ¯=λf1+f2 (31)

which arises from the assembly of (29), represents the discrete equilibrium of the whole reinforced soil structure. The inequalities constraints

g1(σ¯, σ¯r)0 g2(σ¯)0 0σ¯rσo (32)

stem from the evaluation at each node of the first, second and third of inequalities (14), respectively, since it is sufficient to enforce them at each node in order that they are satisfied throughout the mesh.

The nonlinear convex optimization problem (30) is not a SOCP problem. To formulate it as such (Lobo et al., 1998), we introduce the auxiliary variables

v={v1v2v3}=[110cos2θ0002sin2θ0sinϕsinϕ0sinϕ2ccosϕ]{σxσyτxyσr1} (33)

into (7) to state Fs ≤ 0 as the three-dimensional second-order cone

v12+v22v3 v30 (34)

sketched in Figure 3. Now, the problem can be treated as SOCP by just replacing g1(σ¯, σ¯r)0 with

h1(σ¯, σ¯r, v¯)=0 h2( v¯)0 v¯30 (35)

where the vectors v¯ and v¯3 collect the nodal values of v and v3. The new constraints (35) stem from the evaluation of (33) and (34) at each node, and all the problem nonlinearity concentrates in the constraint h2 ≤ 0 related to the cone definition.

Figure 3: Second-order cone (34). 

The optimization problem solution is carried out following the steps: (a) set up of the SOCP problem with YALMIP (Löfberg, 2004) in the MATLAB environment: problem (30) with g1 ≤ 0 replaced by (35); (b) solution by MOSEK using the interior-point optimizer developed for nonlinear convex optimization by Andersen et al. (2003).


Predictions of the collapse load for some well-known stability problems of soil mechanics, namely the bearing capacity of footings and the stability of slopes, are computed after the soil is reinforced. This paper focuses on the limit analysis based on the static theorem as well as on the efficiency and speed with which the problems are solved as SOCP. For comparison, the piecewise linearization of Fs adopted by Yu and Sloan (1997) is accomplished, using 48 sides in the linearized yield polygon presented therein, so that the problems can also be solved as LP by performing a sequence of optimization started with the Newton barrier method and concluded with the simplex method at FICO® Xpress Optimization Suite. The role played by the simplex method in such a procedure is to guarantee that the optimal solution is achieved. The computations are performed on a Sony VAIOTM All-in-One machine (i5 dual core 2.50 GHz CPU, 12 GB RAM), running a 64-bit Windows 7. The reported CPU times refer to the time spent on the optimization iterations.

6.1 Strip Footing on Cohesionless Reinforced Soil

For the strip footing of width B indicated in Figure 4a, only half of the full geometry needs to be considered due to symmetry conditions. To deal with the semi-infinite domain, a rectangular portion is discretized as in Figure 4b extending horizontally from the symmetry plane to the right over a length 3B and vertically downwards to a Depth 2B. This domain is sufficiently extended so that conditions at the far boundary do not have significant effect on the calculated limit load.

Figure 4: Strip footing: (a) half of the adopted domain with horizontal reinforcement; (b) unstructured mesh with 2051 elements. 

In this example, f1 present in (30), denotes all the loads applied on Γt and f2 = 0 also present in (30), accounts for the null body forces. The adopted mechanical boundary conditions are

ts=0on the left(symmetric)edgetn=0 ts=0on the upper edge without strip footingtn=qon the upper edge with strip footing (36)

where q is the force per unit area exerted by the footing, tn and ts are the traction components normal and tangential to the edges.

The weight function components wn and ws are imposed to be null at nodes on the edges where the traction components tn and ts are unknown, respectively:

wn=0on the left(symmetric)edgews=0on the upper edge with strip footing (37)

In addition, wn and ws are made to be null on the far right and bottom edges because tn and ts are taken to be unknown over there.

It is well-known that the bearing capacity of strip footing is null when the soil is purely frictional, weightless and subjected to no surcharge (Terzaghi, 1943; Davis and Booker, 1971). The simple inclusion of horizontal reinforcement, however, rises the bearing capacity to

qeσo=(1+sinϕ)e(π2+ϕ)tanϕ (38)

This is the exact solution, as demonstrated by Kulczykowski (1985), when the reinforced soil is treated as a homogeneous and anisotropic composite material with perfectly rough soil-reinforcement interface. Table 1 summarizes our results obtained by simulating the above property using ci = c = 0, ϕ i = ϕ. It is also depicted the number of iterations and the CPU time required to solve the SOCP and the LP problems. The small difference between the SOCP and LP predictions can be reduced increasing the number of sides in the linearized yield polygon. However, the computational cost of an LP solution may become prohibitive for a larger number of sides because of the amount of iterations and CPU time involved. The SOCP number of iterations agrees with already published numerical experiments, in the sense that an interior-point method demands typically between 5 and 50 iterations to solve any SOCP problem (Lobo et al., 1998). In the LP solutions, the simplex phase answers for 99.2% to 99.9% of the total number of iterations. For instance, out of 10103 iterations for ϕ = 10° the count of 10037 takes place in the simplex phase.

Table 1: Bearing capacity of strip footing on cohesionless reinforced soil. 

q/σo iter CPU q/σo iter CPU
10° 1.3370 36 0.95 1.3362 10103 24
15° 1.8972 30 0.84 1.8953 6864 16
20° 2.5239 37 0.95 2.5210 53023 150
25° 3.4577 40 1.01 3.4511 137467 289
30° 4.8029 26 0.73 4.7904 14094 37
35° 7.0282 33 0.83 7.0000 11473 28

iter: number of iterations CPU: CPU time in seconds

The exact solution (38) is graphically compared in Figure 5 to our SOCP predictions. The rigorous lower bounds provided by Yu and Sloan (1997) are also shown, which are based on the discretization of the whole domain into triangular elements devised by Lysmer (1970) along with extension elements developed by Pastor (1978). In that work, as stated previously, the optimization problem is conceived as linear programming. For a given mesh, the scheme adopted by Yu and Sloan leads to a much larger problem size because of the use of Lysmer element and the piecewise linearization of Fs. The gap to the exact solution displayed by our predictions is almost steady for all angles of internal friction considered, whereas it is noticeable the deterioration of Yu and Sloan lower bounds for higher values of ϕ.

Figure 5: Strip footing on cohesionless reinforced soil: q/σo versus ϕ. 

6.2 Strip Footing on Cohesive-Frictional Reinforced Soil

The exact bearing capacity of strip footing on weightless unreinforced soil is given by

qe=c[eπtanϕtan2(π4+ϕ2)1]cotϕ (39)

as demonstrated by Prandtl (1920). The advantage of having horizontal reinforcement can be estimated by analyzing numerically the change in (39) due to the reinforcement inclusion. The same mesh of Figure 4b is adopted and the mechanical boundary conditions are taken as

ts=0on the left(symmetric)edgetn=0 ts=0on the upper edge without strip footingtn=q ts=0on the upper edge with strip footing. (40)

The weight function component wn is imposed to be null at nodes on the left edge, and wn and ws are made to be null on the far right and bottom edges as before.

Table 2 presents the ratio qr/qe between the bearing capacities of strip footing on reinforced and unreinforced soil for different σo/c and angles of internal friction. It is assumed perfectly rough soil-reinforcement interface (ci = c, , ϕi = ϕ). Although the SOCP and LP results differ from each other by a small amount, the difference of iterations and CPU time between both optimization techniques is remarkably large. Results for σo/c = 0 illustrate the accuracy of our predictions for the bearing capacity of strip footing on unreinforced soil.

Table 2: Ratio qr/qe between the bearing capacities of strip footing on cohesive-frictional soil. 

σo/c ϕ = 10º ϕ = 20º ϕ = 30º
qr/qe iter CPU qr/qe iter CPU qr/qe iter CPU
0.0 1.0002† 20 0.67 0.9989 22 0.70 0.9789 21 0.67
0.9980‡ 2683 5 0.9959 2171 3 0.9742 2489 4
0.5 1.0956 23 0.78 1.0896 28 0.78 1.0572 26 0.75
1.0935 3192 5 1.0865 2090 3 1.0521 5394 8
1.0 1.1908 23 0.70 1.1802 27 0.78 1.1350 24 0.72
1.1887 4268 7 1.1770 2573 3 1.1301 4888 8
1.5 1.2859 23 0.72 1.2707 26 0.75 1.2125 24 0.72
1.2838 4248 7 1.2674 2647 4 1.2073 6519 11
2.0 1.3808 22 0.69 1.3611 25 0.73 1.2896 23 0.73
1.3788 3380 6 1.3577 2587 4 1.2842 4585 8

iter: number of iterations CPU: CPU time in seconds SOCP LP

Figure 6 shows the bearing capacity ratio qr/qu obtained with SOCP, where qu refers to the unreinforced scenario. For a given ϕ, the advantage of the reinforcement increases with σo/c. Indeed, the ratio qr/qu rises almost linearly with σo/c, where σo can be augmented by either increasing the reinforcement yield strength σyield or the reinforcement volumetric fraction d/h. We can observe that the benefit from horizontal reinforcement is less for soils with higher angle of internal friction.

Figure 6: Strip footing on cohesive-frictional soil: qr/qu versus σo/c for different angles of internal friction ϕ. 

6.3 Cohesionless Reinforced Earth Wall

To deal with the reinforced earth wall, the dimensions depicted in Figure 7a are considered sufficiently large to represent adequately the semi-infinite domain and the mesh adopted in the analysis is shown in Figure 7b. In this example, f1, present in (30), represents the vertical body force and f2 = 0, also present in (30), accounts for the null loads acting on Γt. The mechanical boundary conditions are

tn=0 ts=0on the three free edges. (41)

Figure 7: Earth wall: (a) adopted domain with horizontal reinforcement; (b) unstructured mesh with 4147 elements. 

On the far left, right and bottom edges, where tn and ts are supposed to be unknown, wn and ws are set equal to zero.

Table 3 condenses our results for horizontal reinforcement under perfectly rough soil-reinforcement interface (ci = c = 0, ϕi = ϕ). Note again the small difference between SOCP and LP predictions, but a much higher computational cost of an LP solution. Figure 8 shows that our results obtained with SOCP are bracketed by the lower bounds of Yu and Sloan (1997) and the upper bounds

γHupperσo=2tan2(π4+ϕ2) (42)

of Sawicki and Lesniewska (1987).

Table 3: Critical height of cohesionless reinforced earth wall. 

γH/σo iter CPU γH/σo iter CPU
10° 2.0428 24 2.18 2.0415 222199 1223
15° 2.6837 25 1.90 2.6814 348539 1792
20° 3.4463 23 2.04 3.4435 415115 2796
25° 4.3886 24 3.63 4.3832 428370 2439
30° 5.5307 25 2.00 5.5231 469423 2422
35° 6.9380 25 1.90 6.9269 541160 2836

iter: number of iterations CPU: CPU time in seconds

Figure 8: Cohesionless reinforced earth wall: γH/σo versus ϕ. 

6.4 Cohesive-Frictional Reinforced Earth Wall

To estimate the advantage of having horizontal reinforcement in cohesive-frictional earth walls, analyses are performed for cases with and without reinforcement. The same mesh and boundary conditions of the previous example are considered. Assuming perfectly rough soil-reinforcement interface (ci = c, ϕi = ϕ), Table 4 presents the ratio Hr/Hc between the critical heights of reinforced and unreinforced earth wall for different σo/c and angles of internal friction. The values Hc are elaborate upper bounds computed by Chen (1975) for the unreinforced soil using a logarithmic spiral surface discontinuity mechanism passing through the vertical slope toe. Although the difference between the SOCP and LP solutions is small, the computational cost of one technique differs excessively from the other. Results for σo/c = 0 illustrate that our solutions and Chen upper bounds are in close agreement.

Table 4: Ratio Hr/Hc between the critical heights of cohesive-frictional earth wall. 

σo/c ϕ = 10º ϕ = 20º ϕ = 30º
Hr/Hc iter CPU Hr/Hc iter CPU Hr/Hc iter CPU
0.0 0.9938† 23 2.48 1.0004 23 2.40 1.0037 23 2.50
0.9924‡ 231096 1283 0.9990 242228 1279 1.0017 264709 1401
0.5 1.2734 25 3.45 1.3440 25 3.42 1.4284 27 3.26
1.2722 319168 1636 1.3423 403088 2017 1.4259 490291 2449
1.0 1.5477 26 3.63 1.6840 25 3.46 1.8508 27 3.31
1.5462 300305 1559 1.6820 427462 2157 1.8480 515209 2613
1.5 1.8163 24 3.74 2.0215 24 3.43 2.2715 27 3.17
1.8149 296787 1474 2.0192 435428 2207 2.2681 547401 2712
2.0 2.0786 24 3.26 2.3566 23 2.98 2.6909 24 5.51
2.0769 319057 1561 2.3537 461221 2343 2.6871 533469 2637

iter: number of iterations CPU: CPU time in seconds SOCP LP

Figure 9 depicts the ratio Hr/Hu predicted with SOCP, where Hu refers to the critical height of unreinforced soil. The ratio increases almost linearly with σo/c for a given ϕ, similar to what is observed in Figure 6 for strip footing. However, the benefit of reinforcement inclusion in earth wall is greater for soils with higher angle of internal friction and shows to be more dependent on that angle.

Figure 9: Cohesive-frictional earth wall: Hr/Hu versus σo/c for different angles of internal friction ϕ. 

6.5 Cohesionless Reinforced Earth Wall Under Surcharge

Suppose the earth wall treated in the third example is now replaced by one which is weightless and subjected to a vertical load p (force per unit area) applied to the upper edge. The load is uniformly distributed over a length H measured from the vertical edge. Sawicki and Lesniewska (1987) managed to identify

peσo=tan2(π4+ϕ2) (43)

as the exact solution for the collapse load of this problem. Our analysis is carried out using the mesh of Figure 7b and the mechanical boundary conditions

tn=p ts=0on the upper edge with surchargetn=0 ts=0on the upper edge without surcharge and on the two free edges. (44)

The weight function components wn and ws are made to be null on the far left, right and bottom edges as before.

Table 5 outlines our results normalized to the exact solution (43), where excellent accuracy is observed for all angles of internal friction considered. Any attempt to solve the problem as LP has faced significant slowness in the simplex phase (>3600s).

Table 5: Bearing capacity of cohesionless reinforced earth wall under surcharge. 

ϕ 20° 25° 30° 35° 40° 45°
p/pe 0.9986 1.0000 1.0000 1.0000 1.0000 1.0000
iter 32 27 28 27 29 28
CPU 2.87 2.51 2.78 2.73 2.65 2.76

iter: number of iterations CPU: CPU time in seconds


The paper proposes a numerical approach for computing static limit loads of reinforced soil structures by combining discretization into simple triangular elements and SOCP. Examples illustrate that the above scheme provides excellent results and the fact of not satisfying the equilibrium equations and the mechanical boundary conditions rigorously is far from being a severe handicap for the developed element. For comparison, the examples are also treated as LP after a piecewise linearization of the reinforced soil yield condition. All the optimization problems idealized either as SOCP or LP are of large scale in the sense that they involve more than 10000 variables and constraints. The SOCP problems are solved by MOSEK and the LP problems are solved by FICO® Xpress Optimization Suite, both of which are state of the art optimization packages. It is clear from the results that the solution of large scale SOCP problems with the MOSEK interior-point optimizer is highly efficient and fast when compared with the combination of Newton barrier and simplex methods at FICO® Xpress Optimization Suite.


Anderheggen, E., Knöpfel, H. (1972). Finite element limit analysis using linear programming. International Journal of Solids and Structures 8(12):1413-1431. [ Links ]

Andersen, E.D., Roos, C., Terlaky, T. (2003). On implementing a primal-dual interior-point method for conic quadratic optimization. Mathematical Programming, Series B, 95(2):249-277. [ Links ]

Bleyer, J., de Buhan, P. (2014). Lower bound static approach for the yield design of thick plates. International Journal for Numerical Methods in Engineering 100(11):814-833. [ Links ]

Bottero, A., Negre, R., Pastor, J., Turgeman, S. (1980). Finite element method and limit analysis for soil mechanics problems. Computer methods in Applied Mechanics and Engineering 22(1):131-149. [ Links ]

Boyd, S., Vandenberghe, L. (2009). Convex Optimization. Cambridge University Press. [ Links ]

Chapuis, R.P. (1980). La compression triaxiale des sols armés. Canadian Geotechnical Journal 17(2):153-164. [ Links ]

Chen, W.F. (1975). Limit Analysis and Soil Plasticity, Elsevier (Amsterdam). [ Links ]

Ciria, H., Peraire, J., Bonet, J. (2008). Mesh adaptive computation of upper and lower bounds in limit analysis. International Journal for Numerical Methods in Engineering 75(8):899-944. [ Links ]

Davis, E.H., Booker, J.R. (1971). The bearing capacity of strip footings from the standpoint of plasticity theory. In Proceedings of the first Australian-New Zealand Conference on Geomechanics:276-282 (Melbourne). [ Links ]

de Buhan, P., Mangiavacchi, R., Nova, R., Pellegrini, G., Saleçon, J. (1989). Yield design of reinforced earth walls by a homogenization method. Géotechnique 39(2):189-201. [ Links ]

de Buhan, P., Siad, L. (1989). Influence of a soil-strip interface failure condition on the yield strength of reinforced earth. Computers and Geotechnics 7(1-2):3-18. [ Links ]

Drucker, D.C., Prager, W., Greenberg, H.J. (1952). Extended limit design theorems for continuous media. Quarterly Journal of Applied Mathematics 9(4):381-389. [ Links ]

Drud, A.S. (1996). CONOPT: A system for large scale nonlinear optimization, Reference manual for CONOPT subroutine library. ARKI Consulting and Development A/S (Bagsvaerd). [ Links ]

Fair Isaac Corporation (2013). FICO® Xpress-Optimizer reference manual. [ Links ]

Krabbenhoft, K., Damkilde, L. (2003). A general non-linear optimization algorithm for lower bound limit analysis. International Journal for Numerical Methods in Engineering 56(2):165-184. [ Links ]

Kulczykowski, M. (1985). Analysis of bearing capacity of reinforced subsoil loaded by a footing, Ph.D. Thesis (in Polish), Institute of Hydro-Engineering, Poland. [ Links ]

Le, C.V., Nguyen-Xuan, H., Nguyen-Dang, H. (2010). Upper and lower bound limit analysis of plates using FEM and second-order cone programming. Computers & Structures 88(1-2):65-73. [ Links ]

Lobo, M.S., Vandenberghe, L., Boyd, S., Lebret, H. (1998). Applications of second-order cone programming. Linear Algebra and its Applications 284(1-3):193-228. [ Links ]

Löfberg, J. (2004). YALMIP: A toolbox for modeling and optimization in MATLAB. In Proceedings of the IEEE CSS International Symposium on Computer Aided Control System Design:284-289 (Taipei). [ Links ]

Long, N.T., Guégan, Y., Legeay, G. (1972). Étude de la terre armée à l’appareil triaxial. Rapport de Recherches 17, Laboratoire Central des Ponts et Chaussées (Paris). [ Links ]

Lyamin, A.V., Sloan S.W. (2002). Lower bound limit analysis using non-linear programming. International Journal for Numerical Methods in Engineering 55(5):573-611. [ Links ]

Lysmer, J. (1970). Limit analysis of plane problems in soil mechanics. Journal of the Soil Mechanics and Foundations Division 96(SM4):1311-1334. [ Links ]

Makrodimopoulos, A., Martin, C.M. (2006). Lower bound limit analysis of cohesive-frictional materials using second-order cone programming. International Journal for Numerical Methods in Engineering 66(4):604-634. [ Links ]

MOSEK ApS (2016). Modeling cookbook. Available from [ Links ]

Muñoz, J.J., Bonet, J., Huerta, A., Peraire, J. (2009). Upper and lower bounds in limit analysis: adaptive meshing strategies and discontinuous loading. International Journal for Numerical Methods in Engineering 77(4):471-501. [ Links ]

Murtagh, B.A., Saunders, M.A. (2003). MINOS 5.51 user’s guide. Technical Report Sol 83-20R, Department of Management Science and Engineering, Stanford University (Stanford). [ Links ]

Nguyen-Thoi, T., Phung-Van, P., Nguyen-Thoi, M.H., Dang-Trung, H. (2015). An upper-bound limit analysis of Mindlin plates using CS-DSG3 method and second-order cone programming. Journal of Computational and Applied Mathematics 281:32-48. [ Links ]

Pastor, J. (1978). Analyse limite: détermination numérique de solutions statistiques complètes. Application au talus vertical. Journal de Mecánique Appliquée 2(2):167-196. [ Links ]

Pian, T.H.H., Wu, C.C. (2006). Hybrid and Incompatible Finite Element Methods, Chapman (Boca Raton). [ Links ]

Prandtl, L. (1920). Über die härte plastischer körper. Nachr K Ges Wiss, Göttingen, Mathematish-Physikalische Klasse:74-85. [ Links ]

Sawicki, A. (1983). Plastic limit behavior of reinforced earth. Journal of Geotechnical Engineering 109(7):1000-1005. [ Links ]

Sawicki, A., Lesniewska, D. (1987). Failure modes and bearing capacity of reinforced soil retaining walls. Geotextiles and Geomembranes 5(1):29-44. [ Links ]

Sawicki, A., Lesniewska, D. (1989). Limit analysis of cohesive slopes reinforced with geotextiles. Computers and Geotechnics 7(1-2):53-66. [ Links ]

Sturm, J.F. (1999). Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software 11-12:625-653. [ Links ]

Terzaghi, K. (1943). Theoretical Soil Mechanics, John Wiley Sons (New York). [ Links ]

Tütüncü, R.H., Toh, K.C., Todd, M.J. (2003). Solving semidefinite-quadratic-linear programs using SDPT3. Mathematical Programming, Series B, 95:189-217. [ Links ]

Washizu, K. (1982). Variational Methods in Elasticity and Plasticity, 3rd edn., Pergamon Press (Oxford). [ Links ]

Yu, H.S., Sloan, S.W. (1997). Finite element limit analysis of reinforced soils. Computers & Structures 63(3):567-577. [ Links ]

Yu, X., Tin-Loi, F. (2006). A simple mixed finite element for static limit analysis. Computers & Structures 84(29-30):1906-1917. [ Links ]

Received: February 05, 2017; Revised: July 21, 2017; Accepted: September 19, 2017

*Corresponding author

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License