Three-dimensional slope stability analysis using independent cover based numerical manifold and vector method

Engineering Geology 225 (2017) 83–95

Contents lists available at ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

Three-dimensional slope stability analysis using independent cover based numerical manifold and vector method Gaoyang Liu a, Xiaoying Zhuang a,⁎, Zhouquan Cui b a b

Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, China Yunnan Phosphate Chemical Group Co., LTD., Kunming 650600, China

a r t i c l e

i n f o

Article history: Received 1 September 2016 Received in revised form 13 February 2017 Accepted 22 February 2017 Available online 28 February 2017 Keywords: Independent cover Manifold method Slope stability Vector sum method Critical slip surface

a b s t r a c t A three-dimensional (3D) slope stability analysis method is presented in this paper based on three-dimensional independent cover based manifold method (ICMM3D) and vector sum method (VSM). ICMM3D is proposed in the framework of independent cover, which avoid the time-consuming and error-prone cover system generation of convention numerical manifold method (NMM). It is very suitable for the continuous/discontinuous deformation analysis in slope engineering. Then, with the stress field obtained by ICMM3D, VSM is employed to calculate the factor of safety, which is effective in computation and clear in physical meaning. In addition, a new strategy of discretization of the slip surface is proposed, which discretizes the slip surface into a set of calculating points in order to avoid errors caused by the special handling in the boundary columns in limit equilibrium method (LEM) and the tedious triangulation of the slip surface in VSM. Finally, genetic algorithm is employed to search for the critical slip surface. Numerical examples demonstrate the efficiency, accuracy and robustness of the proposed method. © 2017 Elsevier B.V. All rights reserved.

1. Introduction In the stability analysis of landslides, the rock masses are usually heterogeneous and contain large quantities of discontinuities, e.g., faults and joints. Besides recent advantages in remeshing procedure, see e.g. the seminal work of Areias et al. (2016a, 2013b, 2016b, 2015), Ghorashi et al. (2015), Jia et al. (2013), Nanthakumar et al. (2013), and Thai et al. (2014), the finite element method requires that the discontinuities are aligned to the discretization. A more elegant solution has been proposed in the context of so-called partition of unity enriched methods such as the extended finite element method (XFEM) (Belytschko and Black, 1999; Belytschko et al., 2001), generalized finite element method (GFEM) (Strouboulis et al., 2000) or specific improvements such as the smoothed extended finite element method (Nguyen-Xuan et al., 2008; Thai et al., 2012), isogeometric extended finite element methods (Ghasemi et al., 2015; Hughes et al., 2005; Nguyen-Thanh et al., 2014, 2015; Nguyen et al., 2016, 2015; Thai et al., 2015, 2016) or phantom node method (Chen et al., 2012; Song et al., 2006). In XFEM, the discontinuity can evolve arbitrarily in the discretization during crack propagation. However, modeling contact in the XFEM framework which is important for the analysis of landslides, remains a challenge. Meshfree methods (Belytschko et al., 1996; Li and Liu, 2002) offer an alternative to finite element methods. They are ⁎ Corresponding author. E-mail address: [email protected] (X. Zhuang).

http://dx.doi.org/10.1016/j.enggeo.2017.02.022 0013-7952/© 2017 Elsevier B.V. All rights reserved.

also based on continuum mechanics and can handle large deformations and contact (at least when Eulerian kernels are used) easier than finite element methods (Dilts, 1999; Libersky and Petschek, 1991; Rabczuk and Areias, 2006; Rabczuk et al., 2007a; Randles and Libersky, 1996; Ren et al., 2016; Vu-Bac et al., 2013; Zhuang et al., 2012b, 2014b) which is important in the simulation of landslides. Moreover, they are also well suited for modeling discrete fracture (Nguyen-Thanh et al., 2011; Nguyen et al., 2008; Rabczuk and Belytschko, 2007; Rabczuk et al., 2004, 2007b, 2010a, 2010b, 2010c; Zhuang et al., 2012a, 2014a, 2011). An alternative to methods which solve continuum mechanics equations are methods which account for the ‘fine-scale’ structure such as the discrete element method (DEM) (Cundall and Strack, 1979) or discontinuous deformation analysis (DDA) (Shi, 1988). They are capable of accounting for the faults and joints of the rock. However, capturing the macroscopic material behavior with those methods remains a major challenge. Therefore, coupling methods have been proposed that combine the strength of discrete methods and the finite element method (Munjiza et al., 1995); see also the contributions on multiscale methods for fracture (Amiri et al., 2014; Areias et al., 2013a; Budarapu et al., 2014; Quayum et al., 2015; Talebi et al., 2015; Vu-Bac et al., 2015). A very powerful method for the solution of continuumdiscontinuum problems is the numerical manifold method (NMM); see Ma et al. (2010) for a recent survey. The NMM was proposed by Shi (1991) and has been widely developed and extended to various two-dimensional problems, e.g. cover system generation (Cai and Wu,

84

G. Liu et al. / Engineering Geology 225 (2017) 83–95

2016; Chen and Li, 2015), fracture (Ma et al., 2009; Zheng and Xu, 2014), slope engineering (Koyama et al., 2012; Ning et al., 2011; Zheng et al., 2014), to name a few. The NMM has also been implemented in 3D (Cheng and Zhang, 2007, 2008; He and Ma, 2010; Jiang et al., 2010, 2009; Terada and Kurumatani, 2005) and applied to joint rock slope stability analysis (He et al., 2013). However, the efficient 3D implementation remains a challenge since it requires advanced and efficient contact algorithms. One of the most important issues is the cover system generation. For two-dimensional problems, Chen and Li (2015) considered the mathematical mesh as the union of mathematical elements (MEs) rather than mathematical covers (MCs) which made the cover system generation more efficient. Cai and Wu (2016) employed a pre-defined symbol function by which large amounts of computational geometries were replaced by computational algebras facilitating the generation process of the cover system. There are comparatively few papers in three dimensions on the generation of the cover system. Li and Zhang (2014) proposed a 3D manifold cutting method based on the 3D block cutting approach (Jing, 2000) which is able to generate arbitrary 3D MEs based on tetrahedral and hexahedral mesh covers. In order to solve the problems above, Cai and Liu (2015), Cai et al. (2013) proposed an independent cover based manifold method (ICMM). In the ICMM, various high-order cover functions are employed at the independent covers, and the corresponding elements are defined between the adjacent independent covers, which are different from the virtual springs in DDA and DEM. Complex algorithms for the cover system generation used in conventional NMM as well as the rank deficiency due to the linear dependence of the global degrees of freedom present in high-order NMM are completely avoided in ICMM. The continuous and discontinuous deformation analysis can be unified in one framework in the ICMM. Furthermore, the ICCM can be easily extended to three dimensions which is the key contribution of this manuscript. It will be extended to slope stability analysis for validation. The determination of the stability of slopes consists of two parts (Baker, 1980): calculate the factor of safety of a potential slip surface and search for the critical slip surface over all admissible slip surfaces. The definition of safety of factor is not unique. Zheng et al. (2006) suggested there are mainly two kinds of definitions in the abstract. One is the strength reserving definition, where the factor of safety is defined as “the factor by which the shear strength of the soil would have to be divided to bring the slope into a state of barely stability equilibrium” (Duncan, 1996), e.g. the limit equilibrium method (LEM) (Fellenius, 1939) and the strength reduction method (SRM) (Zeinkiewicz et al., 1975). By now, this definition is believed to be most familiar to engineers (Zheng et al., 2006). However, Ge (2010) pointed out that there are two problems for both the LEM and SRM. One is that the rationality of the strength reduction principle seems to be doubtful. For example, the cohesive force c and internal friction coefficient tanϕ of material divided by same F simultaneously is not very reasonable. If c and tanϕ are divided by different factors respectively, it will be complicated and infinite combining solutions could be obtained. The other is that because force is a vector, the algebraic sum of force is hold only for some special cases while the superposition principle of vectors should be hold in the analysis. The other definition of factor of safety is overloading definition. It defines the factor of safety as the scalar ratio of total resisting forces to total driving forces. The physical meaning is clear for a straight line or a circular slip line. However, it is questioned for a non-straight line or a non-circular slip line either, because it is neither the summation of force vectors in space nor the summation of projections of force vectors in a fixed direction. In order to overcome the problems above, Ge (2010) proposed the vector sum method (VSM), which estimates the stability of a slope by comparing the projection of the total resisting force with the total sliding force. The factor of safety is computed based on the real stress state and the vector sum algorithm, so the stress field needs to be calculated only once and the physical meaning is sound and clear. In this paper, the factor of safety is calculated using VSM based on the stress fields obtained by ICMM3D.

Fig. 1. An analysis domain with a joint.

The second step in slope stability analysis is finding the critical surface over all admissible slip surfaces. The aim is to find the best solution among available candidates by minimizing or maximizing an objective function. Optimization techniques are usually employed to determine the critical slip surface, e.g. pattern search methods (Bishop and Morgenstern, 1960; Greco and Gulla, 1985; Prater, 1979), calculus of variations (Baker and Garber, 1978; De Josselin De Jong, 1980; Friedli and Giger, 1978), dynamic programming methods (Baker, 1980; Yamagami et al., 1991), random search methods (Greco, 1996; Malkawi et al., 2001a, 2001b; Mowen, 2004; Yang et al., 2016), heuristic optimization methods (Bolton et al., 2003; Cheng, 2003; Cheng et al., 2008, 2005; Gao, 2015; Kahatadeniya et al., 2009; Kalatehjari et al., 2015), etc. In recent years, genetic algorithm has been widely used to locate the critical slip surface because of its elegance and efficiency. It has been found that genetic algorithm is a robust search technique which often gives global solution (Ahangar-Asr et al., 2012; Li et al., 2010; Manouchehrian et al., 2014; Sengupta and Upadhyay, 2009). In this paper, a search technique based on genetic algorithm and ellipsoidal shape is proposed to search for the critical slip surface. In addition, the potential slip surface during the process is discretized into a set of calculating points in order to avoid errors caused by the special handling in the boundary columns in limit equilibrium method (LEM) and the tedious triangulation of the slip surface in VSM. In this paper, ICMM3D is firstly proposed to calculate the stress field, then VSM is employed to obtain the factor of safety. After that the critical slip surface is determined based on genetic algorithm among the potential ellipsoidal shape slip surfaces. Finally, five examples are investigated to demonstrate the accuracy of ICMM3D and the presented slope stability analysis method.

2. The basic theory of ICMM3D 2.1. Stiffness matrix of the independent cover Fig. 1 illustrates an analysis domain with a joint. The associated finite element mesh based on hexahedra is shown in Fig. 2. The NMM provides a unified framework for both continuous and discontinuous problems (Chen and Li, 2015). In the conventional NMM, a mathematical mesh, which is the union of mathematical covers (MCs), is first assigned to model the problem. MCs have three properties: (1) MCs are arbitrarily defined by users; (2) they are independent of physical features, but their union must completely cover all physical features; and (3) they may overlap. The finite element mesh is used to define the mathematical mesh for the numerical manifold method. Considering any node, all elements having this node form a MC. In Fig. 2, there are 20 nodes and each node has a MC, e.g. MC of node 9 is the region 5-6-8-7-13-14-16-15 as shown in Fig. 3.

Fig. 2. Hexahedron element mesh.

G. Liu et al. / Engineering Geology 225 (2017) 83–95

85

Fig. 3. MC of node 9.

Fig. 5. Independent cover C1 and independent cover C2.

The second basic concept in the NMM is the physical cover (PC), which is the intersection of a MC with certain physical features. As a matter of fact, the PCs can also be understood as the subdivision of MCs by physical features. In Fig. 4 the MC of node 9 splits by the boundary and joint to two PCs: 91 and 92. The third basic concept is the manifold element (ME), which is defined as the common regions of physical covers. In the hexahedron mesh, usually the ME is the common part of 8 PCs (in this paper, we refer this type of ME as “complete overlapped covers”). In the NMM, the weight function is defined on the MC while the cover function is defined on the PC, which can be a constant, linear or higher order polynomial (Ma et al., 2009). With the weight function and cover function the global displacement function can be constructed on each ME based on the partition of unity property. By virtue of the dual cover system, NMM provides a robust numerical solution to both continuous and discontinuous problems in a unified way (He and Ma, 2010). However, the generation of the cover system is quite time-consuming and error-prone, which significantly restricts the development of NMM for complex geometries, especially in 3D. Linear dependence is another inevitable problem because of its partition of unity property (An et al., 2012). In order to solve these problems, an independent cover based manifold method (ICMM) has been proposed by Cai and Liu (2015), Cai et al. (2013). The ICMM will be expanded to 3D in the paper and subsequently named ICMM3D. The basic theory of ICMM3D is explained in this section. Therefore, consider Fig. 5 which illustrates two MEs C1: 9-10-1211-13-14-16-15 and C2: 13-14-16-15-17-18-20-19 from Fig. 2. In the conventional NMM, these two MEs should be the overlapping part of 8 PCs, respectively, which is the “complete overlapped cover”. In ICMM3D, we regard these MEs as “independent covers”. The displacement function in independent cover Ci is approximated as

Here (xc,yc, zc) is the centroid of the independent cover. From Eq. (1), the strain can be computed

uðx; y; zÞ ¼ Φðx; y; zÞ∙a

ð1Þ

where a is the degree of freedoms; Φðx; y; zÞ could be a constant, linear or higher order of polynomial 2

1 0 0 Φðx; y; zÞ ¼ 4 0 1 0 0 0 1

xc 0 0

0 xc 0

0 0 xc

yc 0 0

0 yc 0

Fig. 4. PC 91 and PC 92.

0 0 yc

zc 0 0

0 zc 0

0 0 zc

3 … … 5 ð2Þ …

ε ¼ Lu ¼ LΦa ¼ Ba

ð3Þ

where B is the matrix containing the derivatives of Φ and L is the differential operator matrix. The constitutive relation is given by σ ¼ DBa ¼ Sa

ð4Þ

where D is the elastic constant matrix. Hence the strain energy in independent cover Ci is ΠC ¼

1 T a ∫ BT ∙D∙BdΩa 2 Ω

ð5Þ

2.2. Stiffness matrix of the link element In ICMM3D, a thin-layer link element is applied between adjacent independent covers. Take independent covers C1 and C2 in Fig. 5 for example. Suppose that C1 and C2 are connected by a link element L12 with width b as shown in Fig. 6. From Eq. (1), the displacement functions of C1 and C2 can be expressed as uC 1 ¼ ΦC 1 ∙aC 1

ð6Þ

uC 2 ¼ ΦC 2 ∙aC 2

ð7Þ

where subscripts refer to the region. For example, uC1 and uC2 are the displacement vectors of independent cover C1 and C2 respectively. A local coordinate system ðx; y; zÞ is established in the adjacent plane of independent covers C1 and C2, in which the z-axis is perpendicular to the plane. The local coordinate is established in the adjacent plane, see Fig. 6. The width b is taken as b = l / 10000, where l is the mean value of the edge lengths of the link element L12. For a given point p in L12, there are two points p1 and p2 with the same local coordinate z adjacent to independent covers C1 and C2, respectively. From Eqs. (7) and (8) the local displacement of p1 and p2 is uC 1 ¼ λ∙uC 1

ð8Þ

uC 2 ¼ λ∙uC 2

ð9Þ

Fig. 6. Link element L12.

86

G. Liu et al. / Engineering Geology 225 (2017) 83–95

where 2

l1 λ ¼ 4 m1 n1

l2 m2 n2

3 l3 m3 5 n3

ð10Þ

is the coordinate transformation matrix, in which (l1,m1,n1), (l2,m2,n2), (l3,m3,n3) are the direction cosines of the local coordinate axes. Hence the local displacement of p can be computed by 8 9   < u = 1 z  1 z − uC 1 þ þ uC 2 uL ¼ v ¼ : ; 2 b 2 b w

ð11Þ

The strain in the link element is given by εL ¼ LuL ¼

    1  T 1 1 0 0 wC 2 −wC 1 0 vC 2 −vC 1 uC 2 −uC 1 b b b

ð12Þ

Since the εx , εy and γ xy approximately equal to 0, Eq. (13) can be rewritten in the form 8 9 8 9 < γzx = 1 < uC 2 −uC 1 = 1   v −vC 1 ¼ u −uC 1 εL ¼ γ yz ¼ : ; b : C2 ; b C2 wC 2 −w εz   C1 aC 2 1  ¼ BL ∙aL ¼ λ ΦC 2 −ΦC 1 aC 1 b

ð13Þ

 aC 2 contains the deaC 1 gree of freedoms. The stress tensor in the link element L12 is computed by −ΦC 1  and the vector aL ¼

8 9 < τzx = σL ¼ τyz ¼ DL ∙εL ¼ DL ∙BL ∙aL ¼ S L ∙aL : ; σz

ð14Þ

G DL ¼ 4 0 0

3 0 0 G 05 0 E

ð15Þ

According to Eq. (13) the εx, εy and γxy approximately equal to 0, so εL is rewritten to Eq. (14) for concision and convenience. Here G = E/ [2(1+ ν)] is shear modulus and E is elastic modulus. So the strain energy in link element L12 is ΠL ¼

1 T a ∫ BT ∙D∙BL dΩa 2 Ω L

1 1 Δus ¼ N s as ¼ Bs as b b

ð19Þ

The stress of joint element is obtained by σs ¼ Ds εs ¼ Ds Bs as

ð20Þ

where 2 kxy Ds ¼ 4 0 0

0 kyz 0

2 3 G 0 0 14 5 0 G ¼ 0 b 0 0 kzx

3 0 05 E

ð16Þ

ΠS ¼

1 T a ∫ N T ∙D∙N S dSa 2b S S

Assume the surface 17-18-20-19 of independent cover 13-14-1615-17-18-20-19 denoted as C3 is fixed as shown in Fig. 8. Similar to the link element in adjacent independent covers, a link element Lb with width b is assumed in the fixed surface. The displacement function in C3 is ð23Þ

The independent cover 5-6-7-8-9-10-11-12 is denoted as C3. C3 is divided by the joint into C31 and C32 as shown in Fig. 7. Goodman element is utilized. The displacement difference between the surfaces adjacent to C31 and C32 is i 1 h λ Φc32 −Φc31 as ¼ N s as b

ð17Þ

Here, Ns is the matrix containing the shape functions of the joint element and  aS ¼

ac32 ac31

 ð18Þ

ð22Þ

2.4. Displacement boundary

uc3 ¼ λ∙F e3 ∙ac3

2.3. Stiffness matrix of the joint element

Δus ¼

ð21Þ

kxy, kyz, kzx are the stiffness coefficients of the joint. The strain energy in joint element LS is

Here SL =DL ∙ BL and DL is the elasticity matrix given by 2

is the vector containing the degree of freedoms of the joint element. The strain of joint element is given by εs ¼



with BL ¼ 1b λ½ ΦC 2

Fig. 7. Joint element.

Fig. 8. Link element in displacement boundary.

G. Liu et al. / Engineering Geology 225 (2017) 83–95

F d ¼ ∫ N d T f 0 dS

so the displacement function in link element is 8 9 < u = 1 z  − ue3 ub ¼ v ¼ : ; 2 b w

S

ð24Þ

With Eq. (27), the strain energy in link element Lb can be expressed

87

ð34Þ

Assembling all the stiffness matrices and equivalent load vector, the equilibrium equation can be expressed as Ka ¼ F

ð35Þ

as 1 Π b ¼ aTe3 ∫ BTb ∙D∙Bb dΩae3 2 Ω

ð25Þ

with BTb ¼ 1b λΦe3 ; ae3 contains the degree of freedoms of C3. The strain energy of fixed point boundary is obtained in analogy. 2.5. Stress boundary  T Assume that a distributed load f 0 ¼ f x f y f z is applied on the surface 17-18-20-19 as shown in Fig. 9. As in the previous sections, a link element Ld is assumed on the surface. The strain energy in the link element Ld is then given by 1 Π d ¼ − aTe3 ∫ N Td f 0 dSae3 2 S

ð26Þ

where N d = λΦc3.

where K is the global stiffness matrix containing all sub-matrices Ki (i = c,L,S,b), a is the vector containing all DOFs of the independent covers and F is the equivalent load vector. When the rock mass contains quantities of discontinuities, the conventional NMM needs complex algorithm for cover system generation which involves large amounts of computational geometries. In ICMM3D the independent covers just need to be divided into several minor independent covers, which tremendously facilitates the implementation and improves the computational efficiency. Or in other words with only one physical cover in each independent cover, the implementation is much more natural and convenient. 3. Vector sum method The vector sum method (Ge, 2010) is proposed based on the vector property and the global potential sliding direction of the slope. The definition of factor of safety of the VSM is the ratio of the projection of the total mobilized anti-sliding force vector and the total sliding force vector on the potential sliding direction of the potential sliding mass.

2.6. Equilibrium equation 3.1. Derivation of safety factor with VSM With Eqs. (5), (17), (24), (28), (31), the total strain energy is Π ¼ ∑ðΠ C þ Π L þ Π S þ Π b þ Π d Þ

ð27Þ

Minimum of the total energy leads to

As shown in Fig. 10, σs ,στ ,σn are stress vector, shear stress and nor^ is the mal stress, respectively, at point A in the potential slip surface; n unit normal vector of the tangent plane at the same point on the poten^ is the unit vector of the potential sliding direction and tial slip surface, d S is the slip surface. Then

∂Π ¼0 ∂a

ð28Þ

^ σs ¼ σ∙n

ð36Þ

It can easily be shown that this leads finally to the following discrete system of equations:

^ Þn ^ σn ¼ ðσs ∙n

ð37Þ

∑ðK C aC þ K L aL þ K S aS þ K b ab −F d Þ ¼ 0

ð29Þ

στ ¼ σs −σn

ð38Þ

T

K C ¼ ∫ B ∙D∙BdΩ

ð30Þ

where σ is the stress tensor at the point A on the potential slip surface. The normal stress acting on the sliding mass at the point A can be expressed as

K L ¼ ∫ BTL ∙DL ∙BL dΩ

ð31Þ

σ0n ¼ −σn

1 ∫ N T ∙D∙N S dS bS S

ð32Þ

where Ω

Ω

KS ¼

ð39Þ

The factor of safety of the vector sum method can be expressed as

K b ¼ ∫ BTb ∙D∙Bb dΩ

ð33Þ

Ω

Fig. 9. Link element in stress boundary.

F¼

R T

ð40Þ

Fig. 10. The stress state at point A on the potential slip surface.

88

G. Liu et al. / Engineering Geology 225 (2017) 83–95

where R is the projection of total anti-sliding force on the opposite po^ T is the projection of total sliding force on potential sliding direction d;

^ They can be computed by tential sliding direction d.

 ^ dS R ¼ ∫ σ0s ∙ −d

ð41Þ

 ^ dS T ¼ ∫ σs ∙d

ð42Þ

S

s

In Eq. (41), limited anti-sliding stress vector σ′s can be expressed σ0s ¼ σ0τ þ σ0n

ð43Þ

Based on the Mohr-Coulomb criterion the limited shear stress can be expressed as ^r σ0τ ¼ ðc−σ n tanφÞd

ð44Þ

^ r is the unit direction of limited anti-sliding shear force of each where d section on the slip surface; c is cohesive force; φ is internal friction angle. ^ 3.2. Determination of whole potential sliding direction d According to the principle of friction, i.e. “the directions are always opposite to each other between the static friction force and relative sliding trend among two objects” (Ge, 2010), the opposite direction of the vector sum of all of limited anti-sliding forces on the slip surface is taken as the potential sliding direction. And the whole potential sliding direction can be expressed as ^ r dS − ∫ ðc−σ n tanφÞd ^¼ S d ^ r dS ∫ ðc−σ n tanφÞd

ð45Þ

S

^ r is the opposite projected direction of the In Eq. (45) the direction d ^ on the tangent plane of the point. The whole potential sliding direction d ^ is calculated using the direction d ^ r . So the implicit relationdirection d ship is set up between them. The direction of integral of shear stress along the slip surface is taken as the initial direction of potential sliding direction ∫ στ dS ^ ¼ S d 0 ∫ στ dS

ð46Þ

S

^ is the initial direction of potential sliding direction. With Eqs. where d 0 ^ can be (45) and (46), the final direction of potential sliding direction d

Fig. 11. Projection of the rotating ellipsoid on x-y plane.

Kalatehjari et al., 2015). In this paper a general ellipsoidal shape (Huang and Tsai, 2000) is adopted which can be converted to spherical or cylindrical shapes during the analysis. The ellipsoid can be expressed as  2 x cosθxy þ y sinθxy −X c cosθxy −Y c y sinθxy 2 a  2 −x sinθxy þ y cosθxy þ X c sinθxy −Y c cosθxy ðz−Z c Þ2 þ þ 2 c2 b ¼1

where Xc, Yc, Zc are the coordinates of the center of the ellipsoid in the x-, y-, and z-directions, respectively; a, b, and c are semi-radii of the ellipsoid in the x-, y-, and z-directions, respectively; θxy is the rotation angle of the ellipsoid in the x-y plane as shown in Fig. 11. It should be noted that this ellipsoid can be easily transformed into a spherical or cylindrical slip surface by using three and two equal semi-radii, respectively. 4.2. Discretization of the slip surface The equilibrium equation of LEM is set up based on columns, so the sliding mass needs to be discretized into a number of columns by two series of crossed parallel lines (Ahangar-Asr et al., 2012; Huang and Tsai, 2000; Yamagami and Jiang, 1997). The whole slope is discretized into columns first, and then the active columns defined as the columns inside the boundary of the sliding mass is selected by the locations of the centroids in the x-y plane. The columns should be small enough to guarantee the accuracy of the selection in the boundary. In LEM, the normal direction is determined by the base plane of the column while the slip surface is actually an ellipsoidal shape, which affects the accuracy of the analysis. Ge (2010) discretized the slip surface into triangular patches using VSM in the 3D slope stability analysis. In this paper, several calculating points are utilized for convenience and accuracy. As shown in Fig. 12, a set of points are arranged inside the projection of the slip surface in

obtained by iteration ^ d

^

i−1 −di bε

ð47Þ

where ε is the stopping tolerance. With the direction of potential sliding

^ the factor of safety of 3D slopes can be easily obtained by Eq. direction d, (40). 4. Searching for the critical slip surface 4.1. Shape of the slip surface

The shape of the slip surface is quite complicated and simplifications are needed in the slope stability analysis. Different surfaces such as spherical, ellipsoidal, and irregular shape have been used for the failure surface in 3D analysis (Alkasawneh et al., 2008; Huang and Tsai, 2000;

ð48Þ

Fig. 12. Calculating points on the slip surface.

G. Liu et al. / Engineering Geology 225 (2017) 83–95

x-y plane along the x- and y-directions, and the corresponding points on the slip surface is the calculating points. For each calculating point (x, y, ^ (nxCSS, nyCSS, nzCSS) can be computed by z), the normal direction n  nxCSS ¼

 x cosθxy þ y sinθxy −X c cosθxy −Y c y sinθxy cosθxy

þ  nyCSS ¼

R2y

ð49Þ

Table 1 Material properties for the five examples used in this study. Example no.

Material no.

γ (kN/m3)

Cohesion (kPa)

φ (°)

ν

E (kN/m2)

1 2 3

1 1 1 2 1 1

0 19.2 19.2 19.2 17 18

0 29.3 29.3 0 15 10

0 20 20 0 20 10

0.3 0.3 − − − −

3E+4 1E+6 – − − −

4 5

 x cosθxy þ y sinθxy −X c cosθxy −Y c y sinθxy sinθxy 2

þ

nzCSS ¼

R2x    −x sinθxy þ y cosθxy þ X c sinθxy −Y c cosθxy − sinθxy

89

Rx   −x sinθxy þ y cosθxy þ X c sinθxy −Y c cosθxy cosθxy R2y

z−Z c R2c

ð50Þ

ð51Þ

The stress state σ can be obtained by ICMM3D. With the stress state and normal direction of each calculating points. The factor of safety of the slip surface can be acquired by VSM easily. This method avoids the selection of columns in the boundary and the trial triangulation of the slip surface.

analogy for each parameter would greatly reduce the computational cost. In this paper, the recommended height of the searching region is twice the slope height. The parameters ranges are estimated mainly based on this principle. (2) Evaluate the fitness of each individual in the population. Each binary individual is decoded into 7 decimal parameters. Then slip surface is obtained from Eq. (48) which is needed to compute the corresponding factor of safety Fi. As the GA implementation here allows only for optimization to find a maximum value, the fitness function is transformed to f fit ¼

1

ð53Þ

F 2i

4.3. Searching for the critical slip surface based on genetic algorithm Genetic algorithms (GAs) are stochastic global search and optimization methods that mimic the metaphor of natural biological evolution (Holland, 1975). GAs operate on a population of potential solutions applying the principle of survival of the fittest to produce successively better approximations to a solution. Via the operations of selection, crossover, and mutation the GA will converge over successive generations towards the global optimum. From Eq. (48) the ellipsoidal slip surface is determined by 7 parameters including the centroid coordinates (Xc,Yc,Zc), 3 semi-radiuses a, b, c and rotation angle θxy. The objective function can be expressed in terms of these 7 parameters   F i ¼ f X c ; Y c ; Z c ; a; b; c; θxy

ð52Þ

With the objective function the procedure of search for the critical slip surface based on GA is (1) Generate a population of individuals. Each individual is a combination of 7 pieces of binary code encoded by the 7 parameters in Eq. (48) in given ranges, which means that each individual can be decoded to an ellipsoidal slip surface. It should be noted that suitable value ranges by engineering experience and engineering

Here the square of Fi was used to assist the GA in discriminating the optimal value. (3) Create a new population by selection operator, crossover operator and mutation operator. (4) Repeat step (2)–(3) until fixed number of generations is reached. Afterwards, the critical slip and the corresponding factor of safety are obtained. The procedure of the method for 3D slope stability analysis is illustrated in Fig. 13. 5. Examples Five validation examples are presented in this section. Example 1 is used to verify the accuracy of ICMM3D by performing a cantilever bending test. In example 2 and example 3 the performance of the proposed method is verified by calculating the factor of safety of specified slip surfaces. In example 4 and example 5 the cylindrical and ellipsoidal slip surfaces are searched, respectively, in order to verify the performance of the proposed method in searching for the critical slip surface of different shape. The cover function is second order polynomial in these five examples. The material properties are shown in Table 1. 5.1. Cantilever beam As shown in Fig. 14 a cantilever is analyzed to verify the accuracy of ICMM3D. This problem can actually be modeled in two dimensions but we have used a full 3D since an analytical solution is available. The left end is clamped while point A at the center of the free end is subjected by a concentrated force P = 10 kN. The dimensions of the cantilever

Fig. 13. Procedure of the proposed method for 3D slope stability analysis.

Fig. 14. Geometry of the cantilever.

90

G. Liu et al. / Engineering Geology 225 (2017) 83–95 Table 2 Deflection at point A. Nodes

Analytical solution (mm)

FEM (mm)

ICMM3D (mm)

99 525

17.125 17.125

16.672 16.993

16.841 17.141

Fig. 15. Mesh with 99 nodes for ICMM3D.

are L = 50 m, D = 10 m and B = 10 m. The analytical solution (Han and Atluri, 2004) is " !# 9 > Py D2 > > ð6L−3xÞx þ ð2 þ ν Þ y2 − ux ¼ − > > 6EI 4 > > " #> > = 2 P D x 3νy2 ðL−xÞ þ ð4 þ 5νÞ þ ð3L−xÞx2 uy ¼ > 6EI 4 > ! > > > > > P ðL−xÞy P D2 > > ; σ y ¼ 0; τxy ¼ − −y2 σx ¼ ; I 2I 4

I¼

D3 B 12

ð54Þ

ð55Þ

The cantilever is model with 99 and 525 nodes. Fig. 15 shows the mesh with 99 nodes. The deflection along the neutral axis is shown in Fig. 16 and the deflection at point A is listed in Table 2. It is found that ICMM3D is more accurate than FEM using the same mesh and displacement is getting closer to the analytical solution with finer mesh. Then the convergence test is carried out to further validate ICMM3D. To this end, four different meshes with node numbers as 99, 525, 2562, 6629 are investigated in which cover functions with first, second and third order polynomial are employed respectively. The relative energy norm errors away from the analytical solution are plotted against the node number in Fig. 16. As shown in Fig. 16, with the increase of node number and cover function order, the relative energy norm error decreases rapidly (Fig. 17).

Fig. 17. Convergence results of ICMM3D in example 1.

5.2. Slope of homogeneous material In the second example (Zhang, 1988), the slope is homogeneous and the slip surface is part of ellipsoid. The parameters of the ellipsoidal slip surface are Xc = 100, Yc = −36.6, Zc = 27.4, a = 66, b = 24.4, c = 24.4. The central section of the slip body is illustrated in Fig. 18. Note that the weak layer does not appear in this example and we only refer to the geometry used in Fig. 18. The material of the slope is homogeneous and only the geometry is discretized with 4623, 5822, 7398 nodes, respectively. The mesh with 5822 nodes is shown in Fig. 18. The critical slip surface and calculating points are shown in Fig. 20. The factor of safety obtained by the method proposed in this paper is compared with other results in Table 3. The agreement is excellent.

Fig. 18. Central section of slip surface.

Note the presented method requires less nodes than FEM based VSM (Ge, 2010) to obtain a reasonable factor of safety. The factor of safety calculated with 7398 nodes is 2.9% smaller than Ge (2010).

5.3. Slope with weak interlayer The slope with a weak interlayer, where an ellipsoidal slip surface is cut off by a weak layer (Fig. 19), was introduced by Zhang (1988). The slip surface is a combination of ellipsoid and plane. The critical slip surface and calculating points are shown in Fig. 21. As can be seen in Fig. 21, the slip surface is cut by the weak layer (the dark flatten off plane in the bottom of the slip surface), and the calculating points are distributed on the ellipsoid and plane respectively. The result is shown in Table 4. The factor of safety obtained by the presented method with 7398 nodes is 1.530, which is insignificant (0.2%) less compared to the one obtained by Ge (2010).

Table 3 Comparisons of factor of safety in example 2.

Fig. 16. Deflection along the neutral axis.

Zhang (1988)

Ge (2010)

The proposed method

2.122

7094 nodes 2.042

4623 nodes 2.142

5822 nodes 2.092

7398 nodes 2.061

G. Liu et al. / Engineering Geology 225 (2017) 83–95

91

Fig. 19. Mesh used for cover generation for slope with homogeneous material.

Fig. 20. Critical slip surface and calculating points for slope with homogeneous material.

Fig. 21. Critical slip surface and calculating points for slope with weak interlayer.

5.4. Alkasawneh slope Fig. 22 illustrates the geometry of the slope which has been studied for instance by Alkasawneh et al. (2008) and Kalatehjari et al. (2015). A 3D model was developed based on this 2D section

in which the third dimension was extended by 100 m. Fig. 23 shows the mesh of the slope with 2058 nodes. In this example, a cylindrical slip surface is employed to determine the critical slip surface by set c to infinity or just eliminate the monomial containing c.

92

G. Liu et al. / Engineering Geology 225 (2017) 83–95

Table 4 Comparisons of factor of safety in example 3. Zhang (1988)

Ge (2010)

The proposed method

1.533

7094 nodes 1.545

4623 nodes 1.562

5822 nodes 1.546

7398 nodes 1.530

Fig. 25. Evolution of the factor of safety for Alkasawneh slope.

Fig. 22. Geometry of slope used in Alkasawneh et al. (2008).

Fig. 26. Critical slip surface and calculating points for Alkasawneh slope.

Fig. 23. Mesh for ICMM3D for Alkasawneh slope.

Fig. 24 shows the evolution of the average and best fitness of the population. With the selection operator, crossover operator and mutation operator of GA, the best fitness is improved in the early generations and come to the maximum at last during the evolution. The average fitness of the population continues changing during the evolution mainly due to the crossover operator and the mutation operator, which avoids “Premature Convergence” leading to the local optima. Fig. 25 illustrates the minimum factor of safety during the evolution. The critical slip surface and the calculating points are shown in Fig. 26. As shown in Table 5, the factor of safety obtained by our method is similar to others, 0.6% less than (Alkasawneh et al., 2008) and 0.7% larger than (Kalatehjari et al., 2015).

5.5. Slope subjected to surface load The last example is performed to verify the ability of the presented method to determine the critical slip surface with a general ellipsoidal shape and the minimum factor of safety in comparison to previous studies. This example was initially analyzed by Yamagami et al. (1991) and re-analyzed by Yamagami and Jiang (1997) and Kalatehjari et al. (2015). The example involves a homogeneous slope with a gradient of 2:1 subjected to a square load of 50 kPa on the top. The uniform load is applied on a square surface with 8 m sides at the top center of the slope. Fig. 27 illustrates the geometry of example 5. Fig. 28 shows the mesh with 4880 nodes. The evolution of fitness and the minimum factor of safety are shown in Figs. 29 and 30, respectively. The revolution speed is relatively slow compared with example of Alkasawneh slope employing cylindrical slip surfaces. The main reason is that by employing the ellipsoidal shape, large quantities of disqualified slip surfaces is generated during the evolution, which dramatically decreased the average fitness value. The critical slip surface and calculating points is shown in Fig. 31. The results illustrated in Table 6 shows that the factor of safety obtained by the presented method in this paper is in good agreement with others.

Table 5 Comparisons of factor of safety for Alkasawneh slope.

Fig. 24. Evolution of the fitness of for Alkasawneh slope.

Alkasawneh et al. (2008)

Kalatehjari et al. (2015)

The proposed method

1.80

1.78

1.793

G. Liu et al. / Engineering Geology 225 (2017) 83–95

93

Fig. 30. Evolution of the factor of safety in example 5.

Fig. 27. The plane view and central cross section of the slope subject to surface load.

Fig. 31. Critical slip surface and calculating points in example 5.

Fig. 28. Mesh for ICMM3D in example 5.

of convention numerical manifold method (NMM). It is very suitable for the continuous/discontinuous deformation analysis in slope engineering. 2) VSM is based on the stress field in the natural state and factor sum algorithm. It does not need the assumptions of the inter-slice force in LEM and the trial iterative computations in SRM. It is effective in computation and has a clear physical meaning. 3) The slip surface is discretized into a set of calculating points in order to eliminate errors caused by the special handling in the boundary columns in the limit equilibrium method (LEM) and the tedious triangulation of the slip surface in VSM. 4) Genetic algorithm is employed to search for the critical slip surface. The computer code is able to search for the critical slip surface of different shape, including sphere, cylinder and ellipsoid. The cover functions of independent covers in ICMM3D are polynomial in this paper. It should be noted that the cover functions in ICMM3D could be other kinds of functions in crack propagation analysis. And the critical slip surface could be determined by the crack propagation and the corresponding factor of safety could be obtained by VSM. That will be the topic for the future work.

Fig. 29. Evolution of the fitness in example 5.

6. Conclusions A 3D slope stability analysis method using ICMM3D and VSM is proposed in this paper. Its advantages are: 1) ICMM3D is proposed in the framework of independent cover, which avoids the time-consuming and error-prone cover system generation

Table 6 Comparisons of factor of safety in example 5. Yamagami et al. (1991)

Yamagami et al. (1991)

Kalatehjari et al. (2015)

The proposed method

1.14

1.03

0.95

0.972

94

G. Liu et al. / Engineering Geology 225 (2017) 83–95

References Ahangar-Asr, A., Toufigh, M.M., Salajegheh, A., 2012. Determination of the most probable slip surface in 3D slopes considering the effect of earthquake force direction. Comput. Geosci. 45:119–130. http://dx.doi.org/10.1016/j.cageo.2011.10.024. Alkasawneh, W., Malkawi, A.I.H., Nusairat, J.H., Albataineh, N., 2008. A comparative study of various commercially available programs in slope stability analysis. Comput. Geotech. 35, 428–435. Amiri, F., Millan, D., Shen, Y., Rabczuk, T., Arroyo, M., 2014. Phase-field modeling of fracture in linear thin shells. Theor. Appl. Fract. Mech. 69:102–109. http://dx.doi.org/10. 1016/j.tafmec.2013.12.002. An, X.M., Zhao, Z.Y., Zhang, H.H., Li, L.X., 2012. Investigation of linear dependence problem of three-dimensional partition of unity-based finite element methods. Comput. Methods Appl. Mech. Eng. 233–236:137–151. http://dx.doi.org/10.1016/j.cma.2012. 04.010. Areias, P., Dias-da-Costa, D., Sargado, J.M., Rabczuk, T., 2013a. Element-wise algorithm for modeling ductile fracture with the Rousselier yield function. Comput. Mech. 52: 1429–1443. http://dx.doi.org/10.1007/s00466-013-0885-0. Areias, P., Rabczuk, T., Dias-da-Costa, D., 2013b. Element-wise fracture algorithm based on rotation of edges. Eng. Fract. Mech. 110:113–137. http://dx.doi.org/10.1016/j. engfracmech.2013.06.006. Areias, P., Reinoso, J., Camanho, P., Rabczuk, T., 2015. A constitutive-based element-by-element crack propagation algorithm with local mesh refinement. Comput. Mech. 56, 291–315. Areias, P., Msekh, M.A., Rabczuk, T., 2016a. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Eng. Fract. Mech. 158:116–143. http://dx.doi.org/10.1016/j.engfracmech.2015.10.042. Areias, P., Rabczuk, T., Msekh, M., 2016b. Phase-field analysis of finite-strain plates and shells including element subdivision. Comput. Methods Appl. Mech. Eng. Baker, R., 1980. Determination of the critical slip surface in slope stability computations. Int. J. Numer. Anal. Methods Geomech. 4, 333–359. Baker, R., Garber, M., 1978. Theoretical analysis of the stability of slopes. Geotechnique 28, 395–411. Belytschko, T., Black, T., 1999. Elastic crack growth in finite elements with minimal remeshing. Int. J. Numer. Methods Eng. 45:601–620. http://dx.doi.org/10.1002/ (sici)1097-0207(19990620)45:5b601::aid-nme598N3.0.co;2-s. Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., Krysl, P., 1996. Meshless methods: an overview and recent developments. Comput. Methods Appl. Mech. Eng. 139:3–47. http://dx.doi.org/10.1016/s0045-7825(96)01078-x. Belytschko, T., Moës, N., Usui, S., Parimi, C., 2001. Arbitrary discontinuities in finite elements. Int. J. Numer. Methods Eng. 50, 993–1013. Bishop, A., Morgenstern, N., 1960. Stability coefficients for earth slopes. Geotechnique 10, 129–153. Bolton, H.P.J., Heymann, G., Groenwold, A.A., 2003. Global search for critical failure surface in slope stability analysis. Eng. Optim. 35:51–65. http://dx.doi.org/10.1080/ 0305215031000064749. Budarapu, P.R., Gracie, R., Bordas, S.P.A., Rabczuk, T., 2014. An adaptive multiscale method for quasi-static crack growth. Comput. Mech. 53:1129–1148. http://dx.doi.org/10. 1007/s00466-013-0952-6. Cai, Y.C., Liu, G.Y., 2015. High-order manifold method with independent covers. J. Tongji Univ. 1794–1800. Cai, Y.C., Wu, J., 2016. A robust algorithm for the generation of integration cells in Numerical Manifold Method. Int. J. Impact Eng. 90:165–176. http://dx.doi.org/10.1016/j. ijimpeng.2015.10.015. Cai, Y.C., Zhu, H.H., Zhuang, X.Y., 2013. A continuous/discontinuous deformation analysis (CDDA) method based on deformable blocks for fracture modeling. Front. Struct. Civ. Eng. 7, 369–378. Chen, Y.L., Li, L.X., 2015. An improved numerical manifold method and its application. Eng. Anal. Bound. Elem. 52:120–129. http://dx.doi.org/10.1016/j.enganabound.2014.11.033. Chen, L., Rabczuk, T., Bordas, S.P.A., Liu, G.R., Zeng, K.Y., Kerfriden, P., 2012. Extended finite element method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth. Comput. Methods Appl. Mech. Eng. 209:250–265. http://dx.doi.org/ 10.1016/j.cma.2011.08.013. Cheng, Y.M., 2003. Location of critical failure surface and some further studies on slope stability analysis. Comput. Geotech. 30:255–267. http://dx.doi.org/10.1016/S0266352X(03)00012-0. Cheng, M.Y., Zhang, H.Y., 2007. Formulation of a three-dimensional numerical manifold method with tetrahedron and hexahedron elements. Rock Mech. Rock. Eng. 41: 601–628. http://dx.doi.org/10.1007/s00603-006-0120-9. Cheng, M.Y., Zhang, H.Y., 2008. Three-dimensional numerical manifold method - tetrahedron and hexahedron mesh. Int. Rock Mech. Rock Eng. 41 (4), 601–628. Cheng, Y.M., Liu, H.T., Wei, W.B., Au, S.K., 2005. Location of critical three-dimensional nonspherical failure surface by NURBS functions and ellipsoid with applications to highway slopes. Comput. Geotech. 32:387–399. http://dx.doi.org/10.1016/j.compgeo. 2005.07.004. Cheng, Y.M., Li, L., Lansivaara, T., Chi, S.C., Sun, Y.J., 2008. An improved harmony search minimization algorithm using different slip surface generation methods for slope stability analysis. Eng. Optim. 40:95–115. http://dx.doi.org/10.1080/03052150701618153. Cundall, P.A., Strack, O.D., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Dilts, G.A., 1999. Moving-least-squares-particle hydrodynamics—I. Consistency and stability. Int. J. Numer. Methods Eng. 44, 1115–1155. Duncan, J., 1996. State of the art: limit equilibrium and finite-element analysis of slopes. J. Geotech. Eng. 122:577–596. http://dx.doi.org/10.1061/(ASCE)0733-9410(1996)122: 7(577).

Fellenius, W., 1939. Erdstatisch Berechnungen, Berlin W. Ernst und Sohn revised edition. Friedli, P., Giger, M., 1978. Calculus of Variations Applied to Stability of Slopes. Thomas Telford Services Ltd, Thomas Telford House, 1 Heron Quay, London, England E14 4JD, pp. 222–223. Gao, W., 2015. Determination of the noncircular critical slip surface in slope stability analysis by meeting ant colony optimization. J. Comput. Civ. Eng. 30:06015001. http://dx. doi.org/10.1061/(ASCE)CP.1943-5487.0000475. Ge, X.-R., 2010. The vector sum method: a new approach to calculating the safety factor of stability against sliding for slope engineering and dam foundation problems. In: Chen, Y., Zhan, L., Tang, X. (Eds.), Advances in Environmental Geotechnics: Proceedings of the International Symposium on Geoenvironmental Engineering in Hangzhou, China, September 8–10, 2009. Springer, Berlin, Heidelberg, pp. 99–110. Ghasemi, H., Brighenti, R., Zhuang, X., Muthu, J., Rabczuk, T., 2015. Optimal fiber content and distribution in fiber-reinforced solids using a reliability and NURBS based sequential optimization approach. Struct. Multidiscip. Optim. 51:99–112. http://dx. doi.org/10.1007/s00158-014-1114-y. Ghorashi, S.S., Valizadeh, N., Mohammadi, S., Rabczuk, T., 2015. T-spline based XIGA for fracture analysis of orthotropic media. Comput. Struct. 147:138–146. http://dx.doi. org/10.1016/j.compstruc.2014.09.017. Greco, V.R., 1996. Efficient Monte Carlo technique for locating critical slip surface. J. Geotech. Eng. 122, 517–525. Greco, V., Gulla, G., 1985. Slip surface search in slope stability analysis. Riv. Ital. Geotecnica 19 (4), 189–198. Han, Z.D., Atluri, S.N., 2004. Meshless local Petrov-Galerkin (MLPG) approaches for solving 3D problems in elasto-statics. Comput. Model. Eng. Sci. 6, 169–188. He, L.E.I., Ma, G., 2010. Development of 3d numerical manifold method. Int. J. Comput. Methods 07:107–129. http://dx.doi.org/10.1142/S0219876210002088. He, L., An, X.M., Ma, G.W., Zhao, Z.Y., 2013. Development of three-dimensional numerical manifold method for jointed rock slope stability analysis. Int. J. Rock Mech. Min. Sci. 64:22–35. http://dx.doi.org/10.1016/j.ijrmms.2013.08.015. Holland, J.H., 1975. Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence. U Michigan Press. Huang, C.C., Tsai, C.C., 2000. New method for 3D and asymmetrical slope stability analysis. J. Geotech. Geoenviron. 126:917–927. http://dx.doi.org/10.1061/(asce)10900241(2000)126:10(917). Hughes, T.J.R., Cottrell, J.A., Bazilevs, Y., 2005. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl. Mech. Eng. 194:4135–4195. http://dx.doi.org/10.1016/j.cma.2004.10.008. Jia, Y., Zhang, Y.J., Xu, G., Zhuang, X.Y., Rabczuk, T., 2013. Reproducing kernel triangular Bspline-based FEM for solving PDEs. Comput. Methods Appl. Mech. Eng. 267:342–358. http://dx.doi.org/10.1016/j.cma.2013.08.019. Jiang, Q.H., Zhou, C.B., Li, D.Q., 2009. A three-dimensional numerical manifold method based on tetrahedral meshes. Comput. Struct. 87:880–889. http://dx.doi.org/10. 1016/j.compstruc.2009.03.002. Jiang, Q.-h., Deng, S.-s., Zhou, C.-b., Lu, W.-b., 2010. Modeling unconfined seepage flow using three-dimensional numerical manifold method. J. Hydrodyn. Ser. B 22: 554–561. http://dx.doi.org/10.1016/S1001-6058(09)60088-3. Jing, L.R., 2000. Block system construction for three-dimensional discrete element models of fractured rocks. Int. J. Rock Mech. Min. Sci. 37:645–659. http://dx.doi.org/10.1016/ s1365-1609(00)00006-x. De Josselin De Jong, G., 1980. Application of the calculus of variations to the vertical cut off cohesive frictionless soil. Geotechnique 30, 1–16. Kahatadeniya, K.S., Nanakorn, P., Neaupane, K.M., 2009. Determination of the critical failure surface for slope stability analysis using ant colony optimization. Eng. Geol. 108: 133–141. http://dx.doi.org/10.1016/j.enggeo.2009.06.010. Kalatehjari, R., Arefnia, A., Rashid, A.S.A., Ali, N., Hajihassani, M., 2015. Determination of three-dimensional shape of failure in soil slopes. Can. Geotech. J. 52, 1283–1301. Koyama, T., Ohnishi, Y., Miki, S., Nakai, T., Maruki, Y., Yagi, K., Kobayashi, T., 2012. Application of Manifold Method (MM) to the stability problems for cut slopes along the national roads. Geomech. Geoeng. 7, 39–56. Li, S., Liu, W.K., 2002. Meshfree and particle methods and their applications. Appl. Mech. Rev. 55:1–34. http://dx.doi.org/10.1115/1.1431547. Li, H.F., Zhang, G.X., 2014. Researches on the generation of three-dimensional manifold element under FEM mesh cover. Math. Probl. Eng. 17. http://dx.doi.org/10.1155/2014/ 140180. Li, Y.C., Chen, Y.M., Zhan, T.L., Ling, D.S., Cleall, P.J., 2010. An efficient approach for locating the critical slip surface in slope stability analyses using a real-coded genetic algorithm. Can. Geotech. J. 47:806–820. http://dx.doi.org/10.1139/t09-124. Libersky, L.D., Petschek, A., 1991. Smooth particle hydrodynamics with strength of materials. Advances in the Free-Lagrange Method Including Contributions on Adaptive Gridding and the Smooth Particle Hydrodynamics Method. Springer, pp. 248–257. Ma, G.W., An, X.M., Zhang, H.H., Li, L.X., 2009. Modeling complex crack problems using the numerical manifold method. Int. J. Fract. 156:21–35. http://dx.doi.org/10.1007/ s10704-009-9342-7. Ma, G., An, X., He, L.E.I., 2010. The numerical manifold method: a review. Int. J. Comput. Methods 07:1–32. http://dx.doi.org/10.1142/S0219876210002040. Malkawi, A.I.H., Hassan, W.F., Sarma, S.K., 2001a. An efficient search method for finding the critical circular slip surface using the Monte Carlo technique. Can. Geotech. J. 38:1081–1089. http://dx.doi.org/10.1139/cgj-38-5-1081. Malkawi, A.I.H., Hassan, W.F., Sarma, S.K., 2001b. Global search method for locating general slip surface using Monte Carlo techniques. J. Geotech. Geoenviron. 127:688–698. http://dx.doi.org/10.1061/(asce)1090-0241(2001)127:8(688). Manouchehrian, A., Gholamnejad, J., Sharifzadeh, M., 2014. Development of a model for analysis of slope stability for circular mode failure using genetic algorithm. Environ. Earth Sci. 71:1267–1277. http://dx.doi.org/10.1007/s12665-013-2531-8.

G. Liu et al. / Engineering Geology 225 (2017) 83–95 Mowen, X., 2004. A simple Monte Carlo method for locating the three-dimensional critical slip surface of a slope. Acta Geol. Sin. Engl. Ed. 78, 1258–1266. Munjiza, A., Owen, D.R.J., Bicanic, N., 1995. A combined finite-discrete element method in transient dynamics of fracturing solids. Eng. Comput. 12:145–174. http://dx.doi.org/ 10.1108/02644409510799532. Nanthakumar, S.S., Lahmer, T., Rabczuk, T., 2013. Detection of flaws in piezoelectric structures using extended FEM. Int. J. Numer. Methods Eng. 96:373–389. http://dx.doi.org/ 10.1002/nme.4565. Nguyen, V.P., Rabczuk, T., Bordas, S., Duflot, M., 2008. Meshless methods: a review and computer implementation aspects. Math. Comput. Simul. 79:763–813. http://dx.doi. org/10.1016/j.matcom.2008.01.003. Nguyen, V.P., Anitescu, C., Bordas, S.P.A., Rabczuk, T., 2015. Isogeometric analysis: an overview and computer implementation aspects. Math. Comput. Simul. 117:89–116. http://dx.doi.org/10.1016/j.matcom.2015.05.008. Nguyen, B.H., Tran, H.D., Anitescu, C., Zhuang, X., Rabczuk, T., 2016. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems. Comput. Methods Appl. Mech. Eng. 306, 252–275. Nguyen-Thanh, N., Nguyen-Xuan, H., Bordas, S.P.A., Rabczuk, T., 2011. Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids. Comput. Methods Appl. Mech. Eng. 200:1892–1908. http://dx.doi.org/10. 1016/j.cma.2011.01.018. Nguyen-Thanh, N., Muthu, J., Zhuang, X., Rabczuk, T., 2014. An adaptive three-dimensional RHT-splines formulation in linear elasto-statics and elasto-dynamics. Comput. Mech. 53, 369–385. Nguyen-Thanh, N., Valizadeh, N., Nguyen, M.N., Nguyen-Xuan, H., Zhuang, X., Areias, P., Zi, G., Bazilevs, Y., De Lorenzis, L., Rabczuk, T., 2015. An extended isogeometric thin shell analysis based on Kirchhoff-Love theory. Comput. Methods Appl. Mech. Eng. 284: 265–291. http://dx.doi.org/10.1016/j.cma.2014.08.025. Nguyen-Xuan, H., Rabczuk, T., Bordas, S., Debongnie, J.F., 2008. A smoothed finite element method for plate analysis. Comput. Methods Appl. Mech. Eng. 197:1184–1203. http:// dx.doi.org/10.1016/j.cma.2007.10.008. Ning, Y.J., An, X.M., Ma, G.W., 2011. Footwall slope stability analysis with the numerical manifold method. Int. J. Rock Mech. Min. Sci. 48:964–975. http://dx.doi.org/10. 1016/j.ijrmms.2011.06.011. Prater, E.G., 1979. Yield acceleration for seismic stability of slopes. J. Geotech. Eng. Div. 105, 682–687. Quayum, M.S., Zhuang, X., Rabczuk, T., 2015. Computational model generation and RVE design of self-healing concrete. Front. Struct. Civ. Eng. 9, 383–396. Rabczuk, T., Areias, P., 2006. A Meshfree Thin Shell for Arbitrary Evolving Cracks Based on an Extrinsic Basis. Rabczuk, T., Belytschko, T., 2007. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput. Methods Appl. Mech. Eng. 196:2777–2799. http://dx.doi.org/10.1016/j.cma.2006.06.020. Rabczuk, T., Belytschko, T., Xiao, S.P., 2004. Stable particle methods based on Lagrangian kernels. Comput. Methods Appl. Mech. Eng. 193:1035–1063. http://dx.doi.org/10. 1016/j.cma.2003.12.005. Rabczuk, T., Areias, P.M.A., Belytschko, T., 2007a. A meshfree thin shell method for nonlinear dynamic fracture. Int. J. Numer. Methods Eng. 72:524–548. http://dx.doi.org/ 10.1002/nme.2013. Rabczuk, T., Bordas, S., Zi, G., 2007b. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Comput. Mech. 40:473–495. http://dx.doi.org/10.1007/s00466-006-0122-1. Rabczuk, T., Bordas, S., Zi, G., 2010a. On three-dimensional modelling of crack growth using partition of unity methods. Comput. Struct. 88:1391–1411. http://dx.doi.org/ 10.1016/j.compstruc.2008.08.010. Rabczuk, T., Gracie, R., Song, J.-H., Belytschko, T., 2010b. Immersed particle method for fluid–structure interaction. Int. J. Numer. Methods Eng. 81:48–71. http://dx.doi.org/ 10.1002/nme.2670. Rabczuk, T., Zi, G., Bordas, S., Nguyen-Xuan, H., 2010c. A simple and robust three-dimensional cracking-particle method without enrichment. Comput. Methods Appl. Mech. Eng. 199:2437–2455. http://dx.doi.org/10.1016/j.cma.2010.03.031. Randles, P., Libersky, L., 1996. Smoothed particle hydrodynamics: some recent improvements and applications. Comput. Methods Appl. Mech. Eng. 139, 375–408. Ren, H., Zhuang, X., Cai, Y., Rabczuk, T., 2016. Dual-horizon peridynamics. Int. J. Numer. Methods Eng. 108, 1451–1476. Sengupta, A., Upadhyay, A., 2009. Locating the critical failure surface in a slope stability analysis by genetic algorithm. Appl. Soft Comput. 9:387–392. http://dx.doi.org/10. 1016/j.asoc.2008.04.015. Shi, G.-H., 1988. Discontinuous Deformation Analysis—A New Numerical Model for the Statics and Dynamics of Deformable Block Structures. (PhD Dissertation). University of California, Berkeley.

95

Shi, G.-H., 1991. Manifold method of material analysis. Transaction of the 9th Army Conference on Applied Mathematics and Computing, Minneapolis, Minnesota, USA, pp. 57–76. Song, J.-H., Areias, P.M.A., Belytschko, T., 2006. A method for dynamic crack and shear band propagation with phantom nodes. Int. J. Numer. Methods Eng. 67:868–893. http://dx.doi.org/10.1002/nme.1652. Strouboulis, T., Babuška, I., Copps, K., 2000. The design and analysis of the generalized finite element method. Comput. Methods Appl. Mech. Eng. 181:43–69. http://dx.doi. org/10.1016/S0045-7825(99)00072-9. Talebi, H., Silani, M., Rabczuk, T., 2015. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv. Eng. Softw. 80:82–92. http://dx.doi. org/10.1016/j.advengsoft.2014.09.016. Terada, K., Kurumatani, M., 2005. An integrated procedure for three-dimensional structural analysis with the finite cover method. Int. J. Numer. Methods Eng. 63: 2102–2123. http://dx.doi.org/10.1002/nme.1356. Thai, C.H., Nguyen-Xuan, H., Nguyen-Thanh, N., Le, T.H., Nguyen-Thoi, T., Rabczuk, T., 2012. Static, free vibration, and buckling analysis of laminated composite ReissnerMindlin plates using NURBS-based isogeometric approach. Int. J. Numer. Methods Eng. 91:571–603. http://dx.doi.org/10.1002/nme.4282. Thai, C.H., Ferreira, A.J.M., Bordas, S.P.A., Rabczuk, T., Nguyen-Xuan, H., 2014. Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. Eur. J. Mech. A-Solids 43:89–108. http://dx.doi. org/10.1016/j.euromechsol.2013.09.001. Thai, C.H., Nguyen-Xuan, H., Bordas, S.P.A., Nguyen-Thanh, N., Rabczuk, T., 2015. Isogeometric analysis of laminated composite plates using the higher-order shear deformation theory. Mech. Adv. Mater. Struct. 22:451–469. http://dx.doi.org/10.1080/ 15376494.2013.779050. Thai, T.Q., Rabczuk, T., Bazilevs, Y., Meschke, G., 2016. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Comput. Methods Appl. Mech. Eng. 304:584–604. http://dx.doi.org/10.1016/j.cma.2016.02.031. Vu-Bac, N., Nguyen-Xuan, H., Chen, L., Lee, C.-K., Zi, G., Zhuang, X., Liu, G.R., Rabczuk, T., 2013. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. J. Appl. Math. 2013. Vu-Bac, N., Rafiee, R., Zhuang, X., Lahmer, T., Rabczuk, T., 2015. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Compos. B Eng. 68:446–464. http://dx.doi.org/10.1016/j.compositesb.2014.09.008. Yamagami, T., Jiang, J., 1997. A search for the critical slip surface in three-dimensional slope stability analysis. Soils Found. 37:1–16. http://dx.doi.org/10.3208/sandf.37.3_1. Yamagami, T., Kojima, K., Taki, M., 1991. A search for the three-dimensional critical slip surface with random generation of surface. Proceedings of the 36th Symposium on Slope Stability Analyses and Stabilizing Construction Methods, pp. 27–34. Yang, Y.-c., Xing, H.-g., Yang, X.-g., Zhou, J.-w, 2016. Determining the critical slip surface of three-dimensional soil slopes from the stress fields solved using the finite element method. Math. Probl. Eng. 2016:11. http://dx.doi.org/10.1155/2016/7895615. Zeinkiewicz, O., Humpheson, C., Lewis, R., 1975. Associated and non-associated viscoplasticity in soils mechanics. Geotechnique 25, 671–689. Zhang, X., 1988. Three-dimensional stability analysis of concave slopes in plan view. J. Geotech. Eng. 114, 658–671. Zheng, H., Xu, D.D., 2014. New strategies for some issues of numerical manifold method in simulation of crack propagation. Int. J. Numer. Methods Eng. 97:986–1010. http://dx. doi.org/10.1002/nme.4620. Zheng, H., Tham, L.G., Liu, D., 2006. On two definitions of the factor of safety commonly used in the finite element slope stability analysis. Comput. Geotech. 33:188–195. http://dx.doi.org/10.1016/j.compgeo.2006.03.007. Zheng, W., Zhuang, X., Tannant, D.D., Cai, Y., Nunoo, S., 2014. Unified continuum/ discontinuum modeling framework for slope stability assessment. Eng. Geol. 179, 90–101. Zhuang, X.Y., Augarde, C., Bordas, S., 2011. Accurate fracture modelling using meshless methods, the visibility criterion and level sets: formulation and 2D modelling. Int. J. Numer. Methods Eng. 86:249–268. http://dx.doi.org/10.1002/nme.3063. Zhuang, X., Augarde, C.E., Mathisen, K.M., 2012a. Fracture modeling using meshless methods and level sets in 3D: framework and modeling. Int. J. Numer. Methods Eng. 92:969–998. http://dx.doi.org/10.1002/nme.4365. Zhuang, X., Heaney, C., Augarde, C., 2012b. On error control in the element-free Galerkin method. Eng. Anal. Bound. Elem. 36, 351–360. Zhuang, X., Huang, R., Liang, C., Rabczuk, T., 2014a. A coupled thermo-hydro-mechanical model of jointed hard rock for compressed air energy storage. Math. Probl. Eng. 2014. Zhuang, X.Y., Cai, Y.C., Augarde, C., 2014b. A meshless sub-region radial point interpolation method for accurate calculation of crack tip fields. Theor. Appl. Fract. Mech. 69: 118–125. http://dx.doi.org/10.1016/j.tafmec.2013.12.003.