Yongtao Yang, Wenan Wu, Hong Zheng

a State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences, Wuhan,430071,China

b University of Chinese Academy of Sciences, Beijing,100049, China

c Key Laboratory of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, Beijing,100124, China

Keywords:Enhanced limit method Ant colony algorithm Numerical manifold method (NMM)Slope stability

A B S T R A C T This study first reviews the numerical manifold method (NMM) which possesses some advantages over the traditional limit equilibrium methods (LEMs) in calculating the factors of safety (Fs) of the slopes.Then, with regard to a trial slip surface (TSS), associated stress fields reproduced by NMM as well as the enhanced limit equilibrium method are combined to compute Fs. In order to search for the potential critical slip surface (CSS), the MAX-MIN ant colony optimization algorithm (MMACOA), one of the best performing algorithms for some optimization problems, is adopted. Procedures to obtain Fs in conjunction with the potential CSS are described. Finally, the proposed numerical model and traditional methods are compared with stability analysis of three typical slopes. The numerical results show that Fs and CSSs of the slopes can be accurately calculated with the proposed model.

1. Introduction

Failure in rock/soil slopes frequently results in serious consequences to both properties and lives (Farias and Naylor, 1998).Hence, it is essential to investigate the stability of the slopes and associated risks.In the slope stability analysis,two issues need to be resolved:calculation of the factor of safety(Fs),and determination of the critical slip surface (CSS).

The factor of safety for simple slopes can be obtained easily by the traditional limit equilibrium methods (LEMs), which still prevails in the literature (e.g. Bishop, 1955; Morgenstern and Price,1965; Spencer,1967; Janbu,1973; Fredlund and Krahn,1977; Sun et al., 2017). However, the traditional LEMs have limits as they cannot reveal the failure mechanism of slopes. In addition, the stress-strain relationship and deformation behaviors of the slopes are not taken into account. Hence, the stress distribution modes and the shapes of the CSSs have to be assumed (Farias and Naylor,1998; Griffiths and Lane,1999).

Compared with the traditional LEMs, numerical methods (e.g.Zienkiewicz and Taylor, 2000; Yang et al., 2017; Rabczuk and Belytschko, 2004, 2007; Zhuang et al., 2012; Mohammadnejad and Khoei, 2013; Vu-Bac et al., 2016; Hamdia et al., 2017; Zhou et al., 2018; Anitescu et al., 2019; Chen et al., 2018; Yang et al.,2020a,b,c,d; Xu et al., 2020) which are capable of reproducing the stress and deformation fields of the problem domain may be considered as a better alternative for slope stability analysis.

The numerical manifold method(NMM)has received increasing attention among various numerical methods over the past two decades (Shi,1991; Chen et al., 2019; Yang et al., 2016a, b, 2019c).Unlike the finite element method (FEM), mathematical meshes with regular shapes which do not have to match the physical meshes can be adopted to discretize the problem domain.Furthermore, discontinuities in NMM can be captured in a natural manner without using the Heaviside functions (Mohammadnejad and Khoei, 2013). Hence, the NMM is adopted to cope with many complex problems(e.g.Jiang et al.,2010;Ning et al.,2011;Wu and Wong, 2012; Zheng et al., 2013; Zhao et al., 2014; Wei et al., 2017;Yang et al.,2018a,b,c;Zheng et al.,2019;Zheng and Yang,2017;Wu et al., 2020a,b,c; Wang et al., 2019; Yang and Zheng, 2016, 2017).

In the slope stability analysis, the strength reduction technique(SRT)(Zheng et al.,2005)which was initially used in FEM has also been incorporated into NMM. The strength reduction NMM(SRNMM)has been adopted by Yang et al.(2019a,2020e)to assess the stability of soil-rock-mixture (SRM) slopes. Their results indicated that SRNMM can not only accurately predict the factors of safety for both simple and SRM slopes, but also capture the main failure modes of the slopes.

Fig.1. Basic concepts in numerical manifold method (NMM).

Fig. 2. Details about basic concepts in the numerical manifold method presented in Fig.1:(a)The domain of MP1;(b)The domain of PP1 and its corresponding node (c) The domain of PP3 and its corresponding node (d) The domain of PP2 and its corresponding node (e)The domain of E2;(f)The four nodes corresponding to E2;(g) The domain of E3; and (h) The four nodes corresponding to E3.

Fig. 3. Sketch of the trial slip surface for a simple slope.

However, the computation cost to obtain Fsof the slopes using SRNMM is significant. In SRNMM, it is essential to conduct a number of elasto-plastic analyses by reducing cohesion c and angle of internal friction φ using a factor which will be considered as Fswhen the failure of the slope occurs. As the value of this factor increases, these elasto-plastic analyses become more timeconsuming. Similar conclusions also hold in the case of strength reduction FEM (SRFEM)(also called direct FEM by Naylor (1982)).

In the present work, a numerical model called ELNMM is built by incorporating the enhanced limit method (ELM) (Naylor,1982)into NMM. In ELNMM, the stress fields reproduced by NMM and ELM are combined to calculate Fsof the slopes. Compared with SRNMM, one significant advantage of ELNMM is that only one elastoplastic analysis using the actual cohesion and angle of internal friction is needed. The present paper will seek for the details about the implementation of ELM into NMM. Note that ELM has already been incorporated into FEM(Farias and Naylor,1998;Zheng et al., 2006). However, the incorporation procedures of ELM into NMM are much easier than those of ELM into FEM.

In ELNMM, a series of trial slip surfaces (TSSs) is needed to calculate Fs. In these TSSs, the one with the minimum Fsis considered as the CSS.To search for the CSS,many methods such as the pattern search schemes (Baker,1980), calculus-based method(Zolfaghari et al.,2005)and genetic algorithm(Li et al.,2010)can be used.

The ant colony algorithm is a bionic optimization algorithm,which simulates the collective intelligent behavior of natural ant colony (Dorigo and Stutzle, 2004). It has the characteristics of openness,robustness and parallelism,and is good for dealing with complex discrete optimization problems. This method is mostly adopted for combinatorial optimization problems.Since the search for CSSs of the slopes can be considered as a single path search,the ant colony optimization algorithm (ACOA) which copes with path search problems should be suitable for locating the CSSs of the slopes.Application of the ACOA in searching for CSSs of slopes can be found in Sun et al. (2014). To improve the performance of the traditional ACOA, a set of improvements has been proposed in recent years (e.g. Gao et al., 2012). In all these algorithms, a set of discrete points is arranged in advance, and then ant colony crawls on these points to find out the optimal path. Apart from these algorithms, Stutzle and Hoos (2000) proposed the MAX-MIN ACOA(MMACOA), which is currently among the best performing algorithms for the traveling salesman problem and the quadratic assignment problem. In the present work, the MMACOA will be employed to search for the CSS of the slope.

Fig. 4. Sketch for calculating local coordinates of point P(x, y).

Fig. 5. Distribution of discrete points in the slope.

2. Brief introduction to numerical manifold method

2.1. Basic concepts in numerical manifold method

In order to introduce the basic concepts of NMM, Figs.1 and 2 are taken for examples. Two cover systems (CSs), i.e. mathematical CS and physical CS,are adopted in NMM(Zheng and Xu,2014;Yang et al., 2019d).

The mathematical CS is the union of mathematical patches(MPs).Each MP is formed by a series of rectangular grids sharing a mutual node. For example, as shown in Fig. 1, MP1and MP2are formed by four rectangular grids. Compared with FEM, the mathematical CS in a regular shape in NMM to discretize the problem domain does not have to match the physical meshes.This merit of NMM will not only reduce the burden of pre-processing, but also ensure the accuracy. Note that regular meshes usually bring in smaller errors in computation than distorted meshes. However,distorted meshes are inevitable for problems with complex domains in FEM.

Fig. 6. A trial slip surface.

The physical CS is the union of physical patches (PPs). PPs are produced through slicing MPs with physical meshes (PMs, thick real lines in red plotted in Fig.1).Cutting MP1(Figs.1 and 2a)with PMs generates PP1(the shadow area plotted in Fig.2b);Cutting MP2(Fig.1)with PMs generates PP2(the shadow area plotted in Fig.2d)and PP3(the shadow area plotted in Fig.2c).Note that there must be exactly one NMM node attached to a PP. For instance, the NMM nodeis attached to PP1,and the NMM nodesandare attached to PP2and PP3,respectively.The location of the NMM node can either be inside (Fig. 2b and d) or outside the PP domain(Fig. 2c). In addition, two or even more NMM nodes can overlap with each other, for example,andare overlapping. However, for problems undergoing large deformations or rigid movements, the initially overlapping NMM nodes can also move away from each other.Note that in NMM,the degrees of freedom(DOFs)are attached on nodes.

The overlapped zone of four PPs is a manifold element(ME).For example, E1is the overlapped zone of PP4, PP5, PP6and PP7, as shown in Fig. 1. ME is the basic integration element in NMM. It should be noted that ME can be in any shape,such as a rectangle(E1plotted in Fig.1),or a trapezoid(E2plotted in Fig.2e,and E3plotted in Fig. 2g).

For an ME, the global approximation function is obtained by weighting the four related cover functions:

Fig.7. Flow chart for determination of the critical slip surface using the MAX-MIN ant colony optimization algorithm.

Fig. 8. (a) A sketch for the homogeneous slope and (b) the corresponding discretized model.

where wk（x）and uk（x）are the weight function and cover function,respectively.For E1,the four related cover functions are defined on PP4,PP5,PP6and PP7(Fig.1).For E2,the four related cover functions are defined on PP8, PP9, PP10and PP11(Fig. 2f). For E3, the four related cover functions are defined on PP12, PP13, PP14and PP15(Fig. 2h).

uh（x） can also be rewritten in the following form:

where

where d and N represent the vector of DOFs and matrix of shape function, respectively; and uiand virepresent the x- and ycomponent of displacement vector at node i, respectively.

2.2. System equations in matrix forms

The weak form of static elastoplastic solid problems can be expressed as

Table 1 Properties of material used in the first case.

Fig. 9. Displacement contours (in m) of the homogeneous slope obtained from the numerical manifold method: (a) x- and (b) y-displacement.

where Ω denotes the domain of problem, Depdenotes the elastoplastic matrix, ε denotes the strain vector, δu denotes the virtual displacement vector, b denotes the vector of body force,t denotes the designated traction acting on Γt.

The system equations in the discretized form can be obtained through substituting Eq. (2) into Eq. (5) as

where

3. Implementation of enhanced limit method into numerical manifold method

3.1. Brief introduction to enhanced limit method

A sketch for TSS denoted as AB of a simple slope is plotted in Fig.3.According to the ELM(Naylor,1982;Farias and Naylor,1998),Fsfor this TSS is calculated by

Fig.10. Stress contours(in Pa)of the homogeneous slope obtained from the numerical manifold method: (a) σx, (b) σy, and (c) τxy.

where τnrepresents the component of shear stress along the TSS,and τfrepresents the shear strength of the slope material.

The stress at an arbitrary point lying on AB is defined by the stress components including σx, σyand τxy. The normal stress component σnand the shear stress component τnalong the TSS can be transformed using the stress components, and expressed as

The shear strength τfat this point can be calculated based on the Mohr-Coulomb (M-C) criterion:

3.2. Stress component determination

To calculate the factor of safety Fsin Eq. (9), the stresses of the points on the TSS should be calculated first using the displacement function in NMM analysis. In this case, the location of the TSS should be known in advance.

In the present work,the stresses for a specified point P(x,y)are calculated through interpolating the nodal stresses in the grid to which the point belongs. Let σP= ［σx， σy， τxy］ represents the stresses at point P(x, y). The following interpolation will be used:

where

where T is the interpolating matrix, （ξ，η） represents the local coordinates of point P, andis the matrix including the nodal stresses of four PPs.

To obtain the stresses for a specified point P(x,y)using Eq.(13),three steps have to be performed as follows: (1) For a specified point defined in the Cartesian coordinates,the corresponding ME to which the point belongs should be identified; (2) The local coordinates for this point should be calculated so as to compute the interpolating matrix T（ξ，η）; and (3) The nodal stresses should be computed.

In the first step, the point method (Shi, 1988) is adopted to determine whether a point is within a given polygonal region or not. In the second step, it is very easy to calculate the local coordinates of the specified point in NMM. As illustrated in Fig. 4,since the rectangular grids are adopted to form the MPs in this study,the interpolating mesh for an ME is also in rectangular shape.The interpolating mesh for the manifold element E4is in rectangular shape, although E4is in trapezoid shape. Hence, the local coordinates of P(x, y) can be easily computed using the following equations:

where （xi，yi） represents the Cartesian coordinates of the four interpolating NMM nodes; and lxand lyrepresent the width and height of the interpolating mesh, respectively.

In the third step, the nodal stresses can be computed by averaging the stresses of Gaussian point, and expressed as

Fig.11. Locations of potential critical slip surfaces for the homogeneous slope obtained from different numerical methods:(a)Bishop,(b)Janbu,(c)Spencer,(d)Morgenstern-Price,and (e) Proposed method.

Table 2 Factors of safety obtained from different numerical methods for the first case.

where ne represents the number of manifold elements related to PPi, and Akrepresents the area of Ek(k = 1, 2, …, ne).

Fig. 12. (a) A sketch for the slope with two layers of geo-materials and (b) the corresponding discretized model (Yang et al., 2019b).

4. MAX-MIN Ant colony optimization algorithm

In this section,MMACOA(Stutzle and Hoos,2000)is introduced to locate the TSS with the minimum Fs,which will be considered as the CSS. Similar to the traditional ACOA, a set of discrete points should be deployed first, and then ant colony will crawl on these points to find out an optimal path.

For application of MMACOA to the slope stability analysis, we take Fig. 5 as an example. The slope body is discretized by a set of vertical lines, which are named as “stripe lines” herein. Over each stripe line, a set of discrete points are deployed.

The ant crawls firstly from a point on the START region to the next point on the neighboring stripe lines according to the probability presented in Eq.(19).Then,it keeps crawling until reaching apoint on the END region. The crawling trajectory of this ant is considered as a TSS.

Table 3 Properties of material used in the second case.

Fig. 13. Displacement contours (in Pa) of the slope with two layers of geo-materials obtained from the numerical manifold method: (a) x- and (b) y-displacement.

Let（r， i）represent the i-th discrete point on the r-th stripe line,（r+1， j） represent the j-th discrete point on the (r + 1)-th stripe line, ［（r， i）， （r+1， j）］ denote the path between （r， i） and （r + 1， j）,and m stand for the number of ants in the colony.In the process of movement,the transfer direction for the ant k is predicted using the intensity of pheromone and the heuristic value. Then, at the moment t,the probability for the ant k transferring from point（r， i）to point （r+1， j） can be calculated by

where Mr+1represents the number of discrete points on the(r+1)-th stripe line; τ［（r，i），（r+1，j）］（t） represents the pheromone intensity at path ［（r， i）， （r+1， j）］ at the moment t; η［（r，i），（r+1，j）］（t） represents the heuristic value at path ［（r， i）， （r + 1， j）］;and v and β are the specified parameters determining the relative importance with regard to the pheromone trail and the heuristic information,respectively.

Fig.14. Stress contours (in Pa) of the slope with two layers of geo-materials obtained from the numerical manifold method: (a) σx, (b) σy, and (c) τxy.

As the process of iteration continues, the pheromone intensity in the path of the ants decreases gradually.Let ρ（0 ≤ρ ≤1）denote the persistence of the pheromone intensity in the path, and 1- ρ represent the decay degree of the pheromone intensity.After all the m ants find their own complete paths in a cycle according to Eq.(19), m TSSs will be formed. Based on these m TSSs, a series of feasible solutions, i.e. factors of safety, can be obtained. Then the value of pheromone intensity over a path is adjusted by

Otherwise, following equation is obtained:

In order to avoid prematurity due to overemphasis of some paths during the pheromone update,the pheromone on each path is limited within ［τmin， τmax］, and the pheromone beyond or below this range is set to τmaxor τmin.Additionally,for the purpose of improving the searching efficiency in conjunction with obtaining a reasonable TSS, following conditions are assumed as (Sun et al.,2014; Gao,2016):

(1) The horizontal coordinates of the searching paths formed by a set of discrete points should be in increase or decrease order.

(2) The entrance and exit points of the TSS should be located on the slope surface.

(3) The range of horizontal coordinates of the searching paths should be within the horizontal border of the slope, while the range of vertical coordinates of the searching paths should be within the vertical border of the slope.

(4) According to the Rankine theory, the initial value of θ1in Fig.6 should be in the range of［270°， 45°+φ/2］(Liu and Cai,2017).In this figure,θirepresents the angle between the path segment and the slope crest.

(5) Generally, the TSS should be downward convex, which indicates that:

According to Sun et al. (2014), the parameters associated with MMACOA are set as follows: m = 20, v =1,β = 0.1,ρ = 0.3, Q = 2,τmax= 3 and τmin= 0.1.

The flow chart to search for the CSS using MMACOA is plotted in Fig. 7.

5. Numerical tests

Three typical numerical tests will be conducted using NMM.For comparison, following LEMs are adopted in this section:

(1) Bishop method (Bishop,1955).

(2) Janbu method (Janbu,1973).

(3) Spencer method (Spencer,1967).

(4) Morgenstern-Price method (Morgenstern and Price,1965).

Fig.15. Locations of potential critical slip surfaces for the slope with two layers of geo-materials obtained from different numerical methods:(a)Bishop,(b)Janbu,(c)Spencer,(d)Morgenstern-Price, and (e) Proposed method.

Table 4 Factors of safety obtained from different numerical methods for the second case.

Note that the CSS related to Janbu, Bishop, Spencer and Morgenstern-Price methods will be obtained with SLIDE,which is an interactive slope stability program.

5.1. A homogeneous slope

In the first case, we consider a homogeneous slope under the effect of gravity(Sun et al.,2014).Fig.8a presents a sketch for this slope. The NMM discretized model with 2220 manifold elements and 2351 PPs is plotted in Fig.8b.Table 1 lists the parameters used in this case.

Fig.16. (a) A sketch for the slope with three types of geo-materials and (b) the corresponding discretized model.

Table 5 Properties of material used in the third case.

Fig.17. Displacement contours (in m) of the slope with three types of geo-materials obtained from the numerical manifold method: (a) x- and (b) y-displacement.

Fig.18. Stress contours(in Pa)of the slope with three types of geo-materials obtained from the numerical manifold method: (a) σx, (b) σy, and (c) τxy.

The displacement and stress contours for the slope are presented in Figs. 9 and 10, respectively. The magnitudes of the ydisplacement and σycontours of the slope are all negative. These phenomena are mainly due to the fact that the slope is in a state of settlement caused by gravity.

Fig. 11 presents the potential locations of the CSSs of the homogeneous slope obtained from different numerical methods.The potential CSS predicted with MMACOA is almost in circular arc shape which agrees well with that using SLIDE. Table 2 lists the Fsvalues of the slope predicted with different numerical methods.The Fsvalue of the slope based on the proposed numerical method is 1.035 which also agrees well with the Bishop method(1.006)and the Morgenstern-Price method (1.005).

Fig.19. Locations of potential sliding surfaces based on different limit equilibrium methods for the slope with three types of geo-materials: (a) Bishop, (b) Janbu, (c) Spencer, (d)Morgenstern-Price, and (e) Proposed method.

Table 6 Factors of safety obtained from different numerical methods for the third case.

5.2. A slope with two layers of geo-materials

In the second case, we consider a slope with two layers of geomaterials. Fig.12a presents a sketch for this slope. The NMM discretized model with 917 manifold elements and 1011 PPs is plotted in Fig.12b.The parameters used in this case are listed in Table 3.A reference value of Fsrelated to this slope is 1.58(Zheng et al.,2006).

The displacement and stress contours of the slope are plotted in Figs. 13 and 14, respectively. Similar to the first case, the magnitudes of the y-displacement(Fig.13b)and σycontours(Fig.14b)of the slope are all negative.

Fig.15 presents the locations of CSSs of the slope with two layers of geo-materials obtained from different numerical methods. The potential CSS predicted with MMACOA agrees well with that from SLIDE. In addition, the Fsvalues of the slope computed using different numerical methods are tabulated in Table 4. The Fsvalue of the slope based on the proposed numerical method is 1.572 which is larger than those obtained from the Bishop method(1.41)and the Morgenstern-Price method (1.411), but agrees well with the reference solution from Zheng et al. (2006) (1.58).

5.3. A slope with three types of geo-materials

In the last case, we consider a slope with three types of geomaterials. Fig.16a presents a sketch for this slope. The NMM discretized model with 5112 manifold elements and 5529 PPs is plotted in Fig.16b. The parameters used in this case are listed in Table 5.

The displacement and stress contours of the slope are plotted in Figs.17 and 18, respectively. The stress field presented in Fig.18 is adopted to compute Fsof the slope using ELM.

Fig.19 presents the potential locations of the CSSs of the slope with three types of geo-materials obtained from different numerical methods. The potential CSS predicted with MMACOA agrees well with that by SLIDE. In addition, the Fsvalues of the slope computed using different numerical methods are listed in Table 6.The Fsvalue of the slope based on the proposed numerical method is 1.574, which agrees well with that from the Janbu method(Table 6) (1.562).

6. Conclusions

In this study, a numerical model to search for the CSS of the slope with the stress fields reproduced by NMM is proposed.In the numerical model, the MMACOA which is considered as one of the best performing algorithms for some optimization problems is adopted to search for the CSS. During the searching processes for the CSS, the stress fields of the slope (reproduced by an NMM analysis)and the ELM are combined to compute the factor of safety(Fs)of a series of TSSs.Procedures to compute the Fsvalues and the CSS are described.

With the developed numerical model, Fsand CSS for three different slope cases are investigated. The numerical results indicate that the Fsvalues and the CSSs of the three slopes can be accurately calculated with the proposed numerical model.

Note that the practical engineering problems are all in threedimensional (3D) space, and only very limited problems can be simplified into two-dimensional(2D)case.Therefore,development of 3D numerical model is essential.However,the numerical model proposed in the present work is only suitable for 2D slopes.Extension of the proposed numerical model for 3D slopes will be considered in our future work.

Declaration of competing interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

This study is supported by the Youth Innovation Promotion Association of Chinese Academy of Sciences (Grant No. 2020327)and the National Natural Science Foundation of China (Grant No.51609240).