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

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.


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 Latin American Journal of Solids and Structures 14 (2017) 2497-2517 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 andLesniewska, 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.

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 ( 1 d h  ), as shown in Figure 1.Three stress tensors are defined at every point in the homogenized continuum: tensor s s q s q t q q s s q s q t q q t s s q q t q q = + - where q represents the angle between the horizontal axis x and the reinforcement direction t , meas- ured counterclockwise.Substitution of (1) and relations s s q s q t q q s s q s q t q q t s s q q t q q The reinforcement inside the soil is supposed to act merely as tensile load carrying elements (de Buhan et al., 1989) 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 i c and i f 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 .

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

Equilibrium Equations
Let W be the region occupied by the reinforced soil structure subjected to the body force throughout the domain W , where the differential operator

Mechanical Boundary Conditions
Let G be the boundary of the domain W , with unit outward normal vector denoted by . In addition to (10), the macrostress field must also satisfy the mechanical boundary conditions on the portion t G of G on which the traction (stress vector) contains the components of n.
Latin American Journal of Solids and Structures 14 (2017) 2497-2517   Figure 2: Finite element e W with the unit normal vector n on its boundary e G .

Yield Criterion
To assure that the yield criterion is satisfied it is necessary to impose simultaneously As the function i F defined in ( 9) has a term in modulus, the second of inequalities ( 14) splits into two other conditions given by where

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 e W shown in Figure 2, the definition of the boundary should be extended to include the traction continuity on the interelement portion i G : ( ) ( ) where the superscripts " + " and " -" denote the two sides of i G .Expressions ( 10), ( 12) and ( 17) may now be viewed as the equilibrium within the element and on its boundary t G and i G , 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 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 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 The weight function is also taken to vary linearly, Latin American Journal of Solids and Structures 14 (2017) 2497-2517 with nodal values 1,2, 3 Substitution of ( 21) and ( 24) into (20) yields where x y y ds dx dy Because ( 26) holds for any arbitrary weight function, it follows that Elimination of the integral over u G by choosing w null over there introduces zero components of 1 w , 2 w or 3 w 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).
Latin American Journal of Solids and Structures 14 (2017) 2497-2517 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: 1 lf which is adjusted during the optimization by means of the load factor l , and 2 f which is kept constant.Vectors s and r s collect the nodal values of s and , r s 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 Figure 3: Second-order cone (34).
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 Latin American Journal of Solids and Structures 14 (2017) 2497-2517 into (7) to state 0 s F £ as the three-dimensional second-order cone sketched in Figure 3. Now, the problem can be treated as SOCP by just replacing 1 ( , ) r £ 0 s s g with where the vectors v and 3 v collect the nodal values of v and 3 v .The new constraints (35) stem from the evaluation of ( 33) and ( 34) at each node, and all the problem nonlinearity concentrates in the constraint 2 £ 0 h related to the cone definition.
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 1 £ 0 g replaced by ( 35); (b) solution by MOSEK using the interior-point optimizer developed for nonlinear convex optimization by Andersen et al. (2003).

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 s F 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.

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, 1 , f present in (30), denotes all the loads applied on t G and 2 , = 0 f 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, n t and s t are the traction components normal and tangential to the edges.The weight function components n w and s w are imposed to be null at nodes on the edges where the traction components n t and s t are unknown, respectively: In addition, n w and s w are made to be null on the far right and bottom edges because n t and s t 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 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 f =  the count of 10037 takes place in the simplex phase.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 s F .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 f. o q s versus f .

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  (40) The weight function component n w is imposed to be null at nodes on the left edge, and n w and s w are made to be null on the far right and bottom edges as before.
Table 2 presents the ratio r e q q between the bearing capacities of strip footing on reinforced and unreinforced soil for different o c s and angles of internal friction.It is assumed perfectly rough soil- 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 0 c s = illustrate the accuracy of our predictions for the bearing capacity of strip footing on unreinforced soil.
q q iter CPU r e q q iter CPU r e q q iter CPU 0.  q q between the bearing capacities of strip footing on cohesive-frictional soil.
Figure 6 shows the bearing capacity ratio r u q q obtained with SOCP, where u q refers to the unreinforced scenario.For a given f , the advantage of the reinforcement increases with o c s .Indeed, the ratio r u q q rises almost linearly with o c s , where o s can be augmented by either increasing the reinforcement yield strength yield s 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.q q versus o c s for different angles of internal friction f .

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.
Note again the small difference between SOCP and LP predic- tions, 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)  of Sawicki and Lesniewska (1987).

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

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)  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 n w and s w 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 ( 3600 > s).

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.
u G of the element boundary e G falls on the reinforced soil boundary with pre- scribed displacement.

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

Figure 5 :
Figure 5: Strip footing on cohesionless reinforced soil: is adopted and the mechanical boundary conditions are taken as edge with strip footing.

Figure 9 :
Figure 9: Cohesive-frictional earth wall: r u H H versus o c s for different angles of internal friction f .

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

Table 2 :
Ratio r e

Table 3 :
Critical height of cohesionless reinforced earth wall.
o H g s versus f .

Table 4
Chen (1975)s computed byChen (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

Table 4 :
Ratio r c H H between the critical heights of cohesive-frictional earth wall.
Figure 9 depicts the ratio r u H H predicted with SOCP, where u H refers to the critical height of unreinforced soil.The ratio increases almost linearly with o c s for a given f , similar to what is observed in Figure

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