1 INTRODUCTION

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.

2 FAILURE CONDITIONS

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

where
^{Yu and Sloan, 1997}).

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

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

into (2) yields

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

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}):

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

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

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

in the system xy.

3 STATIC THEOREM

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

throughout the domain Ω, where the differential operator

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

on the portion Γt of Γ on which the traction (stress vector) t = ⌊tx ty⌋^{
T
} is prescribed as being

contains the components of n.

3.3 Yield Criterion

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

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

where

4 FINITE ELEMENT FORMULATION

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:

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}).

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

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

In view of the divergence theorem, one writes

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

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

where

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

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,

with nodal values

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

where

and

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

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}).

5 THE SOCP PROBLEM

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

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
^{
r
} respectively.

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

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

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

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

sketched in Figure 3. Now, the problem can be treated as SOCP by just replacing

where the vectors

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}).

6 RESULTS

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 VAIO^{TM} 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.

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

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:

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

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.

ϕ | SOCP | LP | ||||||
---|---|---|---|---|---|---|---|---|

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 ϕ.

6.2 Strip Footing on Cohesive-Frictional Reinforced Soil

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

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

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.

σ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.

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

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

of ^{Sawicki and Lesniewska (1987}).

ϕ | SOCP | LP | ||||||
---|---|---|---|---|---|---|---|---|

γ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

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.

σ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.

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

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

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).

7 CONCLUSIONS

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.