A two-grid search scheme for large-scale 3-D finite element analyses of slope stability

Computers and Geotechnics 62 (2014) 203–215

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Technical Communication

A two-grid search scheme for large-scale 3-D finite element analyses of slope stability Xi Chen a,⇑, Yongkang Wu b, Yuzhen Yu b, Jiankun Liu a, Xi Frank Xu a, Jun Ren a a b

School of Civil Engineering, Beijing Jiaotong University, China State Key Laboratory of Hydroscience and Engineering, Tsinghua University, China

a r t i c l e

i n f o

Article history: Received 27 February 2014 Received in revised form 20 July 2014 Accepted 20 July 2014

Keywords: Slope stability Shear strength reduction finite element method (SSRFEM) Factor of safety Generalised bisection search algorithm Two-grid scheme Undrained analysis

a b s t r a c t It is well known that the trial process for seeking the safety factor in the shear strength reduction finite element method (SSRFEM) is quite expensive, particularly for large 3-D slope stability analyses. The search algorithm for the safety factor is crucial to the entire solution process for the shear strength reduction finite element method, but few studies have attempted to exploit it. Among search algorithms, the commonly used bracketing and bisection search has not been fully optimised. Consequently, to improve the search scheme for the safety factor associated with the shear strength reduction finite element method, two strategies are suggested. First, a generalised bisection search algorithm is proposed to reduce the possibility of encountering non-convergence from a statistical point of view. To further improve the efficiency, a new two-grid scheme, characterised by a coarse mesh search and followed by a fine mesh search, is developed. Based on the drained or undrained analyses of the 3-D slope examples, the new search algorithm can markedly outperform the commonly used bisection search algorithms based on a single finite element mesh. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Geotechnical stability is a classical and long-standing problem in the geotechnical field, and according to Sloan [1], four typical methods are generally used to perform a geotechnical stability analysis. Conventionally, slope stability analysis has resorted to limit-equilibrium methods (e.g. [2]), but the shear strength reduction finite element method (SSRFEM) has become appealing only recently following advances in computer technologies. The limit equilibrium methods (LEM) possess the advantages of simplicity and effectiveness, but the FE-based slope stability analysis methods may provide more information (such as the stress and deformation fields) apart from the safety factor. One of the earliest ideas for strength reduction is credited to Zienkiewicz et al. [3], but little attention had been given to it until the works by Ugai [4], Matsui and San [5], Dawson et al. [6], and Griffiths and Lane [7]. Because of these developments, the potential for the SSRFEM method has been widely exploited, but it was initially restricted

⇑ Corresponding author. Address: School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China. Tel.: +86 10 51684867; fax: +86 10 51683764. E-mail addresses: [email protected] (X. Chen), [email protected] (Y. Wu), [email protected] (Y. Yu), [email protected] (J. Liu), xixu@m. bjtu.edu.cn (X.F. Xu), [email protected] (J. Ren). http://dx.doi.org/10.1016/j.compgeo.2014.07.010 0266-352X/Ó 2014 Elsevier Ltd. All rights reserved.

to 2-D simulations (e.g. [8]). The 3-D finite element analysis, which had been hindered in past decades by the limitations of computer hardware and solution technologies, only recently became practical for geotechnical applications [9,10]. For example, Cai and Ugai [11], Zheng et al. [12], Griffiths and Marquez [13], Wei et al. [14], Huang and Jia [15] and Wei et al. [16] studied slope stability within the 3-D finite element framework. To satisfy rapidly growing demands, the commercial FEM software PLAXIS extends the facility of shear strength reduction from 2-D to 3-D simulations (e.g. [17]). According to these studies, the advantages of the SSRFEM method for slope stability over the limit-equilibrium approaches have been very well recognised. For instance, in the finite element analysis, no assumptions need to be made about the location or shape of the slip surface and the slice side forces versus the limitequilibrium methods. Some classical elasto-plastic soil models, such as the Mohr–Coulomb model and the Drucker–Prager model, can be readily incorporated. The complicated geometry, boundary and loading conditions can be easily treated, and the stress/ deformation fields can be progressively monitored. The most important benefit is that the SSRFEM method provides a more realistic analysis tool for complicated slopes, which explains why some researchers believe that when combined with the shear strength reduction (SSR) technique, the FEM is a powerful tool for slope stability analysis (e.g. [18]).

204

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Based on our experience, the problems associated with the SSR technique, including the high computational cost of the SSR technique and the very limited soil models that are appropriate for it, remain unsolved, which hinders its applications. Originally, the SSR technique was proposed for the Mohr–Coulomb and Drucker–Prager models, whereas other models, such as the Hoek–Brown model, were to be transformed to suit this technique (e.g. [19–21]). Hence, this work aims to improve the efficiency of the SSR technique to make it realistic and appealing for large-scale 3-D geotechnical applications, such as slope stability analysis. When implementing the SSRFEM method, Smith and Griffiths [22] evaluated a series of trial strength reduction factors to seek the factor of safety (FOS). However, providing this sequence of factors can be difficult for inexperienced users. Hence, an initial search range for the FOS should be prescribed. Within the predefined search range, the question of how to search for the FOS effectively becomes crucial because the search algorithm can significantly influence the accuracy of the FOS and the entire computational cost. Among search algorithms, the bracketing and bisection search algorithm is most frequently used (e.g., [6,23– 26]). Another search method, which was mentioned by Won et al. [27], is called the incremental refining search. These methods have been frequently utilised but have been subjected to further improvement because their inherent characteristics have not been fully exploited. For example, the nonlinear iterations’ distribution in the search range is uneven; that is, the number of nonlinear iterations in the convergence range is usually small, while the number of nonlinear iterations in the non-convergence range is the maximum allowable number (such as 1000 or 500). Without considering this characteristic, the bracketing and bisection search method sets equal weights on the convergence and nonconvergence ranges. This paper is thus organised as follows. Section 2 briefly introduces the SSRFEM method using the Mohr–Coulomb material model. In Section 3, the bracketing and bisection search algorithm is reviewed first, and then a generalised bisection search algorithm is proposed with a bias towards the convergence range so the proposed search method can appropriately reduce the likelihood of encountering non-convergence. In Section 4, several schemes used for the preliminary search bound of the FOS are introduced, and a new two-grid scheme is emphasised. Section 5 introduces several 3-D slope examples upon which the newly proposed search algorithm is assessed, and the proposed algorithm is compared with the commonly used ones. Finally, concluding remarks are given in Section 6. 2. Shear strength reduction finite element method In numerical simulations, geotechnical failures can be triggered by either an increase in the external loads or a decrease in the geomaterial strength, which lead to the incremental loading method and shear strength reduction method, respectively, as

τ

σ3 Real Mohr circle

(a)

sf ¼ c0 þ rn tan /0

σ1

Mohr circle at j-th step of gravity loading

ð1Þ 0

0

where the cohesion c and the internal friction angle / are two effective shear strength parameters. With the SSR technique, the shear strength is decreased by reducing the values of the two shear strength parameters, that is,

cr ¼ c0 =SRF c

ð2aÞ

tan /r ¼ tan /0 =SRF /

ð2bÞ

where SRFc and SRF/ are the strength reduction factors corresponding to the cohesion and the internal friction angle, respectively, and it is usually assumed that SRFc = SRF/ = SRF. When applying the shear strength reduction method to geotechnical problems such as slope stability analyses, three behaviours can be utilised to identify the occurrence of slope failure, namely, sudden substantial changes in the displacement of various marked nodes, a connected plastic shear band, and a non-equilibrium state. Because all three slope behaviours result in the non-convergence of the nonlinear iterative method, assessing the convergence or non-convergence of the nonlinear iterative method, which may be utilised independently or in conjunction with other behaviours/criteria, is an effective approach for judging slope stability in finite element analyses. Consequently, given the predefined maximum number of iterations maxit_nl and the convergence tolerance tol_nl, the non-convergence of the nonlinear iterative method is employed in the present study to identify slope failures, which must be carefully assessed to ensure that they are produced by a genuine geotechnical failure instead of artificial numerical effects (e.g. [1]). When an iterative linear solver is employed, the non-convergence of the nonlinear iterative method may be attributed to the non-convergence of the iterative linear solver. In that case, the direct linear solver can be employed to rule out this possibility. Evidently, SRF at the transition point from convergence to non-convergence can be defined as the FOS, that is, FOS = SRFf. When increasing the SRF, the transition point from convergence to non-convergence is defined as the FOS. The bisection algorithm was originally proposed for finding the roots of a nonlinear equation, f(x) = 0 (e.g. [29]), as shown in Fig. 2a. Nonetheless, the commonly utilised bisection search algorithm is not suitable for the SSRFEM method. From Fig. 2a and b, a significant difference can be observed when the bisection algorithm is applied to the two problems. That is, the evaluation cost for f(x) > 0 or f(x) < 0 is almost the same in the root finding problems, while in the SSRFEM method, the evaluation cost at the convergence side is remarkably smaller than that at the non-convergence side. Consequently, the

τ

MC failure envelope

φ0 c0

demonstrated in Fig. 1a and b. Compared with the incremental loading method, the shear strength reduction method has received a wider range of applications (e.g. [28]). When the Mohr–Coulomb failure criterion is used, the shear strength is expressed as

Real MC failure envelope

φ0

σ

c0

σ3

σ1

Reduced MC failure envelope

Mohr circle at failure Real Mohr circle

(b) Fig. 1. Evolution of stress circle (a) incremental loading method; (b) shear strength reduction method.

σ

205

Nonlinear iterations in each trial

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Fig. 2. Utilization of bisection algorithm for (a) root finding of f(x); (b) FOS search in SSRFEM.

conventional bisection algorithm for the safety factor associated with the SSRFEM method can be improved in two aspects: (1) A bias on the convergence side should be reflected in the new search algorithm. (2) Because the evaluation of the equilibrium state is usually expensive, a cost-effective approach should be designed to assess the equilibrium state, which is equivalent to finding a function g(x) close to f(x), but g(x) is apparently less costly to evaluate than f(x). In the following sections, the search algorithm for the safety factor in the SSRFEM method is designed based on these two aspects. 3. Generalised bisection algorithm for the shear strength reduction finite element method 3.1. Conventional bracketing and bisection search algorithm As mentioned in [29], there are two distinct phases in finding the roots of a nonlinear equation: bounding the roots and refining the roots. In the slope stability analysis, the search scheme is used to find FOS = s in a range [a0, b0] that either the user prescribes or is automatically generated so that the adopted nonlinear iterative method converges in the range [a0, s] but not in the range (s, b0], according to the user-specified maximum number of iterations. To find the switching point s, the bracketing and bisection search algorithm is defined concisely as (

n1 an ¼ an1 þb ; bn ¼ bn1 ; if nonlinear iteration converges at 2 n1 an ¼ an1 ; bn ¼ an1 þb ; if nonlinear iteration diverges at 2

an1 þbn1 2 ðn P 1Þ an1 þbn1 2

ð3Þ

where the trial is indexed by n. Because the strength reduction factor is defined as SRFn = (an1 + bn1)/2, the reduced shear strength is expressed as

cn ¼ c0 =SRF n

ð4aÞ

tan /n ¼ tan /0 =SRF n

ð4bÞ

Convergence is achieved once the following criterion is satisfied,

bn  an ¼ ðb0  a0 Þ

1 6 2n



it nl tot  maxit nl  n  50% þ it nl av e  n  50%

ð5Þ

ð6Þ

Here, ceil() denotes rounding a number towards positive infinity. When Eq. (5) is satisfied and the trial process is terminated, the

ð7Þ

where maxit_nl is the prescribed maximum number of nonlinear iterations, and it_nl_ave is the average number of nonlinear iterations for the convergence cases. 3.2. A generalised bisection search algorithm As shown in Fig. 2b, the evaluation cost at the convergence side (i.e. SRF 2 [a0, s]) is substantially smaller than that at the non-convergence side (i.e. SRF 2 (s, b0]) in the SSRFEM method, indicating that taking equal weights for the two search regions, as occurs in the conventional bisection algorithm, is not optimal. To have a bias on the convergence region, a generalised bisection search algorithm is proposed:

8 an ¼ an1 ð1  aÞ þ bn1 a; bn ¼ bn1 ; > > > < if nonlinear iteration converges at a ð1  aÞ þ b a n1 n1 ðn P 1Þ > a ¼ a ; b ¼ a ð1  a Þ þ b a ; n n1 n n1 n1 > > : if nonlinear iteration diverges at an1 ð1  aÞ þ bn1 a ð8Þ where a 2 (0, 1/2] is a weighting factor; when a = ½, it reduces to the conventional bisection algorithm. In the generalised bisection algorithm, there are two limit cases:

ðb0  a0 Þan 6  for the fastest shrinking

ð9aÞ

ðb0  a0 Þð1  aÞn 6  for the slowest shrinking

ð9bÞ

Obviously, FOS  a0 corresponds to the fastest shrinking case, and FOS  b0 corresponds to the slowest shrinking case. The required numbers of trials for the two bounds are thus estimated, respectively, by

  ðb0  a0 Þ ðaÞ for the fastest shrinking n P n ¼ ceil log1a



ð1aÞ

n P n

where  is the prescribed convergence tolerance. Obviously, to satisfy the convergence criterion, the total number of trials has to be

  ðb0  a0 Þ n P n ¼ ceil log2

factor of safety is obtained by setting FOS = SRFn. As mentioned previously, for each trial, the probability of encountering convergence is the same as the probability of encountering non-convergence for the nonlinear iterative method. Hence, the total consumed nonlinear iterations is estimated as

 ¼ ceil log

ðb0  a0 Þ 1 1a



ð10aÞ

 for the slowest shrinking ð10bÞ

After exiting the iteration, the factor of safety is obtained by setting FOS = SRFn. It is also clear that the slowest and fastest shrinking cases correspond to the upper and lower bounds of the required total nonlinear iteration number, respectively, and the two bounds are estimated as ð1aÞ

it nl tot  maxit nl  n

for the upper bound

ð11aÞ

206

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215 ðaÞ

it nl tot  it nl av e  n

for the lower bound

ð11bÞ

As a result, the required total nonlinear iteration number is estimated as ð1aÞ

it nl tot  maxit nl  n

ðaÞ

 a þ it nl av e  n  ð1  aÞ

ð12Þ

The key issue is thus finding an optimal value of a that can be used to minimise the total number of required nonlinear iterations. A small a (e.g. 0.1) obviously increases the number of trails. The conventional bisection search method has the fastest speed to shrink the search range, but it puts equal weights on the convergence and non-convergence sides. In the generalised bisection search algorithm, a may not remain constant during the search process. From Fig. 2b, a close observation of the characteristics of the required nonlinear iterations versus the SRF indicates that a sudden increase in the nonlinear iterations only appears in a narrow range of the SRF, indicating that the search process can be divided into two phases. Hence, a, which is utilised in Phase 1 of the search process, can be defined differently from that utilised in Phase 2; that is, a is allowed to switch to the value of 0.5 once Phase 2 is accessed. However, for simplicity, a remains constant during the search process. For more details, see Algorithm 1 in Appendix A. 3.3. Bracketing and incremental refining search algorithm Differing from the bracketing and bisection search algorithm, it is possible to use a fixed division, which leads to the following bracketing and incremental refining search algorithm (e.g. [27]). In the incremental refining search method, if the nonlinear iteration converges at ai1 + (bi1  ai1)jb, then the search range [ai, bi] is not changed but the trial point is moved to the next point:

ai1 þ ðbi1  ai1 Þðj þ 1Þb;

ðj ¼ 1; . . . ; 1=b  2Þ

ð13aÞ

If the nonlinear iteration does not converge at ai1 + (bi1  ai1)jb, then the two limits are updated by

ai ¼ ai1 þ ðbi1  ai1 Þðj  1Þb; bi ¼ ai1 þ ðbi1  ai1 Þjb; ðj ¼ 1; . . . ; 1=b  1Þ

ð13bÞ

where b 2 (0, 1/2]. Index i represents the updating level of the two limits, and index j denotes the uniform division in the search range. It is interesting to note, however, that the incremental refining search algorithm in which b = 1/2 reduces to the bisection algorithm. 4. A two-grid scheme for the FOS search algorithm Because the evaluation of the equilibrium state in a large 3-D finite element framework is usually expensive, an efficient approach should be designed to assess the equilibrium state. Here, a new search scheme is proposed by decomposing the search process into two phases (or stages), as shown in Fig. 2b. It is well known that the factor of safety can be affected by the mesh size, element type, and dimensions. To rapidly determine the preliminary bounds for the safety factor, three possible approaches can be attempted: (1) Switching from a 3-D finite element mesh to a 2-D mesh. (2) Switching from high-order finite element types to low-order element types. (3) Switching from a fine finite element mesh to a coarse mesh (i.e. a two-grid scheme). As we observe, all three approaches are intended to rapidly find the preliminary bounds for the safety factor. The basic idea behind

the three approaches is analogous to solving a cost-effective function g(x) to approximately replace f(x). Although Griffiths and Marquez [13] compared the results from the 2-D and 3-D SSRFEM analyses and noted that the FOS obtained from the 2-D models is usually more conservative than that obtained from the 3-D models, the observation may not be guaranteed for complex slope examples. Hence, using a 2-D mesh to determine the lower bound of the safety factor for the 3-D FEM analyses may not be wise. There is a certain degree of similarity between the multi-grid iterative method and the present two-grid (or coarse–fine mesh) method. In the multi-grid iterative method, the coarse mesh is utilised to find nodal solutions whose interpolations are employed as the initial guess of an iterative scheme based on a fine mesh (e.g. [30]), whereas in the present two-grid scheme, the coarse mesh is utilised to narrow the search range of the FOS so that fewer trials will be involved in the fine mesh search. Here, the slope example in Griffiths and Marquez [13] is employed to investigate the change of the FOS with mesh refining. Fig. 3 shows a 3-D finite element mesh and its geometry. For this slope, the FOS computed using the finite element method decreases with mesh refining as shown in Fig. 4, but in [13], the FOS is 1.734, based on the coarse finite element mesh that is composed of 550 20-node hexahedral elements. In this example, the FOS that varies with mesh refining is obtained by implementing the SSRFEM method in ABAQUS (e.g. [31]). A decrease in the FOS with 3-D mesh refining has also been observed in other studies (e.g. [13]). Because the two-grid scheme is more general than the scheme of switching from a 3-D mesh to a 2-D mesh, it is employed to determine the search range of the safety factor for a later fine mesh search. As a result, the renewed search range of the safety factor is defined as [FOS determined by coarse mesh  dFOS, FOS determined by coarse mesh], where dFOS is a predefined tolerance to ensure that the FOS computed by the 3-D fine mesh is located in the range. 5. Numerical examples The shear strength reduction finite element program for the 3-D slope stability analysis is coded in FORTRAN, and the two-grid scheme for the FOS is automatically accomplished with the coarse mesh search followed by the fine mesh search. To apply the SSRFEM with two-grid search scheme, the user-defined parameters include maxit_nl and tol_nl for terminating the nonlinear iterations, the initial search range [a0, b0] for bracketing FOS, the tolerance  for terminating the bisection search process and dFOS for updating the search range of FOS. All of these parameters are generally prescribed by the program and can be changed only by experienced users before executing FE computations. In addition, it must be emphasised that all of the numerical experiments conducted in this work are based on ordinary computer architectures. 5.1. Mathematical simulations To investigate the performance of the search algorithms, a function of SRF is assumed to be (

  f ðSRFÞ ¼ INT maxit nl  ½ðSRF  a00 Þ=ðFOS  a00 Þp ; a00 < a0 ; SRF 6 FOS SRF > FOS f ðSRFÞ ¼ maxit nl; ð14Þ

where INT() is an operator for rounding numbers, f(SRF) is a function of SRF for the number of nonlinear iterations, as shown in Fig. 2b, and p is the power, which is assumed to be 100 according to the curve of nonlinear iterations versus the SRF from the practical slope stability analysis. Fig. 5 shows the number of accumulative nonlinear iterations of the three search algorithms with a search tolerance of  = 0.001

207

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

10m

10m

20m

H=10m 26.57° 5m

(b) Soil shear strength:φ u= 0°, cu = 0.2γ H.

(a)

Fig. 3. Finite element mesh and geometry of Griffiths and Marquez’s slope example [13].

FOS

achieves time savings of 25%, 32%, 28%, and 25% for FOS = 0.6, 1.2, 2.2 and 3.2, respectively, compared with the conventional bisection algorithm. 5.2. Three-dimensional slope stability analysis using the SSRFEM technique In the finite element computational framework, the Newton– Raphson nonlinear iteration is expressed as

Rðuði1Þ Þ þ

Fig. 4. FOS computed by finite element method varying with mesh refining.

Accumulative nonlinear iterations

@ UðuÞ ¼ Fext  Fint ðuÞ @u

ð16Þ

where Fint(u) and Fext are the internal and external force vectors, respectively. The gradient of R about u is the tangential stiffness matrix. To terminate the iterative process of the Newton–Raphson method, an appropriate convergence criterion, such as the displacement-based criterion or the out-of-balance force based criterion,

Generalized bisection (α ) Incremental refining search ( β) Coventional bisection ( α=0.5 )

20000 16000

12000

12000

8000

8000

FOS=0.6

4000 0

0.1

(a) Accumulative nonlinear iterations

RðuÞ ¼

16000

0.2 0.3 α (or β)

0.4

FOS=1.2

4000 0

0.5

0.1

(b)

20000

20000

16000

16000

12000

12000

8000

0.2 0.3 α (or β)

0.4

0.5

8000

FOS=2.2

4000 0

(c)

ð15Þ

where the superscript i indexes the nonlinear iteration in each load increment, u is the displacement vector, and R is the residual force or the out-of-balance force vector derived from the total potential energy U:

versus the a (or b) value for the four FOS. Apparently, the proposed generalised bisection search algorithm is generally superior to the incremental refining algorithm, and it outperforms the conventional bisection search algorithm for a wide range of a. In the range a 2 [0.06, 0.3], the generalised bisection search algorithm leads to savings of 8%, 20%, 15% and 14% in solution time for FOS = 0.6, 1.2, 2.2 and 3.2, respectively, compared with the conventional bisection search algorithm. When the search tolerance is set as  = 0.01, as shown in Fig. 6, the advantage of the generalised bisection search algorithm over the other two algorithms is more evident, and in the range a 2 [0.06, 0.3], the proposed method

20000

@R ði1Þ ðu ÞDuðiÞ ¼ 0 @u

0.1

0.2 0.3 α (or β)

0.4

FOS=3.2

4000

0.5

0

(d)

0.1

0.2 0.3 α (or β)

0.4

0.5

Fig. 5. Accumulative nonlinear iterations for generalised bisection search and incremental refining search with  = 0.001 versus a (or b) when (a) FOS = 0.6; (b) FOS = 1.2; (c) FOS = 2.2; (d) FOS = 3.2.

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Accumulative nonlinear iterations

208

10000

12000 10000

8000

8000

6000

6000

4000

4000

FOS=0.6

2000

(a) Accumulative nonlinear iterations

Generalized bisection (α ) Incremental refining search ( β) Coventional bisection ( α=0.5 )

12000

0

0.1

0.2 0.3 α (or β)

0.4

0.5

(b)

16000

16000

12000

12000

8000

8000

4000

0

0.1

0.2 0.3 α (or β)

0.4

0.5

4000

FOS=2.2

0

(c)

FOS=1.2

2000

0

0.1

0.2 0.3 α (or β)

0.4

FOS=3.2

0 0

0.5

0.1

(d)

Fig. 6. Accumulative nonlinear iterations for generalised bisection search and incremental refining search with FOS = 2.2; (d) FOS = 3.2.

should be selected. In this work, the Newton–Raphson iteration is terminated when the out-of-balance force is sufficiently small:

RðuðiÞ Þ 6 tol nl or maxit nl is exceeded Fext

ð17Þ

Here, tol_nl and maxit_nl are the prescribed tolerance and the allowable maximum nonlinear iteration count, respectively. To solve the linear system of equations arising from each Newton–Raphson iteration, a direct linear solver or an iterative linear solver can be employed. However, for a large 3-D finite element analysis, iterative linear solvers are usually more efficient. When only associated soil plasticity is involved, the symmetric successive over-relaxation (SSOR) preconditioned conjugate gradient (PCG), in conjunction with the Eisenstat trick (e.g. [32]), is employed to solve the resultant symmetric positive definite (SPD) linear system of equations, but a non-symmetric iterative linear solver must be utilised if non-associated soil plasticity is considered (e.g. [33]). When the iterative linear solver is applied, the relative residual criterion, with a tolerance of 1.e  6, is used to terminate its iterative process (e.g. [34]). 5.2.1. Three-dimensional slope stability analysis with/without an external surface load According to [17], the safety factor of a slope without an external surface load is 1.57 (1.56 when using the Bishop slip circle method), but it is 1.25 when the external load p = 30 kPa is applied. To investigate the single mesh scheme and the two-grid scheme, a series of finite element meshes with 270, 680, 1587 or 3744 20node hexahedral elements is generated (Fig. 7 presents the finest finite element mesh and its geometry for the slope example). First, maxit_nl = 50 and tol_nl = 0.001 are adopted for the Newton–Raphson nonlinear iteration, and the generalised bisection search algorithm is assessed with and without the two-grid scheme, as shown in Table 1. a = 0.5, which is based on a single mesh (nels = 3744), corresponds to the conventional bisection search algorithm; its solution time (i.e. 183.9 min and 163.1 min for the cases without and with external load, respectively) is used here as a benchmark. Table 1 indicates that the generalised bisection search algorithm with a = 0.2 leads to a savings in solution time of 33% and 4% for the cases without and with external

0.2 0.3 α (or β)

0.4

0.5

 = 0.01 versus a (or b) when (a) FOS = 0.6; (b) FOS = 1.2; (c)

load, respectively, and confirms the observation made in Section 5.1. To evaluate the slope stability with the finest mesh (nels = 3744), three coarse meshes with 270, 680, and 1587 elements are provided, leading to the two-grid scheme with mesh ratios Rm = 7.2%, 18.2% and 42.4%, respectively. When the two-grid scheme is utilised in the search algorithm, a further improvement of the total solution time can be achieved. As noticed, the performance of the two-grid scheme may generally degenerate with an increase in Rm. For example, when Rm = 42.4%, the solution time of the SSRFEM method with the two-grid scheme becomes even longer than that observed for a single mesh in the case without an external load. When the two-grid scheme is utilised, the total solution time of the SSRFEM method comprises two parts: the solution time for the fine mesh search and the solution time in the bracket for the coarse mesh search. With the smallest mesh ratio, Rm = 7.2%, a more remarkable savings in the total solution time, 62% and 47%, is attained using the two-grid scheme. In Table 1, the accumulative nonlinear iterations and the ratio of times from non-convergence to total trials, i.e., the ratio of the non-convergence trial, are also investigated to show their effects on the search algorithm. The ratio of the non-convergence trial may further be used to evaluate the performance of the search algorithms because many non-convergence trials may significantly increase the accumulative nonlinear iterations. The accumulative nonlinear iterations can be viewed as a good indicator in a computational cost comparison because their accumulation is proportional to the total solution time required for the SSRFEM method. To be complete, these indicators for the coarse mesh search are also given in brackets. As can be observed, the low ratio of times from non-convergence to total trials contributes substantially to the efficiency of the two-grid scheme in the search algorithms for the FOS. The renewed search range of the FOS is provided to show that the use of [FOS determined by coarse mesh  dFOS, FOS determined by coarse mesh] with dFOS = 0.1 should be adequate for the fine mesh search. To be more certain, however, a larger value (such as 0.2) can be used. Furthermore, maxit_nl = 500 is used to evaluate the generalised bisection search algorithm with other parameters kept unchanged, as shown in Table 2. The solution time spent by the single mesh (nels = 3744) with a = 0.5 is 1335.8 min and 978.4 min for the cases

209

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Fig. 7. The finest mesh with 3744 elements and geometry of embankment investigated in [17].

Table 1 Generalised bisection search algorithms with and without two-grid scheme, respectively, for slope stability analysis (maxit_nl = 50, tol_nl = 0.001) providing that [a0, b0] = [0.5, 4.0],  = 0.01 and dFOS = 0.1. Two-grid scheme, Rm

3744/3744

3744/3744

270/3744

680/3744

1587/3744

a

0.5

0.2

0.2/0.2

0.2/0.2

0.2/0.2

Without external load FOS Renewed search range of FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.5544 – 6/12 421 1.000

1.5536 – 3/12 288 0.669

1.5513 [1.4906, 1.5906] 1/7 (3/13) 169 (301) 0.363 (0.022)

1.5527 [1.4688, 1.5688] 1/11 (2/16) 235 (290) 0.514 (0.072)

1.5543 [1.4649, 1.5649] 2/14 (3/16) 365 (396) 0.786 (0.294)

With external load FOS Renewed search range of FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.2332 – 5/12 378 1.000

1.2325 – 5/10 363 0.961

1.2336 [1.1588, 1.2588] 1/9 (3/10) 207 (241) 0.509 (0.020)

1.2312 [1.1473, 1.2473] 1/11 (3/10) 202 (298) 0.487 (0.086)

1.2336 [1.1404, 1.2404] 0/13 (5/10) 228 (374) 0.543 (0.329)

Table 2 Generalised bisection search algorithms with and without two-grid scheme, respectively, for slope stability analysis (maxit_nl = 500, tol_nl = 0.001) providing that [a0, b0] = [0.5, 4.0],  = 0.01 and dFOS = 0.1. Two-grid scheme

3744

3744

270/3744

680/3744

1587/3744

a

0.5

0.2

0.2/0.2

0.2/0.2

0.2/0.2

Without external load FOS Renewed search range of FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.5613 – 5/12 2936 1.000

1.5604 – 4/13 2343 0.800

1.5564 [1.5039, 1.6039] 1/7 (3/17) 696 (2615) 0.229 (0.026)

1.5613 [1.4779, 1.5779] 2/12 (3/11) 1743 (1906) 0.583 (0.068)

1.5619 [1.4688, 1.5688] 0/13 (2/16) 903 (1343) 0.293 (0.147)

With external load FOS Renewed search range of FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.2349 – 4/11 2161 1.000

1.2358 – 3/8 1717 0.843

1.2350 [1.1665, 1.2665] 1/8 (4/11) 690 (2975) 0.319 (0.041)

1.2367 [1.1473, 1.2473] 2/14 (3/10) 1791 (1644) 0.846 (0.079)

1.2368 [1.1436, 1.2436] 0/13 (4/10) 636 (2405) 0.290 (0.361)

without and with external load, respectively. Similarly, in the generalised bisection search algorithm, a = 0.2 is still better than a = 0.5 from the perspective of the solution time. Compared with the results obtained using maxit_nl = 50, a more remarkable savings in the solution time can be achieved using the two-grid scheme with Rm = 7.2%; in the cases without and with external

load, the time savings are 75% and 64%, respectively. According to Eq. (6), the required number of trials for the conventional bisection search is at least n ¼ 9, which is generally less than the realistic trials presented in Table 1 or Table 2 because the condition of ðb  aÞ <  and the convergence trial for the Newton–Raphson must be satisfied concurrently in Algorithm 1.

210

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

5.2.2. Slope stability analysis without and with a weak layer An interesting slope example with a thin weak layer was studied by Griffiths and Lane [7], but this example was in 2-D space. Here, the same problem is evaluated in 3-D space, with its geometry and finite element meshes (for L/H = 2) consisting of 336 and 3784 20-node hexahedral elements, as shown in Fig. 8. As a result, the mesh ratio in the two-grid scheme is Rm = 9.7% when the two-grid scheme is applied. The undrained clay slope has the shear strength /u = 0°, cu1 = 0.25 csH = 25 kPa (cs = 20 kN/m3 and H = 5 m), and the embedded thin layer is assessed in three cases: cu2/cu1 = 1.0, cu2/cu1 = 0.6 and cu2/cu1 = 0.2. The undrained analyses of the geotechnical problems are generally more challenging than the drained analyses and should be treated carefully because the numerical difficulty that is associated with the undrained analysis is more severe for the iterative linear solvers than it is for the direct linear solvers (e.g. [35]). In the present study, two approaches were attempted and compared with each other in the undrained analysis. The first is called the penalty method, which produces the 1  1 block system:

Kud u ¼ RðuÞ with Kud ¼

XZ e

Ve

BTu Dud Bu dV



(b)

H H

(c)

2H

cu2 cu1

2H

2

cu1

1

2 1

1.2H

0.6H

0.4H

1 1 0.4H Shear strength of clay: φ u = 0°, cu1 = 0.25γ s H.

Fig. 8. Griffiths and Lane’s slope example [7]. (a) Fine mesh with 3784 elements; (b) coarse mesh with 366 elements; (c) geometry and profile of the slope.

ð19Þ

Kw Dev n

with Dev ¼ mT De

ð20Þ

where e = {e11, e22, e33, e12, e23, e31}T and ev is the volumetric strain. The second approach is called the Lagrange multiplier method or the mixed-form method and results in the 2  2-coupled block system



K

L

LT

0



Du Dpw



 ¼

Df 0

ð21Þ

where

K¼

XZ Ve

e

XZ e

(a)

2H 1.2H

pw;i ¼ pw;i1 þ

L¼

L/2 = 5m

Kw mmT n

where Dep is the elasto-plastic constitutive matrix; Kw is the bulk modulus of water; n is the porosity of soil; and m = {1 1 1 0 0 0}T is the Kronecker delta in vector form. At each integration point, the excess pore water pressure pw is updated by

ð18Þ

where Bu is the strain–displacement matrix, and

0.6H

Dud ¼ Dep þ

Ve

BTu Dep Bu dV

BTu mNp dV

 ð22aÞ

 ð22bÞ

where Np is the shape function vector for the fluid elements. Obviously, the 2  2-coupled block system in Eq. (21) is readily derived from Biot’s coupled linear system of equations by setting Dt = 0, which is reasonable because setting Dt = 0 is equivalent to a very fast loading rate. In the following studies, both approaches are utilised and evaluated. The numbers of unknowns in the 1  1 block system and the 2  2-coupled block system arising from the fine mesh are 47,095 and 51,847 (i.e. 47,095 + 4752), respectively. As shown in Table 3, the FOS computed by Rocscience Inc. [24] are given as a reference, but these FOS were attained in a 2-D plane strain analysis with a 6-node triangular element type (i.e. T6) and an 8-node quadrilateral element type (i.e. Q8), respectively. In Table 3, 3-D SSRFEM computations are performed based on the penalty method and the 1  1 block system in Eq. (18), and all of the FOS computed are generally higher than those given by [24]. The single grid scheme with the conventional bisection search algorithm may lead to a solution time of 587.8 min, 801.1 min and 580.6 min for cu2/cu1 = 1.0, 0.6 and 0.2, respectively. The generalised bisection search algorithm with a = 0.2 can result in a slight savings in the total solution time compared with the conventional algorithm (which corresponds to a = 0.5), but significantly greater savings are achieved by the two-grid scheme depending on the ratio of cu2/cu1. Because of the preliminary search range determined by the coarse mesh search, the trials associated with the fine mesh search can be reduced by approximately 40–60%. Nevertheless, when comparing the single-grid search algorithm (a = 0.2) with the two-grid search algorithm, it is found that the two-grid search algorithm may not always outperform the single-grid search algorithm in this study case. A close observation of the indicators, such as the ratio of non-convergence trial and the accumulative nonlinear iterations, discloses that for cu2/cu1 = 1.0 and cu2/cu1 = 0.6, the slight savings attained in the solution time because of slightly less accumulative nonlinear iterations can be overwhelmed by the cost of the coarse mesh search. Fig. 9 illustrates the search process of the generalised bisection algorithm for the FOS. It is clearly observed that in the two-grid scheme, the search process for the FOS is decomposed into two stages: the coarse mesh search and the fine mesh search. The non-convergence trial points are denoted by filled symbols to illustrate the possibility of non-convergence trials in each search algorithm. The probability that non-convergence trials will occur in the conventional bisection search is approximately 50%, while it is gener-

211

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Table 3 Performance of generalised bisection search algorithms with and without two-grid scheme, respectively, for slope stability analysis based on the penalty method (given maxit_nl = 50, tol_nl = 0.001, [a0, b0] = [0.2, 3.0],  = 0.01, dFOS = 0.1 and Kw/n = 500,000 kPa). cu2/cu1

1.0

0.6

0.2

1.45 (T6) 1.47 (Q8)

1.35 (T6) 1.35 (Q8)

0.62 (T6) 0.59 (Q8)

Single grid scheme, a = 0.5 FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.5084 3/12 311 1.000

1.4318 6/12 430 1.000

0.7804 7/13 536 1.000

Single grid scheme, a = 0.2 FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.5068 3/15 292 0.910

1.4309 3/13 310 0.717

0.7706 3/9 372 0.657

Two-grid scheme (Rm = 366/3784), a (0.2/0.2) FOS Renewed search range of FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.5081 [1.4946, 1.5946] 1/7 (4/12) 271 (407) 0.903 (0.059)

1.4329 [1.4067, 1.5067] 2/7 (3/15) 284 (349) 0.685 (0.037)

0.7736 [0.7529, 0.8529] 2/5 (4/9) 228 (353) 0.399 (0.034)

Rocscience Inc. [24] FOS

2.0

FOS

1.6 1.2 Two-grid scheme (fine mesh search) Two-grid scheme (coarse mesh search) Generalized bisection search, α = 0.2 Bisection search, α = 0.5

0.8 0.4 0.0

(a)

0

5

10

15

20

Trial step

1.6

FOS

1.2 0.8 0.4 0.0

(b)

0

5

10

15

20

25

Trial step

1.6

FOS

1.2 0.8 0.4 0.0

(c)

0

4

8

12

16

Trial step

Fig. 9. FOS versus trial step (a) cu2/cu1 = 1.0; (b) cu2/cu1 = 0.6; (c) cu2/cu1 = 0.2.

ally less than 50% for the generalised bisection search with a < 0.5. Fig. 10 demonstrates the slope failure modes for cu2/cu1 = 1.0, cu2/cu1 = 0.6 and cu2/cu1 = 0.2. When cu2/cu1 = 0.6, the slope slip surface is partially dominated by the weak soil layer, and when cu2/cu1 = 0.2, the slope fails along the weak soil layer. To further assess the generalised bisection search algorithms, maxit_nl is increased to 200, and the other parameters are not changed. Because more nonlinear iterations are required for maxit_nl = 200, the Lagrange multiplier approach is applied, and the resultant 2  2 linear system of equations is solved using the modified SSOR (MSSOR) preconditioned symmetric quasi-minimal

Fig. 10. Slope failure mode for (a) cu2/cu1 = 1.0; (b) cu2/cu1 = 0.6; (c) cu2/cu1 = 0.2.

212

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Table 4 Performance of generalised bisection search algorithms with and without two-grid scheme, respectively, for slope stability analysis based on the Lagrange multiplier method (given maxit_nl = 200, tol_nl = 0.001, [a0, b0] = [0.2, 3.0],  = 0.01 and dFOS = 0.1). cu2/cu1

1.0

0.6

0.2

Single grid scheme, a = 0.5 FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.4769 4/11 1088 1.000

1.4089 6/14 1706 1.000

0.6676 5/11 1683 1.000

Single grid scheme, a = 0.2 FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.4768 3/14 970 0.884

1.4081 3/15 1190 0.677

0.6662 4/14 1431 0.803

Two-grid scheme (Rm = 366/3784), a (0.2/0.2) FOS Renewed search range of FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.4727 [1.4202, 1.5202] 1/7 (2/15) 369 (796) 0.331 (0.025)

1.4060 [1.3374, 1.4374] 1/8 (3/15) 524 (1471) 0.303 (0.030)

0.6670 [0.5922, 0.6922] 1/9 (2/15) 1181 (1435) 0.670 (0.035)

Table 5 Performance of generalised bisection search algorithms with and without two-grid scheme, respectively, for 2-D slope stability analysis based on the Lagrange multiplier method (given maxit_nl = 200, tol_nl = 0.001, [a0, b0] = [0.2, 3.0],  = 0.01 and dFOS = 0.2). cu2/cu1

1.0

0.6

0.2

Single grid scheme, a = 0.5 FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.4527 7/14 1725 1.000

1.3662 5/12 1282 1.000

0.5076 7/12 1803 1.000

Single grid scheme, a = 0.2 FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.4499 3/12 714 0.415

1.3649 3/13 838 0.651

0.5073 4/13 1548 0.839

Two-grid scheme (Rm = 168/5629), a (0.2/0.2) FOS Renewed search range of FOS Times of non-convergence/total trials Accumulative nonlinear iterations Normalised total solution time

1.4464 [1.3187, 1.5187] 1/10 (3/16) 336 (1104) 0.193 (0.006)

1.3616 [1.2208, 1.4208] 1/10 (3/11) 358 (985) 0.273 (0.007)

0.5079 [0.4349, 0.6349] 2/6 (2/15) 819 (1717) 0.445 (0.009)

residual (SQMR) solver (e.g. [36]). Table 4 indicates that in the generalised bisection search algorithm, a = 0.2 is still better than a = 0.5 from the perspective of the solution time, which results in times of 857.9 min, 1341.1 min and 782.9 min for cu2/cu1 = 1.0, 0.6 and 0.2, respectively. Compared with the conventional bisection search algorithm, the savings in the solution time achieved by the two-grid scheme are 64.4%, 66.7% and 29.5% for cu2/cu1 = 1.0, cu2/cu1 = 0.6 and cu2/cu1 = 0.2, respectively. Hence, based on the results from Tables 3 and 4 (or Tables 1 and 2), it is postulated that an increase of maxit_nl (to 500 or 1000) results in a more remarkable savings in solution time because of the twogrid scheme, except when cu2/cu1 = 0.2. From Tables 3 and 4, it is observed that the difference in the FOS comes from two sources: different magnitudes in the maxit_nl and different approaches used for the undrained analysis. The adopted approach obviously plays a major role in the difference: when the Lagrange multiplier approach is applied in conjunction with maxit_nl = 50, the computed FOS is approximately 1.467, 1.398 and 0.593 for cu2/cu1 = 1.0, cu2/cu1 = 0.6 and cu2/cu1 = 0.2, respectively, and when the penalty approach is applied, it is 1.508, 1.431 and 0.775, respectively, as shown in Table 3. To be more accurate for the penalty approach, a penalty number larger than Kw/n = 500,000 kPa may be used, but at the expense of increased iterative counts for the iterative linear solvers and the computational cost.

The numerical results also indicate that the Lagrange multiplier approach is evidently superior to the penalty approach, although the former may produce more unknowns. To be more specific, solving the coupled 2  2 linear Eq. (21) using the MSSOR preconditioned SQMR solver leads to a savings of approximately 50% compared with solving the 1  1 linear Eq. (18) using the SSOR preconditioned PCG solver, which confirms the results obtained in [35]. In addition, the computed FOS tends to be larger with an increase in maxit_nl; the FOS is 1.467, 1.398 and 0.593 at maxit_nl = 50 and 1.476, 1.408 and 0.667 at maxit_nl = 200 for cu2/cu1 = 1.0, cu2/cu1 = 0.6 and cu2/cu1 = 0.2, respectively. In the two-grid scheme, the preconditioned iterative solvers (i.e. SSOR preconditioned PCG or MSSOR preconditioned SQMR) are adopted for all linear systems of equations that arise from the 3-D finite element analyses. To demonstrate the effectiveness of the proposed search scheme, a 2-D finite element analysis is conducted for the same slope problem, where all of the parameters associated with the search scheme remaining unchanged, except dFOS (which is defined as dFOS = 0.2 in this example), as shown in Table 5. For the 2-D finite element analyses, the multi-frontal direct solver MA47 from the Harwell subroutine library (HSL) is adopted to solve the resultant symmetric indefinite linear systems of equations because of the Lagrange multiplier method (e.g. [37]). For cu2/cu1 = 1.0,

213

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Table 6 Performance of generalised bisection search algorithms with and without two-grid scheme, respectively, for three 2-D slope examples investigated in [39] (given maxit_nl = 50, tol_nl = 0.001, [a0, b0] = [0.2, 3.0],  = 0.01 and dFOS = 0.1). Analysis method

Lower bound solution [39]

Upper bound solution [39]

Single grid scheme, a = 0.5

Two-grid scheme, a = 0.2

Example 1 FOS Solution time (s)

0.4000 –

0.4500 –

0.4160 283.3

0.4175 206.0 (12.0)

Example 2 (without water table) FOS 1.2500 Solution time (s) –

1.3700 –

1.2838 366.2

1.2857 315.2 (8.6)

Example 3 FOS Solution time (s)

1.2700 –

1.2090 593.6

1.2039 421.2 (12.5)

1.1300 –

Fig. 11. Slope examples investigated in [39] (a) example 1; (b) example 2; (c) example 3.

cu2/cu1 = 0.6 and cu2/cu1 = 0.2, the computational time based on the single grid approach with a = 0.5 is 102.7 min, 76.5 min and 93.2 min, respectively, which are employed here as a benchmark. Using the generalised bisection search scheme, the solution time may be reduced by 16–58%. Compared with the benchmark examples, a more significant savings in solution time, ranging from 55% to 80%, may be achieved using the two-grid search scheme with the mesh ratio set as 168/5629.

5.2.3. Application to various 2-D slope stability examples For geotechnical stability analyses, LEM and limit analysis methods (LAM) have enjoyed great popularity (e.g. [38]). According to Yu et al. [38] and Kim et al. [39], LEM is most likely the most popular method due to its simplicity and wide acceptance in geotechnical practice, whereas LAM may lead to more mechanically rigorous solutions. Due to advances in computer technologies, SSRFEM has received considerable attention in recent years. To further validate the effectiveness and superiority of the proposed solution scheme, various classical slope examples, including a layered soil slope, a slope with a weak layer and a slope with linearly increased

soil cohesion (as shown in Fig. 11), are employed. Because the influences of elastic parameters (that is, the effective Young’s modulus, E0 , and the Poisson’s ratio, m0 ) are not significant, E0 = 10 MPa and m0 = 0.3 are taken in all of the investigated slope examples. Based on the above examples, the SSRFEM approaches are compared with the lower- and upper-bound limit analysis methods. Utilising the prescribed parameters (i.e. maxit_nl = 50, tol_nl = 0.001, [a0, b0] = [0.2, 3.0],  = 0.01 and dFOS = 0.1, as shown in Table 6) for the nonlinear iteration and the search scheme, slope stability analyses are performed using the single-grid and two-grid SSRFEM methods. The numerical results reveal that for the three slope examples, the proposed solution scheme may produce accurate solutions because the computed FOS are all located within the range of [lower-bound solution, upper-bound solution] (e.g. [39]). When the two-grid scheme is applied to these slope examples, the mesh ratio is Rm = 6.6% (296/4487), 3.8% (185/4868) and 4.1% (189/4582); after the coarse mesh search, the search range of FOS is updated to be [0.4178, 0.4278], [1.2018, 1.3018] and [1.1616, 1.2616] for examples 1, 2 and 3, respectively. Table 6 clearly indicates that due to the two-grid search scheme, SSRFEM may lead to a savings of approximately 23.0%, 11.6% and 26.9% in examples 1, 2 and 3, respectively, in the solution time. Moreover, it may be postulated that with the increase in the prescribed maximum number of nonlinear iterations maxit_nl, more remarkable savings in solution time could be achieved. For example 1, because the intermediate soil layer is comparatively weak compared with the other two layers, the plastic zone is primarily located in the weak layer, and yielding may occur at some stress points along the interface between the weak layer and the third soil layer, as shown in Fig. 12a and b, validating that the sliding surface is dominated partially by the weak layer and partially by the interface. From the displacement field and the plastic zone of example 2, as shown in Fig. 12c and d, the maximum slope deformation occurs at the slope crest, and the sliding surface of the slope is composed of the horizontal part in the weak layer and the slant part in the upper layer, implying that the slope may eventually fail and slide along the thin weak layer. In the third slope example, the upper soil layer has the drained soil behaviour, but the two lower layers have undrained soil behaviours. As observed in Fig. 12e and f, the slope fails as if the undrained soils had been squeezed out of the interior of the slope, which may be signified by the larger portion of the plastic zone in the undrained soil layers and the substantial relative movement between the drained layer and the undrained layers. With the stress/displacement field and plastic zone, the picture of slope sliding behaviour can be clearly painted. Although the SSR approach is generally more time-consuming than LEM and LAM, SSRFEM may produce the progressive stress and deformation fields, which could be more crucial in understanding the mechanism of slope failure. More importantly, the SSRFEM approach can be readily extended to

214

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Fig. 12. Displacement field and plastic zone at critical slope failure state for the slope examples investigated in [39] (a) displacement field for example 1; (b) plastic zone for example 1; (c) displacement field for example 2; (d) plastic zone for example 2; (e) displacement field for example 3; (f) plastic zone for example 3.

solve geotechnical problems that involve non-associated soil plasticity, whereas LEM and LAM are generally incapable of conducting such analyses. 6. Conclusions Because the equilibrium state of a geotechnical system with reduced soil strength must be repeatedly evaluated, the SSRFEM method is generally more expensive than the traditional limit equilibrium method, which significantly limits the former’s range of acceptance and its applications. Therefore, the principal aim of this study is to accelerate the SSRFEM computation using several advanced search algorithms for the slope stability analysis. Based on a number of slope examples, several search algorithms used to find the FOS of the slope stability are investigated and compared with each other, and a few concluding remarks are summarised below: (1) The conventional bracketing and bisection search method is inefficient because it takes equal weights on the two search regions, whereas taking unequal weights on the two search regions could be more appropriate from the perspective of efficiency. Hence, a generalised bisection search algorithm is proposed to reduce the possibility of encountering nonconvergence for the nonlinear iterative method by having a bias (that is, a < 0.5) on the convergence region. Based on several 3-D slope examples, it appears that the generalised bisection search algorithm with a small a (such as 0.2) generally encounters fewer instances of non-convergence so that a savings in the solution time is achieved. (2) To further accelerate the search process for the FOS, a new two-grid scheme is proposed, where the coarse mesh search is used to determine the narrowed search range for the FOS

so that [FOS determined by coarse mesh  dFOS, FOS determined by coarse mesh] can be defined for the fine mesh search. Based on the numerical experiments, it is observed that the new two-grid scheme significantly accelerates the entire solution process, particularly when maxit_nl is very large. In several slope examples, compared with the conventional bisection search, 60–70% or even greater savings in solution time is achieved using the generalised bisection search with the two-grid scheme. (3) Two approaches, i.e., the penalty approach and the Lagrange multiplier approach, are investigated in an undrained analysis of a slope example. It is observed that the undrained analyses of geotechnical problems are generally more difficult to solve than the drained analyses, particularly when the penalty approach is applied and the 1  1 block linear system of equations is involved. A numerical comparison also reveals that solving the coupled 2  2 linear Eq. (21) using the MSSOR preconditioned SQMR solver is superior to solving the 1  1 linear Eq. (18) using the SSOR preconditioned PCG solver. First, a savings of approximately 50% in solution time can be achieved. Second, small penalty numbers (Kw/n) may lead to a deviation from the correct solution, whereas large penalty numbers may produce very ill-conditioned linear systems for which iterative linear solvers converge slowly. (4) Based on various classical slope examples, the SSRFEM approaches are evaluated and compared with those obtained using LEM and LAM. The numerical results reveal that although the SSRFEM approaches are very time-consuming, the detailed stress/displacement and yielding fields achieved may lead to a better understanding of the mechanism of slope failure. Hence, greater efforts are needed to improve the efficiency of SSRFEM approaches.

X. Chen et al. / Computers and Geotechnics 62 (2014) 203–215

Conflict of interest The authors declare that there is no conflict of interest. Acknowledgements The research is supported in part by the Research Fund Program of State Key Laboratory of Hydroscience and Engineering, No. 2013-KY-4, the scientific research foundation for the National Basic Research Program of China, No. 2012CB026104, and the National Natural Science Foundation of China, No. 51379103. Appendix A Algorithm 1: Generalised bisection search algorithm for the safety factor of slope stability Given the search range [a0, b0], a0 2 (0, 1/2] and  (such as 0.001 or 0.01), (1) Initialise a = a0, b = b0 and SRF_trial = a; a = a0 (2) For n = 1:10,000 (3) Apply a nonlinear iterative method, such as the Newton– Raphson iteration; nnl records the required nonlinear iterations for each trial (4) If (b  a) <  and (Newton–Raphson may converge) (5) FOS = SRF_trial; Return (6) End (7) SRF_trial = a⁄(1.0  a) + b⁄a; (8) If the nonlinear iterative method does not converge (9) b = SRF_trial; (10) Else if the nonlinear iterative method converges (11) a = SRF_trial; (12) End (13) End

References [1] Sloan SW. Geotechnical stability analysis. Géotechnique 2013;63(7):531–72. [2] Duncan JM. State of the art: limit equilibrium and finite-element analysis of slopes. J Geotech Eng, ASCE 1996;122(7):577–96. [3] Zienkiewicz OC, Humpheson C, Lewis RW. Associated and nonassociated viscoplasticity and plasticity in soil mechanics. Géotechnique 1975;25(4):671–89. [4] Ugai K. A method of calculation of total factor of safety of slopes by elastoplastic FEM. Soils Found 1989;29(2):190–5. [5] Matsui T, San KC. Finite element slope stability analysis by shear strength reduction technique. Soils Found 1992;32(1):59–70. [6] Dawson EM, Roth WH, Drescher A. Slope stability analysis by strength reduction. Géotechnique 1999;49(6):835–40. [7] Griffiths DV, Lane PA. Slope stability analysis by finite elements. Géotechnique 1999;49(3):387–403. [8] Brinkgreve R, Broere W. PLAXIS 2-D. Version 8.5 finite-element code for soil and rock analyses. Complete set of manuals; 2007. [9] Chen X, Phoon KK, Toh KC. Performance of zero-level fill-in preconditioning techniques for iterative solutions with geotechnical applications. Int J Geomech 2011;12(5):596–605. [10] Chen X, Jie Y, Liu J. Robust partitioned block preconditioners for large-scale geotechnical applications with soil–structure interactions. Int J Numer Anal Meth Geomech 2014;38(1):72–91.

215

[11] Cai F, Ugai K. Numerical analysis of the stability of a slope reinforced with piles. Soils Found 2000;40(1):73–84. [12] Zheng H, Liu DF, Li CG. Slope stability analysis based on elasto-plastic finite element method. Int J Numer Meth Eng 2005;64(14):1871–88. [13] Griffiths DV, Marquez RM. Three-dimensional slope stability analysis by elasto-plastic finite elements. Géotechnique 2007;57(6):537–46. [14] Wei WB, Cheng YM, Li L. Three-dimensional slope failure analysis by the strength reduction and limit equilibrium methods. Comput Geotech 2009;36(1):70–80. [15] Huang M, Jia CQ. Strength reduction FEM in stability analysis of soil slopes subjected to transient unsaturated seepage. Comput Geotech 2009;36(1– 2):93–101. [16] Wei WB, Cheng YM, Li L. Three-dimensional slope failure analysis by the strength reduction and limit equilibrium methods. Comput Geotech 2009;36(1–2):70–80. [17] Plaxis Inc. Phi-c reduction and comparison with Bishop’s method. Validation & Verification, Plaxis; 2012. [18] Rocscience Inc. A new era in slope stability analysis: shear strength reduction finite element technique. Toronto, Ontario, Canada; 2004. [19] Kumar P. Shear failure envelope of Hoek-Brown criterion for rockmass. Tunn Undergr Space Technol 1998;13(4):453–8. [20] Hammah RE. The shear strength reduction method for the generalized HoekBrown criterion. In: Alaska rocks 2005, the 40th US symposium on rock mechanics (USRMS); 2005. [21] Shen J, Karakus M. Three-dimensional numerical analysis for rock slope stability using shear strength reduction method. Can Geotech J 2014;51(2):164–72. [22] Smith IM, Griffiths DV. Programming the finite element method. 4th ed. John Wiley & Sons; 2004. [23] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in FORTRAN: The art of scientific computing. 2nd ed. Cambridge: Cambridge University Press; 1992. [24] Rocscience Inc. Application of the finite element method to slope stability. Toronto; 2002. [25] Chang YL, Huang TK. Slope stability analysis using strength reduction technique. J Chin Inst Eng 2005;28(2):231–40. [26] Chen X, Liu C. Search algorithms for safety factor in finite element shear strength reduction method. Chin J Geotech Eng 2010;32(9):1443–7. [27] Won J, You K, Jeong S, Kim S. Coupled effects in stability analysis of pile-slope systems. Comput Geotech 2005;32(4):304–15. [28] Zheng H, Tham LG, Liu D. On two definitions of the factor of safety commonly used in the finite element slope stability analysis. Comput Geotech 2006;33(3):188–95. [29] Hoffman JD. Numerical methods for engineers and scientists. 2nd ed. New York, Basel: Marcel Dekker Inc.; 2001. [30] Briggs WL, Henson VE, McCormick SF. A multigrid tutorial. 2nd ed. Philadelphia, PA: Society for Industrial and Applied Mathematics; 2000. [31] HKS Inc. ABAQUS theory and user’s manual, version 6.9.1; 2009. [32] Chen X, Phoon KK, Toh KC. Partitioned versus global Krylov subspace iterative methods for FE solution of 3-D Biot’s problem. Comput Methods Appl Mech Eng 2007;196(25–28):2737–50. [33] Wang YJ, Yin JH, Lee CF. The influence of a non-associated flow rule on the calculation of the factor of safety of soil slopes. Int J Numer Anal Meth Geomech 2001;25(13):1351–9. [34] Chen X, Phoon KK. Some numerical experiences on convergence criteria for iterative finite element solvers. Comput Geotech 2009;36(8):1272–84. [35] Phoon KK, Chan SH, Toh KC, Lee FH. Fast iterative solution of large undrained soil–structure interaction problems. Int J Numer Anal Meth Geomech 2003;27(3):159–81. [36] Chen X, Toh KC, Phoon KK. A modified SSOR preconditioner for sparse symmetric indefinite linear systems of equations. Int J Numer Meth Eng 2006;65(6):785–807. [37] Duff IS, Reid JK. MA47, a FORTRAN code for direct solution of indefinite sparse symmetric linear systems. Report RAL, 1995: 95-001. [38] Yu HS, Salgado R, Sloan SW, Kim JM. Limit analysis versus limit equilibrium for slope stability. J Geotech Geoenviron Eng 1998;124(1):1–11. [39] Kim J, Salgado R, Lee J. Stability analysis of complex soil slopes using limit analysis. J Geotech Geoenviron Eng 2002;128(7):546–57.