Slope stability analysis based on a multigrid method using a nonlinear 3D finite element model

Abstract

The rigid-body limit equilibrium method cannot reflect the actual stress distribution in a rock mass, and the finite-element-based strength reduction method also has some problems with respect to convergence. To address these problems, a multi-grid method was adopted in this study to establish a structural grid for finite element computation and a slip surface grid for computing slope stability safety factors. This method can be used to determine the stability safety factor for any slip surface or slide block through a combination of nonlinear finite element analysis and limit equilibrium analysis. An ideal elastic–plastic incremental analysis method based on the Drucker–Prager yield criterion was adopted in the nonlinear finite element computation. Elasto-plastic computation achieves good convergence for both small load steps and large load steps and can increase computation precision to a certain extent. To increase the scale and accuracy of the computation, TFINE, a finite element parallel computation program, was used to analyze the influence of grid density on the accuracy of the computation results and was then applied to analysis of the stability of the Jinping high slope. A comparison of the results with results obtained using the rigid-body limit equilibrium method showed that the slope stability safety factors determined using finite element analysis were greater than those obtained using the rigid-body limit equilibrium method and were in better agreement with actual values because nonlinear stress adjustment was considered in the calculation.

Keywords

slope / stability / multi-grid method / nonlinear / finite element method

Introduction

The methods most commonly used at present for slope stability analysis are the rigid-body limit equilibrium method and the finite element method (FEM) [1-4]. The former yields a safety factor determined by analyzing the limit equilibrium status of a block. The method is characterized by clear concepts and simple calculations. However, it cannot take nonlinear structural deformation into account (indeed, no deformation is considered), and the method assumes that sliding surfaces reaches an ultimate state of failure simultaneously, which does not reflect the actual stress status of slip surfaces [5]. FEM can be used to determine the stress field and displacement field of the slope but cannot yield a specific value for the slope stability safety factor [6]. Although many researchers have obtained slope stability safety factors using the strength reduction method together with finite element analysis [7-9], this method requires certain failure criteria by which to judge whether a system enters into a limit equilibrium state. Examples of such criteria are limits on the number of iterations or on the ratio of unbalanced nodal forces and external loads. Inevitably, use of such criteria introduces a subjective element into the solution process. The strength reduction method is limited to a certain extent in practice.

Rock and soil materials are typical nonlinear materials, and numerical computation should consider the influence of nonlinearity. This article describes a solution approach based on three-dimensional nonlinear finite element analysis with a multi-grid method adopted to calculate slope safety factors. This approach combines nonlinear finite element analysis and limit equilibrium analysis and fully considers the impact of deformation and rock mass elasto-plastic stress adjustment on slope stability [10,11]. To address the problems of the large computation scale and volume due to fine simulation of the slope, a parallel 3D nonlinear finite element computation program, TFINE, was used to conduct three-dimensional nonlinear finite element analysis of slope stability. Based on an analysis of the influence of grid density on the precision of the results and a comparison with result obtained with the rigid-body limit equilibrium method, the multi-grid method was applied to an analysis of the stability of the Jinping high slope.

Slope stability factor determined using nonlinear finite element analysis

The factor of safety

The rigid-body limit equilibrium method typically used for slope stability analysis is illustrated in Fig. 1. A block is detached by three structural planes (the sliding surface F₁, the sliding surface F₂ and the tension surface F₃). The forces on the block include the force of gravity on the block, W, the osmotic pressures U₁, U₂ and U₃ acting on the sliding surface and the resisting force of the block. When the block is in sliding mode along intersection OM, it is assumed that the tension surface NOL is pulled apart, so no normal opposite force R₃ and tangential force S₃ exist, but the osmotic pressure U₃ is present. Suppose that F₁ and F₂ simultaneously attain their limit states and that the tangential forces S₁ and S₂ parallel the intersection. These assumptions mean that there are only three unknown forces: R₁, R₂ and S₁ + S₂. The three unknown forces can be determined by their equilibrium conditions, and the safety factor K_c can be determined as follows:

(1)K_{\mathrm{c}}={\displaystyle \frac{f{}_{\mathrm{1}}(R_{\mathrm{1}}-U_{\mathrm{1}})+c_{\mathrm{1}}{A}_{\mathrm{1}}+f_{\mathrm{2}}(R_{\mathrm{2}}-U_{\mathrm{2}})+c_{\mathrm{2}}{A}_{\mathrm{2}}}{S_{\mathrm{1}}+S_{\mathrm{2}}}},

where R₁, S₁, U₁, A₁, f₁ and c₁ are the normal force, tangential force, osmotic pressure, area, friction coefficient and cohesion, respectively, of sliding surface F₁, and R₂, S₂, U₂, A₂, f₂ and c₂ are the normal force, tangential force, osmotic pressure, area, friction coefficient and cohesion, respectively, of sliding surface F₂.

When using FEM or the finite difference method to calculate the forces on the block, the action of osmotic pressure is considered, so that U₁ and U₂ in the above equation are not taken into account. In addition, the tangential forces S₁ and S₂ of the sliding surfaces, determined as described above, typically do not parallel the intersection OM, as illustrated in Fig. 2. Therefore, the inclination of the two sliding surfaces is projected onto OM and then the safety factor is determined. The safety factor K determined on the basis of FEM is given by the following equation:

(2)K={\displaystyle \frac{f_{\mathrm{1}}+R_{\mathrm{1}}+c_{\mathrm{1}}{A}_{\mathrm{1}}+f_{\mathrm{2}}{R}_{\mathrm{2}}+c_{\mathrm{2}}{A}_{\mathrm{2}}}{S_{\mathrm{1}}\cos {\alpha }_{\mathrm{1}}+S_{\mathrm{2}}\cos {\alpha }_{\mathrm{2}}}}.

If subzones with different materials of the sliding surfaces are taken into account, the formula for computing the block’s factor of safety should be modified as follows:

(3)K={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^n(f_{\mathrm{1}i}{R}_{\mathrm{1}i}+c_{\mathrm{1}i}{A}_{\mathrm{1}i})}+{\displaystyle \sum _{i=\mathrm{1}}^m(f_{\mathrm{2}i}{R}_{\mathrm{2}i}+c_{\mathrm{2}i}{A}_{\mathrm{2}i})}}{{\displaystyle \sum _{i=\mathrm{1}}^n(S_{\mathrm{1}i}\cos {\alpha }_{\mathrm{1}i})}+{\displaystyle \sum _{i=\mathrm{1}}^m(S_{\mathrm{2}i}\cos {\alpha }_{\mathrm{2}i})}}},

where n is the number of material subzones of sliding surface F₁; m is the number of material subzones of sliding surface F₂; f_{\mathrm{1}i}and c_{\mathrm{1}i} are the friction coefficient and cohesion, respectively, of the material subzones of sliding surface F₁; f_{\mathrm{2}i} and c_{\mathrm{2}i} are the friction coefficient and cohesion, respectively, of the material subzones of sliding surface F₂; A_{\mathrm{1}i} and A_{\mathrm{2}i} are the areas of the material subzones of F₁ and F₂, respectively; R_{\mathrm{1}i} and S_{\mathrm{1}i} are the normal force and tangential force, respectively, on the material subzone area of F₁; R_{\mathrm{2}i} and S_{\mathrm{2}i} are the normal force and tangential force, respectively, on the material subzone area of F₂; and {\alpha }_{\mathrm{1}i} and {\alpha }_{\mathrm{2}i} are the angles of tangential force and intersection of F₁ and F₂, respectively.

Similarly, the formula for computing the safety factor of the sliding surfaces can also be determined. For example, the safety factor for sliding surface F₁ is as follows:

(4)K_{\mathrm{1}}={\displaystyle \frac{f_{\mathrm{1}}{R}_{\mathrm{1}}+c_{\mathrm{1}}{A}_{\mathrm{1}}}{S_{\mathrm{1}}}}.

Multi-grid method

In three-dimensional finite element analysis, stress results are located on Gaussian points. How to map the stress results of a finite element computation onto a slip surface is a key problem in limit analysis based on nonlinear finite element analysis. This paper presents a multi-grid method by which the stresses and forces on a slip surface can be calculated by interpolation. In this method, a slip surface (either a flat surface or a curved face) is divided into a slip surface grid that is independent of the structural grid, as illustrated in Fig. 3, and stresses at the nodes on the slip surface grid are obtained by interpolation of Gaussian points on structural grid. For example, for node A on slip surface P, we search for the nearest Gaussian point B first; then we take all Gaussian points of element M (on the structural grid) to which Gaussian point B belongs for stress interpolation (assuming that each element has 8 Gaussian points):

(5)\begin{array}{ll}{\sigma }_{ij}={\displaystyle \sum _{k=\mathrm{1}}^{\mathrm{8}}{\sigma }_{ij}^{(k)}{S}_k},& S_k=\mathrm{1}/\left(L_k^n{\displaystyle \sum _{k=\mathrm{1}}^{\mathrm{8}}{\displaystyle \frac{\mathrm{1}}{L_k^n}}}\right)\end{array},

where {\sigma }_{ij} is a stress component acting at node A on the slip surface grid, {\sigma }_{ij}^{(k)}is a stress component of the kth Gaussian point of element M; S_k is the weight function of the kth Gaussian point of unit M, and L_k is the distance between A and the kth Gaussian point. Because n = 0, the stress at node A is the arithmetic mean of stresses at Gaussian points; as n\to \infty, the stress at node A approaches the stress at the nearest Gaussian point.

To obtain the resultant force for slip surface P, we divide each quadrilateral element into two triangles, as illustrated in Fig. 4, and determine the element force vector as the sum of the force vectors on the two triangles. The force vector of any triangle in an element is calculated as follows:

(6)F_i=\overline{{\sigma }_{ij}}{A}_j,

where F_i is the force vector on a triangle of an element; {\overline{\sigma }}_{ij} is the mean of all nodal stress vectors in the element, that is,

\overline{{\sigma }_{ij}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}({\sigma }_{ij}^{(A)}+{\sigma }_{ij}^{(B)}+{\sigma }_{ij}^{(C)})

; andA_j is a triangular area vector of the element, A_j=\overrightarrow{BC}\times \overrightarrow{BA}/\mathrm{2}.

We calculate the sum of the normal vectors of each element on the slip surface, as well as the mean, which is the mean normal vector of the slip surface. We then project the resultant vector obtained by summing up the force vectors of the units onto the mean normal direction of the slip surface, and we obtain the normal resultant of the slip surface. Another component is the tangential resultant. In this way, the factor of safety of each slip surface and that of the whole block can be determined.

Nonlinear finite element method

Slope stability analysis is a limit analysis problem. In the category of elastic–plastic analysis methods, limit analysis corresponds to ideal elastic–plastic analysis that employs an associated flow rule. In nonlinear rock mass analysis, the Drucker–Prager yield criterion agrees well with the features of rock and soil materials and parameter values determined from current geologic surveys; therefore, the Drucker–Prager yield criterion is widely used for rock and soil materials, and still appeals to many researchers today [12,13]. Therefore, for the nonlinear computation described in this paper, an ideal elastic-plastic model based on the Drucker–Prager yield criterion was adopted. One of the key problems of elastic–plastic numerical analysis is constitutive relation integration, which directly influences the precision and convergence of the elastic–plastic computation.

The solution approach described in this paper employs ideal elastic–plastic incremental analysis based on the Drucker–Prager criterion, which produces an analytical solution of transitional stress in accordance with the associated flow rule, without forming an elastic–plastic increment matrix [14-16]. The incremental analysis method yields good convergence no matter how large or small the load step size is. Because it uses a fine load step division, the method is a strictly ideal elastic–plastic incremental computation. In the convergence field, the maximum load that may be imposed is close to the limit load and tends to be safe.

The Drucker–Prager yield criterion can be expressed as follows:

(7)f=\alpha I_{\mathrm{1}}+\sqrt{J_{\mathrm{2}}}-k\le \mathrm{0},

where

\{\begin{array}{l}{I}_1={\sigma }_1+{\sigma }_2+{\sigma }_3,\hfill \\ J_2={\displaystyle \frac{1}{6}}\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\right],\hfill \end{array}

where {\sigma }_{\mathrm{1}}, {\sigma }_{\mathrm{2}} and {\sigma }_{\mathrm{3}} are principal stresses and \alpha and k can be determined from the internal friction angle \phi and the cohesive strength c of materials using Mohr–Coulomb fitting criteria.

Elastic–plastic finite element computation is conducted by a series of iterative programs approaching the genuine solution. Consider a typical iterative step, as shown in Fig. 5. Suppose that the initial stress at a certain Gaussian point in the initial iterative step is {\mathbf{\sigma }}_{\mathrm{0}} and that f({\mathbf{\sigma }}_{\mathbf{0}})\le \mathrm{0}. Using the displacement method, the strain increment is \mathrm{\Delta }\mathbf{ϵ}, which corresponds to the elastic stress {\mathit{\sigma }}_{\mathrm{1}}={\sigma }_{ij}^{\mathrm{1}}=\mathit{D}:({\mathbf{ϵ}}_{\mathrm{0}}+\mathrm{\Delta }\mathbf{ϵ}), where \mathit{D} is an elastic tensor; if f({\mathbf{\sigma }}_{\mathrm{1}})>\mathrm{0}, it is necessary to adjust the stress. If the plastic strain increment in this iterative step is \mathrm{\Delta }{\mathbf{ϵ}}^p, then the adjusted stress is \mathit{\sigma }={\mathit{\sigma }}_{ij}=\mathit{D}:({\mathbf{ϵ}}_{\mathrm{0}}+\mathrm{\Delta }\mathbf{ϵ}-\mathrm{\Delta }{\mathbf{ϵ}}^p)={\mathit{\sigma }}_{\mathrm{1}}-\mathit{D}:\mathrm{\Delta }{\mathbf{ϵ}}^p. The approximate associated flow rule in incremental form yields \mathrm{\Delta }{\mathbf{ϵ}}^p=\mathrm{\Delta }\lambda {\displaystyle \frac{\partial f}{\partial \mathbf{\sigma }}}, and we determine the representative value of {\displaystyle \frac{\partial f}{\partial \mathbf{\sigma }}} by {\mathbf{\sigma }}_{\mathrm{1}}. Given the following conditions,

(8)\begin{array}{ll}f(\mathit{\sigma })=\mathrm{0},& \mathit{\sigma }={{\mathit{\sigma }}_{\mathrm{1}}-\mathrm{\Delta }\lambda \mathit{D}:{\displaystyle \frac{\partial f}{\partial \mathit{\sigma }}}|}_{\mathrm{\sigma }={\mathrm{\sigma }}_{\mathrm{1}}}\end{array}.

The adjusted stress can be determined as follows:

(9)\mathbf{\sigma }={\sigma }_{ij}=(\mathrm{1}-q){\sigma }_{ij}^{\mathrm{1}}+p{\delta }_{ij},

where,

q={\displaystyle \frac{w\mu }{\sqrt{J_{\mathrm{2}}}}},p=-mw+\mathrm{3}n{I}_{\mathrm{1}},m=\alpha (\mathrm{3}\lambda +\mathrm{2}\mu ),w={\displaystyle \frac{f}{\mathrm{3}\alpha m+\mu }},

in which J_{\mathrm{2}}, I_{\mathrm{1}} and fare all determined from {\mathbf{\sigma }}_{\mathrm{1}}, and \lambda and \mu are Lame constants.

Therefore, for each iterative step, the stress incremental is as follows:

(10)\mathrm{\Delta }\sigma ={\sigma }_{\mathrm{1}}-\sigma =q{\sigma }_{ij}^{\mathrm{1}}-p{\delta }_{ij}.

In the process of stress adjustment, the above analytical solution for the transitional stress based on the Drucker-Prager criterion is equivalent to the linear prediction–radial return method [17]; as with the constitutive relation integration strategy, it is equivalent to the closest-point projection method (CPPM) [18]. The closest-point projection method has first-order precision and is unconditionally stable. As a special instance of the generalized midpoint rule (GMR), the CPPM can be adapted to a large strain increment. In elastic–plastic finite element computation in geotechnical engineering, load steps are typically very large. However, the above-stated method is suitable for limit analysis, and it can maintain a very high level of numerical stability and precision if a large load step is adopted. Actually, with a large load step, if the above-described method is adopted, the largest load in the convergence field is smaller than the real ultimate bearing capacity of the structure, and the corresponding stress field is the static allowable stress field. At the same time because the associated flow rule can be satisfied under the meaning of average, the largest load in the convergence field is close to the real ultimate bearing capacity of the structure, and the result obtained approaches the real solution and tends to be safer.

Numerical results, comparison and discussion

On the basis of the slope stability analysis method described above, a parallel computation program, TFINE, was compiled using C/C⁺⁺ language and MPICH, implemented in a clustered computer system [19] and successfully applied to slope stability computation and analysis for the Jinping Hydropower Station.

Calculation model

The left bank slope of the Jinping Hydropower Station is 600 m high, as illustrated in Fig. 6. At the dam site, faults, fault crush movements, joint fissures and structural planes are well developed. In the left bank slope in particular, there are deep fissures that are parallel with the bank slope and unfavorable to slope stability.

The finite element computation model for the left bank slope of Jinping is illustrated in Fig. 7. The model consists of 139536 elements and 149585 nodes. The Drucker–Prager yield criterion is adopted in the computation, which simulates the entire process of excavation. From excavation of the cable conveyer platform to the completion of the spandrel groove, a total of 17 steps are involved in the excavation. A quasi-static method is used, in consideration of a transverse seismic load with a peak acceleration of 0.1 g.

There are two key sliding blocks for stability of the left bank slope, as illustrated in Fig. 8:

● Sliding block 1 NNW fractures sets+ X + f42-9

● Sliding block 2 NNW fractures sets+ SL48-3+ f42-9

For the various structures involved, the material parameter values adopted in the computation are given in Table 1.

Interpolation precision of stress for the multi-grid method

In theory, under the self-weight of the stress field, the composite force of the upstream slip surface, the downstream slip surface and the rear pull-apart surface should be balanced with the self-weight of the block in the vertical direction. This can be used to verify the precision of the multi-grid method. The computation for the typical sliding block on the Jinping left bank shows that, letting n in Eq. (5) be 2, the direction of the composite force of the three planes is approximately 2 degrees from the vertical, and the composite force and its self-weight differ by approximately 2%, which means that the computation precision does meet requirements. Therefore, we take n = 2 in the subsequent analysis.

Effect of FEM mesh

To compute the composite force of the block under the effect of self-weight, we set the composite force in the X (across the river) and Y (along the river) directions to be T_xy and that in the Z (vertical) direction to be T_z. Because the only force acting on the block is gravity, we have T_xy/ T_z = 0 for ideal conditions. However, due to various errors in actual computation, it is impossible for the ratio to be zero, although the nonzero value can be used to assess the precision of the computation. The computation results for different grid densities were compared. For the initial grid (with 17442 elements), the result of the numerical computation is T_xy/ T_z = 5%; for the refined grid (with 139536 elements), the result of the numerical computation is T_xy/ T_z = 1%.

It can be observed that refining the grid has a significant influence on the precision of the numerical computation. It can also have an influence on the distribution of normal stress and shear stress on the sliding plane as well as safety factor of the sliding block. Therefore, using a refined grid and large-scale numerical computation technology can effectively improve the precision of the algorithm described in this paper.

The safety factor of the slip surface

In three-dimensional finite element computation, the stress distribution and thus the forces acting on the slip surface can be obtained using the multi-grid method. The safety factor of the sliding block can be obtained from Eq. (6), as illustrated in Table 3.

Based on the results obtained, the following conclusions can be drawn:

1) The safety factor of sliding block 1 is 1.35 for unexcavated conditions, 1.39 for excavation to the dam crest elevation and 1.36 for full excavation, which indicates that excavation to the dam crest elevation can significantly reduce the gravity load of rock and that full excavation may reveal the effect of toe excavation to some extent, as well as significant load reduction.

2) The safety factor of sliding block 2 is 1.434 for unexcavated conditions, 1.67 for excavation to the top of the dam and 1.74 for full excavation, which indicates that the sliding block is significantly influenced by excavation, due to its position adjacent to the outside, and that the effect of load reduction is even more obvious in the excavation process.

3) With a horizontal earthquake load, the normal force on the sliding surfaces decreases, slip resistance decreases, and tangential force increases. As the slip forces increase, the block’s degree of security decreases to some extent. Take sliding block 1 as an example. Its degree of security decreases from 1.35 to 1.199; after full excavation, it decreases from 1.36 to 1.22. For sliding block 2, its degree of security decreases from 1.43 to 1.29 before excavation; it decreases from 1.74 to 1.58 after full excavation.

Comparison with other methods

The major differences between the algorithm described in this paper and the rigid-body limit equilibrium method are as follows:

1) Differences in computation approach. The algorithm described in this paper uses nonlinear finite elements in the stress calculation and considers the influences of nonlinear rock behavior, which makes it agree well with actual stress states. The algorithm also considers the influence of nonlinear deformation on stability, and the slope stability safety factor determined from the corresponding computation also agrees well with actual conditions. The rigid-body limit equilibrium method is used to calculate a limit failure state, and the counterforce of the sliding surface is obtained by equilibrium calculation of a slip block, which considers a limit sliding state and is relatively conservative.

2) Differences in consideration of seepage load. In the algorithm, the seepage load is considered in the stress computation, and seepage pressure is left uncalculated in determining the degree of safety for the purpose of stability. In the limit equilibrium method, the seepage load is incorporated by action on the pull-apart surface.

Using the rigid-body limit equilibrium method, the safety factor of sliding block 1 is 0.975 (for the unexcavated case) or 0.938 (for natural conditions plus earthquake loading), and the safety factor of sliding block 2 is 0.911 (for the unexcavated case) or 0.876 (for natural conditions plus earthquake loading). All these results are smaller than those obtained using the method described in this paper because the rigid-body limit equilibrium method does not consider the influence of rock’s nonlinearity in the calculation process, while the method described in this paper does take nonlinear stress adjustment into consideration in the calculation process and considers that better-quality rock can bear a larger load. In this respect, the method proposed in this paper agrees better with actual conditions.

Determination of slope stability based on finite element analysis mainly involves judging whether a slope is unstable by calculating convergence or a yield region or determining whether plastic strain occurs. Since the calculation of convergence using finite element analysis is influenced by many factors, such as the calculation grid, it is not sufficient to judge destabilization only by convergence. Yield region and plastic strain are also related to the grid status to some extent. Solutions obtained with the method presented in this paper are even more objective, and the method fully considers the effect of stability on nonlinear deformation.

Rock can be both continuous and discontinuous, and this can be reflected in the values assigned to rock parameters: there are continuous body parameters, such as the deformation modulus, Poisson’s ratio and friction angle, including the weakening effect of short joints, and discontinuous parameters, such as the connectivity and shear stiffness of long and large joint sets, which are hard to estimate. Finite element analysis is not well suited to simulating the influence of long and large joint sets in rock, and has difficulty handling cross-elements’ connectivity effect. Therefore, it is necessary to check the results of finite element stress computation for stability against sliding of several key sliding blocks and sliding surfaces when analyzing deformation stability using the finite element method. The multi-grid method takes advantage of the results of three-dimensional nonlinear finite element computation and considers joint connectivity and similar factors in calculating the stability safety factors of sliding surfaces. In this way, the multi-grid method can simulate discontinuous effects in rock. This is in accordance with the principle of multi-dimensional design and analysis.

Conclusions

The multi-grid method was used to establish two grids, a structural grid for finite element computation and a sliding surface grid for calculating a sliding surface’s stability safety factor. This combination of grids makes it easy to determine the stability safety factor of any sliding surface or sliding block, and it also considers the influences of nonlinear deformation and elastic–plastic stress adjustment on the stability safety factor. Using an FE parallel calculation procedure makes it possible to enlarge the finite element computation scale to obtain accurate normal forces, shear force distributions and changes on sliding surfaces, which improves computation accuracy for sliding blocks and provides more practical insights into reinforcement measures required. The results of this study demonstrate that the deformation and stabilization mechanisms of high slopes can be better understood, providing a sound basis for implementing alarms for unsafe conditions and preventing high slope accidents.

References

[1] Bishop A W. The use of the slip circle in the stability analysis of slopes. Geotechnique, 1955, 5(1): 7-17

[2] Duncan J M, 0. State of the art: limit equilibrium and finite-element analysis of slopes. Journal of Geotechnical Engineering, 1996, 122(7): 577-596

[3] Zuyu C, Xiaogang W, Jian Y, Zhixin J, Yujie W. Rock slope stability analysis: Theory methods and programs. Beijing: China Water Power Press, 2005 (in Chinese)

[4] Griffiths D V, Lane P A. Slope stability analysis by finite elements. Geotechnique, 1999, 49(3): 387-403

[5] Lenchman J B, Griffiths D V. Analysis of the progression of failure of the earth slopes by finite elements. In: Proceedings of Sessions of Geo-Denver 2000-Slope Stability 2000. GSP 101, 289: 250-265

[6] Liu Y R. Yang Q, Zhu L. Abutment stability analysis of arch dam based on 3D nonlinear finite element method. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(1): 3222-3228 (in Chinese)

[7] Jiang G L, Magnan J P. Stability analysis of embankments: Comparison of limit analysis with methods of slices. Geotechnique, 1997, 47(4): 857-872

[8] Dawson E M, Roth W H, Drescher A. Slope stability analysis by strength reduction. Geotechnique, 1999, 49(6): 835-840

[9] Zhao S Y, Zheng Y R, Shi W M, Wang J L. Analysis on safety factor of slope by strength reduction FEM. Chinese Journal of Geotechnical Engineering, 2002, 24(3): 343-346 (in Chinese)

[10] Yang Q, Zhu L, Xue L J. Application of limit equilibrium method to stability analysis of Jinping high slope based on 3D multi-grid method. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(Supp.2): 5313-5318 (in Chinese)

[11] Liu Y R, Yang Q, Xue L J, Zhou W Y. Rock Slop Stability Analysis with Nonlinear Finite Element Method. In: Cai M F, Wang J, eds. Boundaries of Rock Mechanics. London: Taylor & Francis Group, 2008, 503-507

[12] Hjiaj M, Fortin J, de Saxce G. A complete stress update algorithm for the non-associated Drucker-Prager model including treatment of the apex. International Journal of Engineering Science, 2003, 41(10): 1109-1143

[13] Yang Q, Yang X J, Chen X. On integration algorithms for perfect plasticity based on Drucker-Prager criterion. Engineering Mechanics, 2005, 22(4): 15-19, 47 (in Chinese)

[14] Yang Q, Chen X, Zhou W Y. A practical 3D elastic-plastic incremental method in FEM based on D-P yield criteria. Chinese Journal of Geotechnical Engineering, 2002, 24(1): 16-20 (in Chinese)

[15] Chen X. Yang Q, Huang Y S. Sub-incremental method for perfect elasto-plastic material based on D-P yield criteria. Chinese Journal of Rock Mechanics and Engineering, 2002, Suppl 2: 2465-2469 (in Chinese)

[16] Yang Q. Chen X, Zhou W Y. Elastic-plastic basis of geotechnical engineering reinforcement analysis. Rock and Soil Mechanics, 2005, 26(4): 553-557 (in Chinese)

[17] Schreyer H L, Kulak R F, Kramer J M. Accurate numerical solutions for elastoplastic models. Journals of Pressure Vessel Technology, 1979, 101(2): 226-234

[18] Ortiz M, Popov E P. Accuracy and stability of integration algorithms for elastoplastic constitutive relations. International Journal for Numerical Methods in Engineering, 1985, 21(9): 1561-1576

[19] Liu Y R, Zhou W Y, Yang Q. A distributed memory parallel element-by-element scheme based on Jacobi-conditioned conjugate gradient for 3-D finite element analysis. Finite Elements in Analysis and Design, 2007, 43(6-7): 494-503

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg