Lower Bound Limit Analysis of Anisotropic Soils

Chunguang Li, Cuihua Li, Cong Sun and Hong Zheng

Nomenclature

1 Introduction

The question of slope stability is a well-known classical problem in soil mechanics, and it has been studied by numerous authors using Limit equilibrium methods [Cheng,Lansivaara and Wei (2007); Duncan (1996); Huang (2013); Lam and Fredlund (1993);Stolle and Guo (2008); Yu, Salgado, Sloan et al. (1998)], limit analysis [Bandini (2003);Chen, Yin, Li et al. (2003); Drescher and Christopoulos (1988); Han, Chen, Xia et al.(2014); Kim (1998); Li, Lyamin, Merifield et al. (2009); Tschuchnigg, Schweiger, Sloan et al. (2015); Yu, Salgado, Sloan et al. (1998)], finite element strength reduction method[Griffiths and Lane (1999); Griffiths and Marquez (2007); Tschuchnigg, Schweiger,Sloan et al. (2015); Zheng, Liu and Li (2005)], slip line method [Cheng (2003); Chenot,Felgeres, Lavarenne et al. (1978)] and variational method [Baker (2005); Baker and Garber (1977); Castillo and Luceno (1982)]. The first two methods are classical numerical methods and they are currently used widely in in practical engineering.However, almost all the references are focused on isotropic materials.

It is well known that natural soils and sedimentary rocks are usually deposited and consolidated under one-dimensional conditions, and hence most naturally occurring clays are inherently anisotropic. In 1966 Bishop [Bishop (1966)] performed a lot soil tests and noticed that in-situ undrained strength of a soil is influenced by the anisotropy, and the behavior of soils. The anisotropy of the soil structure, are likely to play a more significant part in the behavior of clays. By triaxial and plane strain tests of clayey soils, it is shown that the anisotropy of undrained shear strength of over-consolidated clay is primarily caused by the direction-dependent (anisotropic) shear strength [Kurukulasuriya, Oda and KAzAMA (1999)]. The anisotropy influence not only the strength of soils, but also on other properties of soils, such as CO2emission [La Scala Jr, Panosso, Pereira et al.(2009)], Hydraulic Conductivity [Soracco, Lozano, Sarli et al. (2010)], etc. Recently,using the Small-Strain Hollow Cylinder Apparatus, Yang had studied the drained anisotropic behavior of sand under generalized stress conditions, and the inherently anisotropic behavior of sands is clearly illustrated in his studies.

It has been suggested that the variation of soil cohesion with direction due to inherent anisotropy is much more significant than the anisotropy of material friction angle [Arthur and Menzies (1972)].

During the past several decades, the variation of shear strength with direction has been taken into account in stability analysis by a number of researchers using different methods. The methods include the limit equilibrium method, finite element method, finite difference method, and the limit analysis method. Aghajani et al. [Aghajani, Salehzadeh and Shahnazari (2015)] investigated the effect of anisotropy of shear strength parameter on the stability of a sandy slope by the limit equilibrium method. Al-Karniand Al-Shamrani using the method of slices, also investigated the influence of cohesion anisotropy on the stability of slopes in homogeneous soils. Since anisotropy can be considered in most limit equilibrium software, such as SLOPE/W [Krahn (2004)], the slope stability problem with anisotropic materials can be easily solved by commercial software. Based on Ubiquitous Joint Rock Mass (UJRM) model in FLAC3D(Itasca Consulting Group, 2012), various slope stability problems have been analyzed [Li, Dai,Li et al. (2011); Sainsbury and Sainsbury (2013); Sainsbury, Sainsbury and Sweeney(2016); Soren, Bodi and Sen (2014)]. Similarly, friction anisotropy also is taken into account by a lot of finite element software, such as ABAQUS, PLAXIS, etc. The stability of slope with fiber-reinforced soils was analyzed with limit analysis by Michalowski, he employed a fiber distribution function to characterize the fiber orientation anisotropy, and demonstrated the application of the kinematic approach of limit analysis to the anisotropic frictional materials [Michalowski (2008)]. Also using upper bound technique of limit analysis, Chen et al. [Chen, Snitbhan and Fang (1975)] established an expression under the condition of variation of cohesion with direction in undrained materials; Using an modified anisotropic yield criterion, Yu and Sloan [Yu and Sloan (1994)] generalized the conventional isotropic Mohr-Coulomb yield criterion to include the effect of variation of cohesion with direction in cohesive-frictional materials. Considering the failure mechanism as a series of rigid blocks, Yang et al. [Yang and Du (2016)] analyzed influence of anisotropy on the bearing capacity of surface strip footings.

It has been found that neglecting the anisotropy feature may lead to the overestimation of the safety factor for geotechnical earthworks and slope stability and thus lead to an unsafe design [Chen, Snitbhan and Fang (1975); Yu and Sloan (1994)]. Because the shear strength varies with direction, it is more realistic to use values of strength appropriate to each direction. However, very little work has been done on the effects of simultaneous variation of cohesion and internal friction angle.

The lower bound theorem is attractive because it gives a conservative result, and this paper aims to develop a general numerical method which can be used to calculate rigorous lower bound solutions for soils whose cohesion and internal friction varies with direction. The modified anisotropic yield surface is then developed to formulate the numerical lower bound limit analysis. The computational results obtained from the formulations were compared with those available in literature.

2 Lower bound theorems of limit analysis

The lower bound theorem provides a safe estimate of the collapse loads for a rigid plastic solid.

Consider a rigid plastic solid, subjected to some distribution of tractionsitand body forcesib. Collapse occurs under the loading,iLbλ, where the factorLλis in fact the safety factor.

To estimateLλ, we will denote the guess for the stress distribution in the solid byijσ. The stress distribution must

1.Satisfy the stress boundary conditions

in whichitis the prescribed traction.

2.Satisfy the equations of equilibrium

within the solid, whereib is the body force per unit volume.

3.Must not violate the yield criterion anywhere within the solid,

The lower bound theorem states that a statically admissible collapse load factor is always less than or equal to the exact collapse load factor

With these three sets of constraints, equilibrium equations, boundary conditions and yield criteria inequalities, the lower bound load optimization problem can be written as

maximizeλ

3 Brief review of lower limit analysis

To obtain a rigorous lower bound collapse load for a discretized continuum, three-node triangles are usually used to construct a piecewise continuous stress field in which statically admissible discontinuities are permitted across element boundaries. This approach was first applied to plane strain problems by Lysmer [Lysmer (1970)],following this early work, Sloan [Sloan (1988)] introduced finite element and mathematical programming formulation that permit large two-dimensional problems to be solved efficiently on a standard personal computer.

Fig. 1 shows the three-noded triangular elements for the lower bound formulation. Each node is associated with three unknown stresses as the nodal variables.

Thus, the stresses vary linearly throughout an element according to:

Where Ni, are linear shape functions.

Figure 1: Three-node linear stress triangle for lower bound limit analysis

3.1 Element equilibrium

Combination of Eqs. (2) and (6) leads to the following matrix form of element equilibrium equations

Thus, the equilibrium condition for each triangular element generates two equality constraints on the nodal stresses.

3.2 Boundary conditions

Substituting Eq. (6) into (1) gives

Thus each boundary edge generates four equality constraints on the nodal stresses.

3.3 Stress compatibility conditions

In order to permit statically admissible discontinuities, it is necessary to enforce the normal and shear stresses must be continuous across an admissible stress discontinuity at the edges of adjacent triangles. As shown in Fig. 2, element a and b share the side defined by nodal pairs (1, 2) and (3, 4), the stress compatibility condition requires

Figure 2: Stress discontinuity between adjacent triangles

3.4 Yield conditions

The Mohr-Coulomb yield criterion in the plain strain condition is stated as:

in which tensile stresses are taken as positive.

This yield function can be linearized by an inscribed polygon of p sides, as shown in Fig.3, which guarantee rigorous lower bound solutions. The new linearized yield function can be written as:

Figure 3: Linearization of Mohr-Coulomb yield function

3.5 Assembly of complete optimization problem

After the assembly of all the constraints equations and objective function, the complete problem formulation can then be written as:

4 Linearization of the mohr-coulomb yield criterion inspace

The Mohr-Coulomb yield criterion may be written as

Where τ and σnare the shear stress and the normal stress on a plane with normal unit vector n, c is cohesion and f=tanφ is friction coefficient and φ is angle of internal friction.As stated by the Mohr-Coulomb yield criterion, failure from shear will occur when the shear stress on a plane reaches a value given by Eq. (14).

Assume that Eq. (14) holds true for any direction at point P, then shear failure will not take place at this point because on any plane passing through point P, the available shear strength is greater than the maximum shearing stress.

The Mohr's circle is approximated by an interior polygon with p vertices starting from (σx,τxy), as shown in Fig. 4. The normal unit vector to the plane corresponding to the kth vertices is given by Mohr's circle is approximated by

where the plane corresponding to the kth vertices makes an angle of kπ/p counterclockwise to the plane perpendicular to x axis.

Adding more vertices, we can get the polygon closer to the Mohr-Coulomb circle and make the lower bound better. Generally speaking, it can be found that choosing p greater than 24 gives a sufficiently accurate approximation to the Mohr-Coulomb yield function.Substituting Eqs. (15-17) into (14) gives the linearized yield function as

It is worth noting that because Mohr’s circle and Coulomb failure line is symmetric aboutaxis, only the positive part ofare considered for Eq. (14).

And we can see that Eq. (18) is a linear equation ofand.

Figure 4: Discretization of Mohr’s circle in τ-σ space

5 Method to ensure the lower bound

Since the number of vertices used to approximate the Mohr circle is finite, not all of the stress states can be ensured to satisfy the Coulomb yield criterion, as is shown in Fig. 5, if we let the vertices on the Mohr circle satisfy Eq. (5), in the worst case, there are stresses on some slip directions (the bolded line on the Mohr circle) will violate the Mohr-Coulomb yield criterion.

Figure 5: The region lies outside the yield surface

Figure 6: Pseudo cohesion c'

To make the solution be a rigorous lower bound, considering circle and the dashed line which is parallel to the Mohr-Coulomb failure line and passes through two adjacent discretized points, as is shown in Fig. 6, when the Mohr's circle of stress touches the Coulomb failure line, the radius of the Mohr's circle is

and the distance from the center of Mohr's circle to the dashed line is

The dashed line that intersects the shear stress axis, gives a cohesion intercept, c'.

If we denote the intercept of the dashed line on the shear axis by c', and call it pseudo cohesion, we have

Using Eqs. (21) and (22), pseudo cohesion c' can be expressed as

Replacing the cohesion c in Eq. (18) with pseudo cohesion c' yields

Eq. (24) ensures that the solution obtained is a rigorous lower bound on the true collapse load.

6 Lower limit analysis for anisotropic materials

For a given anisotropic Mohr-Coulomb material, the cohesion strength c and the internal friction coefficient f can be given by

where θ represents the angle between the direction of a plane where the cohesion c or is internal coefficient f measured and the horizontal direction.

Figure 7: Variation of strength with directions in vertical-horizontal plane

For example, based on the earlier studies [Casagrande and Carillo (1944); Lo (1966)], the cohesion, c, with the major principal stress inclined at an angle θ with the horizontal may be written in the form:

in which chand cv, are the cohesion strengths in the horizontal and vertical directions,respectively.

A simple approach to computing rigorous bounds for anisotropic materials is to just replace the cohesion strength c and the internal friction coefficient f in Eq. (26) with

The advantage of this method is that it is suitable for any type of anisotropic Mohr-Coulomb material, not just for Eq. (27).

7 Numerical examples

7.1 The ACADS test example EX1(c)

The first example, named as Ex1(c), is a study from ACADS (Association for Computer Aided Design, Australia) of a non-homogeneous slope with three significantly different isotropic material parameters, whose geometry is shown in Fig. 8. The property parameters of the problem are shown in Tab. 1.

The slope is meshed uniformly with global element size specified as 1.0 m，as is shown in Fig. 9. Tab. 2 and Fig. 10 shows the comparison of safety factor obtained using c and c'for the nonhomogeneous slope. From Tab. 2 we see that the safety factors using pseudo cohesion c' have relatively good agreement with Sloan’s method [Sloan (1988)], and both the methods using c and c' will converge to a fixed value when p goes to infinity, their only difference between the two methods is that the curve of safety factor-p using c' is bounded from above (has an upper bound), however, the curve using c is bounded from below (has an lower bound).

Compared with the recommended result Fs=1.390, the relative error is 3.96%, a better result will be achieved if the geometry is meshed with more elements.

Figure 8: Slope example EX1(c) from ACADS. (unit: m)

Figure 9: Mesh of example EX1(c)

Table 1: The material property of example EX1(c)

Figure 10: A comparison of strength reduction factor using c and c' for example EX1(c)

7.2 A smooth trapdoor in anisotropic clay

The geometry of the problem is shown in Fig. 11, the stability factor for the trapdoor can be shown to be N=(γH+σs－σt)/ch, where γ is the unit weight of the soil [Sloan, Assadi and Purushothaman (1990)].

The lower bound mesh used to model a trapdoor with H/B=5 [Yu and Sloan (1994)]. Due to symmetry, only one half of the problem needs to be considered and the mesh is as shown in Fig. 12. To compute the stability number for a fixed H/B value, the surcharge and soil unit weight are set to zero [Yu and Sloan (1994)].

A summary of the stability lower bounds for various ratios of cv/ch, using an 18-sided approximation of the yield criterion, is shown in Tab. 3 and Fig. 13. It can be seen that the lower bounds agree well with that obtained by Yu’s method.

Figure 11: Trapdoor problem

Figure 12: Mesh used for analysis of the trapdoor problem

Table 3: Safety factor for trapdoor problem

Figure 13: Stability factor against the value of cv/ch

7.3 Anisotropic slope

The anisotropic slope shown in Fig. 14 has the following soil properties:

The slope is inclined at an angle of 26.57° (2:1) to the horizontal and the boundary conditions are given as vertical rollers on the left boundary and full fixity at the base[Griffiths and Lane (1999)].

With the uniform element size specified as 0.06 H, the mesh of the slope is shown in Fig.15.The lower bound of safety factor is 1.364, which closely agrees with the safety factor 1.380 given for the same problem by the charts of Bishop and Morgenstern [Bishop and Morgenstern (1960)].

Suppose that chand fhsatisfy Eq. (35), let’s consider the following three cases:

(a) Varying cv/chand fv/fh=1.

(b) Varying fv/fhand cv/ch=1.

(c) Varying cv/ch= fv/fh.

Fig. 16 shows the stability factors of the anisotropic slope for the three cases, and the stability factor solutions are listed in Tab. 4. For case (a), not surprisingly, the stability factor increases with increasing cv/chwhen fv/fh=1. Similar behavior to the above can also be seen for case (b) and (c). From the figure it can be seen that the effect on the stability factor of fv/fhis here more pronounced than that of cv/ch.

Figure 14: Geometry and boundary conditions of anisotropic slope

Figure 15: Mesh for anisotropic slope

Table 4: Safety factor for anisotropic slope

Figure 16: Stability factor against the value of cv/ch, fv/fh

8 Conclusions

This paper presents a general finite element formulation of the lower bound theorems for soil whose cohesion and friction coefficient varies with direction. To develop the numerical formulations for anisotropic soil, the Mohr’s circle is discretized inspace and pseudo cohesion is defined to ensure rigorous lower bound solutions. Several numerical examples given in the paper illustrate that the proposed numerical procedure can be used to compute rigorous bound solutions of anisotropic soils.

Aghajani, H. F.; Salehzadeh, H.; Shahnazari, H. (2015): Stability analysis of sandy slope considering anisotropy effect in friction angle. Sadhana, vol. 40, pp. 1955-1974.

Arthur, J.; Menzies, B. (1972): Inherent anisotropy in a sand. Geotechnique, vol. 22, pp.115-128.

Baker, R. (2005): Variational Slope Stability Analysis of Materials with Nonlinear Failure Criterion. Electronic Journal of Geotechnical Engineering, vol. 10(Bundle A).

Baker, R.; Garber, M. (1977): Variational approach to slope stability. TECHNION-ISRAEL INST OF TECH HAIFA.

Bandini, P. (2003): Numerical limit analysis for slope stability and bearing capacity calculations. Purdue University.

Bishop, A. (1966): The strength of soils as engineering materials. Geotechnique, vol. 16,pp. 91-130.

Bishop, A.; Morgenstern, N. (1960): Stability coefficients for earth slopes.Geotechnique, vol. 10, pp. 129-153.

Casagrande, A.; Carillo, N. (1944): Shear failure of anisotropic materials. Journal of Boston Society of Civil Engineers, vol. 31, pp. 74-81.

Castillo, E.; Luceno, A. (1982): A critical analysis of some variational methods in slope stability analysis. International Journal of Rock Mechanics & Mining Sciences &Geomechanics, vol. 6, pp. 195-209.

Chen, J.; Yin, J. H.; Lee, C. (2003): Upper bound limit analysis of slope stability using rigid finite elements and nonlinear programming. Canadian Geotechnical Journal, vol.40, pp. 742-752.

Chen, W. F.; Snitbhan, N.; Fang, H. Y. (1975): Stability of slopes in anisotropic,nonhomogeneous soils. Canadian Geotechnical Journal, vol. 12, pp. 146-152.

Cheng, Y. (2003): Seismic lateral earth pressure coefficients for c-φ soils by slip line method. Computers and Geotechnics, vol. 30, pp. 661-670.

Cheng, Y.; Lansivaara, T.; Wei, W. (2007): Two-dimensional slope stability analysis by limit equilibrium and strength reduction methods. Computers and Geotechnics, vol. 34,pp. 137-150.

Chenot, J.; Felgeres, L.; Lavarenne, B.; Salencon, J. (1978): A numerical application of the slip line field method to extrusion through conical dies. International Journal of Engineering Science, vol. 16, pp. 263-273.

Drescher, A.; Christopoulos, C. (1988): Limit analysis slope stability with nonlinear yield condition. International Journal for Numerical and Analytical Methods in Geomechanics, vol. 12, pp. 341-345.

Duncan, J. M. (1996): State of the art: Limit equilibrium and finite-element analysis of slopes. Journal of Geotechnical Engineering-Asce, vol. 122, pp. 577-596.

Griffiths, D.; Lane, P. (1999): Slope stability analysis by finite elements. Geotechnique,vol. 49, pp. 387-403.

Griffiths, D. V.; Marquez, R. M. (2007): Three-dimensional slope stability analysis by elasto-plastic finite elements. Geotechnique, vol. 57, pp. 537-546.

Han, C. Y.; Chen, J. J.; Xia, X. H.; Wang, J. H. (2014): Three-dimensional stability analysis of anisotropic and non-homogeneous slopes using limit analysis. Journal of Central South University, vol. 21, pp. 1142-1147.

Huang, Y. H. (2013): Slope stability analysis by the limit equilibrium method:fundamentals and methods. American Society of Civil Engineers (ASCE).

Kim, J. (1998): Limit analysis of soil slope stability using finite elements and linear programming. Purdue University.

Krahn, J. (2004): Stability modeling with Slope/W. An Engineering Methodology.Calgary, Canada, Geo-Slope/W international LTD.

Kurukulasuriya, L. C.; Oda, M.; KAzAMA, H. (1999): Anisotropy of undrained shear strength of an over-consolidated soil by triaxial and plane strain tests. Soils Found, vol.39, pp. 21-29.

La Scala Jr, N.; Panosso, A.; Pereira, G.; Gonzalez, A.; Miranda, J. (2009): Fractal dimension and anisotropy of soil CO2 emission in an agricultural field during fallow.International Agrophysics, vol. 23, pp. 353-358.

Lam, L.; Fredlund, D. (1993): A general limit equilibrium model for three-dimensional slope stability analysis. Canadian Geotechnical Journal, vol. 30, pp. 905-919.

Li, A.; Lyamin, A.; Merifield, R. (2009): Seismic rock slope stability charts based on limit analysis methods. Computers and Geotechnic, vol. 36, pp. 135-148.

Li, H.; Dai, F.; Li, W.; Xu, L.; Min, H. (2011): Stability assessment of a slope under a transformer substation using numerical modelling. Bulletin of Engineering Geology and the Environment, vol. 70, pp. 385-394.

Lo, K. (1966): Stability of slopes in anisotropic soils. Journal of Soil Mechanics &Foundations Divsion, vol. 92(Closure).

Lysmer, J. (1970): Limit analysis of plane problems in soil mechanics. Journal of Soil Mechanics & Foundations Division, vol. 96(SM4), pp. 1311-1334.

Michalowski, R. (2008): Limit analysis with anisotropic fibre-reinforced soil.Geotechnique, vol. 58, pp. 489-501.

Sainsbury, D.; Sainsbury, B. (2013): Three-dimensional analysis of pit slope stability in anisotropic rock masses. Slope Stability, vol. 2013, pp. 683-696.

Sainsbury, D.; Sainsbury, B.; Sweeney, E. (2016): Three-dimensional analysis of complex anisotropic slope instability at MMG’s Century Mine. Mining Technology, vol.125, pp. 212-225.

Sloan, S. W. (1988): lower bound limit analysis using finite-elements and linear-programming.International Journal for Numerical and Analytical Methods in Geomechanics, vol. 12,pp. 61-77.

Sloan, S. W.; Assadi, A.; Purushothaman, N. (1990): Undrained Stability of A Trapdoor.Geotechnique, vol. 40, pp. 45-62.

Soracco, C. G.; Lozano, L. A.; Sarli, G. O.; Gelati, P. R.; Filgueira, R. R. (2010):Anisotropy of saturated hydraulic conductivity in a soil under conservation and no-till treatments. Soil and Tillage Research, vol. 109, pp. 18-22.

Soren, K.; Budi, G.; Sen, P. (2014): Stability analysis of open pit slope by finite difference method. Journal of Engineering and Advanced Research Technology, vol. 3,pp. 326-334.

Stolle, D.; Guo, P. J. (2008): Limit equilibrium slope stability analysis using rigid finite elements. Canadian Geotechnical Journal, vol. 45, pp. 653-662.

Tschuchnigg, F.; Schweiger, H. F.; Sloan, S. W. (2015): Slope stability analysis by means of finite element limit analysis and finite element strength reduction techniques.Part I: Numerical studies considering non-associated plasticity. Computers and Geotechnics, vol. 70, pp. 169-177.

Yang, X. L.; Du, D. C. (2016): Upper bound analysis for bearing capacity of nonhomogeneous and anisotropic clay foundation. KSCE Journal of Civil Engineering,vol. 20, pp. 2702-2710.

Yu, H.; Sloan, S.; (1994): Limit analysis of anisotropic soils using finite elements and linear programming. Mechanics Research Communications, vol. 21, pp. 545-554.

Yu, H. S.; Salgado, R.; Sloan, S. W.; Kim, J. M. (1998): Limit analysis versus limit equilibrium for slope stability. Journal of Geotechnical and Geoenvironmental Engineering, vol. 124, pp. 1-11.

Zheng, H.; Liu, D.; Li, C. (2005): Slope stability analysis based on elasto-plastic finite element method. International Journal for Numerical Methods in Engineering, vol. 64,pp. 1871-1888.