A semi-spring and semi-edge combined contact model in CDEM and its application to analysis of Jiweishan landslide

Chun Feng, Shihai Li, Xiaoyu Liu, Yanan Zhang

InstituteofMechanics,ChineseAcademyofSciences,Beijing100190,China

A semi-spring and semi-edge combined contact model in CDEM and its application to analysis of Jiweishan landslide

Chun Feng*, Shihai Li, Xiaoyu Liu, Yanan Zhang

InstituteofMechanics,ChineseAcademyofSciences,Beijing100190,China

A R T I C L E I N F O

Articlehistory:

Received 12 August 2013

Received in revised form

25 November 2013

Accepted 9 December 2013

Continuum-based discrete element method (CDEM)

Contact detection method

Semi-spring

Semi-edge

Landslide

Continuum-based discrete element method (CDEM) is an explicit numerical method used for simulation of progressive failure of geological body. To improve the eff i ciency of contact detection and simplify the calculation steps for contact forces, semi-spring and semi-edge are introduced in calculation. Semispring is derived from block vertex, and formed by indenting the block vertex into each face (24 semisprings for a hexahedral element). The formation process of semi-edge is the same as that of semi-spring (24 semi-edges for a hexahedral element). Based on the semi-springs and semi-edges, a new type of combined contact model is presented. According to this model, six contact types could be reduced to two, i.e. the semi-spring target face contact and semi-edge target edge contact. By the combined model, the contact force could be calculated directly (the information of contact type is not necessary), and the failure judgment could be executed in a straightforward way (each semi-spring and semi-edge own their characteristic areas). The algorithm has been successfully programmed in C++ program. Some simple numerical cases are presented to show the validity and accuracy of this model. Finally, the failure mode, sliding distance and critical friction angle of Jiweishan landslide are studied with the combined model.

? 2013 Institute of Rock and Soil Mechanics, Chinese Academy of Sciences. Production and hosting by Elsevier B.V. All rights reserved.

1. Introduction

Continuum-based discrete element method (CDEM) (Li et al., 2004, 2007; Wang et al., 2005) is an explicit approach to simulate the progressive failure of geological mass, which is the combination of fi nite element method (FEM) and discrete element method (DEM). Due to its small deformation assumption, false contact and block embedding will take place when large translation and rotation of blocks occur. To solve the above-mentioned problem, a semi-spring and semi-edge combined contact model is introduced.

Three-dimensional (3D) block contact detection method is the key topic in 3D block DEM. The direct method, common plane method, penetration edge approach and incision body scheme are the 4 typical ways to detect the contact relationship between

blocks.

There are about 6 contact types between blocks, i.e. vertex to vertex, vertex to edge, vertex to face, edge to edge, edge to face and face to face. Direct method judges the contact type by the geometry information. This method is easy to realize in programming with a high accuracy of judgment, but it is time-consuming for searching. However, some improvement in direct method has been done to reduce the cost of time. Beyabanaki et al. (2008) presented a point-face contact algorithm used in 3D discontinuous deformation analysis (DDA). According to the relationship between point and face, the contact type and contact normal vector could be easily obtained. Wu et al. (2005) and Wu (2008) presented a new edge-to-edge contact calculating algorithm for 3D DDA, in which the edge-to-edge contacts are simply transformed to be vertexto-face ones. Keneti et al. (2008) introduced a new algorithm for detection of all contact patterns between any two convex blocks. In Keneti’s method, according to the “main plane” and the number of gathered points in another approaching plane, different algorithms for searching real contact points and the type in global coordinate system are applied.

Common plane method (Cundall, 1988) is based on the knowledge that two contact blocks will be completely separated by a plane. By translating and rotating the plane, the vertex number of each block contacting with common plane is obtained, by which the contact type could be determined (Cundall, 1988; Hart et al., 1988; Jing and Stephansson, 2007). Compared to direct method, the search eff i ciency of common plane method has been signif icantly improved, and it has been the main contact search methodin 3DEC (Itasca, 1987). In 3D DDA, this method is also well used (Liu et al., 2004; Yeung et al., 2007). However, it is diff i cult to get the real common plane (CP) for the two polyhedral blocks, and the existing methods for CP identif i cation always need a great amount of calculation and terminate on a saddle-point sometimes. To improve the searching eff i ciency, Nezami et al. (2004) proposed a fast common plane (FCP) method, and the algorithm can be up to 40 times faster than available search methods for fi nding the CP.

Penetration edge approach (Cheng et al., 2006) is based on the theory that if two blocks contact with each other, there will be a nonempty set between them. This approach establishes the mapping from the property of penetration edge (edge type, number of original vertex) to the contact type. This approach judges the contact type according to local topological feature, the substance of which is to enlarge the region of each block and obtain the nonempty set.

Fig.1.Four steps to calculate the contact force.

Incision body scheme (Wang et al., 2006) combines the penetration edge approach and common plane method together. By cutting blocks by 3 typical planes, 16 mapping relationships are obtained.

The fi nal destination of contact detection is to calculate the contact force. For this purpose, 4 steps should be taken, if any above-mentioned method is adopted (shown in Fig. 1, where V–V means the vertex–vertex contact, V–E presents the vertex–edge contact E–E demotes the edge–edge contact): (1) contact state judgment (if has contact or not); (2) contact type determination; (3) contact area calculation; and (4) normal and tangent contact force solving. Each step costs lots of calculating time, so it is important to fi nd an ef fi cient contact model and simplify the computation steps.

2. Prof i le of CDEM and modif i cation of rotation

CDEM is an explicit time history-analysis FEM and DEM approach on fi nite difference principles, and forward-difference approximation is adopted to calculate the progressive process through a time marching scheme. During calculation, the dynamic relaxation method is used to achieve convergence in a reasonable period of time with small time steps, and the convergence is reached when the unbalance ratio is small enough (＜1 × 10-5). Fig. 2 shows the main process to solve a classic geological problem.

CDEM contains two kinds of elements, blocks and contacts (Ma et al., 2011) (Fig. 3). A discrete block consists of one or more FEM elements, all of which share the same nodes and faces. A contact contains several normal and tangent springs, and each spring owns two nodes which belong to two different blocks. Inside a block, the FEM is used, while for contact, the DEM is adopted.

Total content method based on original coordinate is used to calculate the deformation force (in fi nite element) and contact force (in contact element) in traditional CDEM (see Fig. 2). The method to calculate deformation force can be written as

Fig.2.The process to solve a geological problem.

where{F}eis the node force vector of element, [K]eis the stiffness matrix of element, and{u}erepresents the node displacement vector of element.

With large rotation, the distortion of element will happen. To solve the problem, incremental method should be adopted and strain matrix [B] should be used to calculate deformation force instead of stiffness matrix [K]. Besides, strain matrix [B] should be renewed at each time step.

Fig.3.Blocks and interfaces for 8 nodes hexahedron.

Fig.4.Block rotation at given time based on incremental method.

The main steps to calculate node force by strain matrix with incremental method are written as

where [B]i,{Δε}i,{Δσ}i, wi,Jiare the strain matrix, incremental strain, incremental stress, integral coeff i cient, and Jacobi determinant in Gaussian pointi, respectively; {σn}iand {σo}irepresent the new stress and old stress in Gaussian pointi, respectively; [D],{Δu}e, {Fn}eare the elastic matrix, incremental displacement vector and new node force vector of element, respectively;Nis the total number of Gaussian point.

To test the rotation precision of the new method, the numerical case about a cube rotating around a fi xed point under gravity has been done. The model’s size is 10 m × 10 m, with the bottom left corner fi xed. The rotation of the cube in typical time is shown in Fig. 4, and the comparison with the theoretical solution is shown in Fig. 5. Figs. 4 and 5 show that the numerical result and the analytical solution are almost the same, suggesting that it is a good choice to solve rotation problem with the above-mentioned method.

3. Main idea about the combined contact model

The vertex or edge in traditional contact model presents the vertex of the block (8 vertexes for a cube) or the edge of block (12 edges for a cube), which is shared by some faces, without characteristic area. For calculating contact forces by traditional methods, four steps should be taken (shown in Fig. 1). If the vertex or edge is one part of the contact pairs, it is diff i cult to calculate the contact area because vertex or edge has no characteristic area in traditional method. If the contact type is face to face, the overlapping area should be calculated.

The main idea of semi-spring and semi-edge combined contact model is getting the contact force directly, without calculating the contact type and contact area. For this purpose, some skills should be adopted. For fi nding the target face or target edge easily, the vertexes and edges of the element should indent to each face, as shown in Fig. 6. After that, the semi-springs and semi-edges are created subsequently. For a cube case, there are 24 semi-springs and 24 semi-edges. The indentation distance is expressed as

wheredis the indentation distance;Lrepresents the distance between face center and block vertex; and α is the indentation coeff i cient, which would be 0.1%–1% (1% is adopted in this paper). To fi nd the correct contact, the indentation distancedshould be larger than the searching tolerance.

Because the semi-spring and semi-edge are on the face of each element, so they own their characteristic areas as

Fig.5.Comparison between numerical and analytical results.

whereASS,ASEare the area of semi-spring and semi-edge, respectively;Afacepresents the area of mother face which semi-spring and semi-edge are derived from;Nvmeans the vertex number in mother face; andASS-i,ASS-jmean the semi-spring areas. Since the semi-edge consists of two semi-springs, the area of semi-edge is the sum of the area of two corresponding semi-springs.

Fig.6.Semi-springs and semi-edges in a cubic block.

By this combined method, 6 contact types would be reduced to 2 types, which are semi-spring and target face contact, semi-edgeand target edge contact, as shown in Fig. 7. During contact searching, the main loops are semi-spring and semi-edge, so it does not need to determine the contact type and calculate the contact area. Due to the small indentation when creating the semi-spring and semi-edge, once the block vertex or block edge approaches the target face or target edge, the semi-spring and target face contact or semi-edge and target edge contact will be created immediately.“Semi” here means that the indented springs and edges could not form the complete contacts, only if they fi nd the other part (target face or target block edge).

Fig. 8 shows the necessity of indentation. Without indentation, since the block vertexes of elements A–D are in the same position, there are lots of target faces for one semi-spring, but only one of them is real. For semi-edge, the same problem will happen. By indentation, the target face or target edge will be unique.

Fig.7.Semi-spring and semi-edge contact model.

4. Contact state judgment

There are two main linked lists in the program, semi-spring linked list and semi-edge linked list. For searching contacts, the loops of these linked lists are executed respectively, and the details of each loop are shown in Fig. 9.

4.1.Searchmethod

To search the target face or target edge eff i ciently, static cell method is adopted. For traditional cell approach, the main loop should be cells, so it is very expensive for large simulations where the spatial distribution of objects is sparse and irregular (large number of empty cells). In the combined contact model, the main loop is semi-spring and semi-edge linked lists, thus looping the cells without block can be avoided.

To improve the search eff i ciency, another skill called potential objects array is adopted. The array is set in each semi-spring (semi-edge could also obtain the potential objects according to semi-springs) with fi xed length (8 potential objects). In each time step, the semi-spring and semi-edge seek the target face and target edge from the potential objects array, and the array will be renewed in some steps (for example, 100–1000 steps). The search method is shown in Fig. 10.

Fig.8.Necessity of indentation.

Fig.9.Two linked lists in the combined contact model.

4.2.Geometryalgorithmsinsemi-springcontact

There are two cases that the target face would be found: (a) the distance between semi-spring and target face is smaller than the searching tolerance, and the projective point of semi-spring locates on the target face or along the edge of target face; (b) semi-spring locates inside of the block (Fig. 11).

In case (a), the distance between semi-spring and face should be calculated fi rst:

In case (a), if the distance is smaller than the searching tolerance (0.01%–0.1% of the characteristic length), the relationship between semi-spring and target face should be checked. The followingequation is used to check if the semi-spring locates inside the face (including along the edge of face) or not:

Fig.10.Static cells and potential objects array.

Fig.11.Two cases that target face would be found.

To avoid embedding, case (b) should be considered. If the semispring locates inside a block, the face nearest to the semi-spring should be found, and it is considered as the target face. The method to determine whether the semi-spring has embedded or not is

4.3.Geometryalgorithmsinsemi-edgecontact

To search the target edge, the target face where the target edge is located should be found fi rst. The loop of potential objects array should be taken, and for each object, the most probable face would be found fi rst by

After fi nding the most probable face in each block, each edge of the right face should be judged by

Eq. (11) is used to judge whether the two edges (semi-edge and target edge) are parallel or not. The following equation is used to get the distance between these two edges:

Fig.12.Contact types that semi-spring could solve.

4.4.Contacttypesthemodelcouldsolve

The semi-spring and semi-edge combined contact model could solve all the contact types, as shown in Figs. 12 and 13, although it does not need to check the contact type.

In Fig. 12, for vertex to vertex contact type, since each vertex belongs to three faces, the semi-springs in blocks A and B could fi nd their target faces respectively from the three potential faces in opposite part. For vertex to edge contact type, due to the edge in block B located in two faces, the semi-springs in block A could fi nd the target face from the two potential faces. For vertex to face contact type and face to face (face vertex in face) contact type, the semi-springs in block A could fi nd the target face in block B easily according to the formulas in Section 4.2.

In Fig. 13, for edge to edge contact type, the semi-edge in block A could fi nd the target edge in block B, and simultaneously, the semiedge in block B could fi nd the target edge in block A. For edge to face contact type, as the face in block B contains four edges, the semiedges in block A could fi nd the target edge from the four potential edges. For face to face (without face vertex in face) contact type, each face consists of four edges, thus the semi-edges in blocks A and B could fi nd their target edges respectively from these potential edges in opposite side.

5. Contact force calculation

Fig.13.Contact types that semi-edge could solve.

Fig.14.Real contact area during initial failure.

Penalty function will be used to calculate the contact force. Once the semi-spring and semi-edge contact is established, the normal and tangent numerical springs will be created. According to incremental method, the elastic force of the contact can be calculated by

whereFnandFsare the normal and tangent contact forces, respectively;KnandKsare the normal and tangent stiffnesses, respectively; Δdnand Δdsare the relative normal and tangent incremental displacements of spring, respectively. Once contact force is calculated, it should be added to the global node force array.

To simulate the progressive failure process of geological body, failure model should be adopted, thus the maximum tensile criterion and Mohr–Coulomb criterion are used:

whereTmeans the tensile strength,cis the cohesion, φ means the friction angle andAis the characteristic area gotten by Eqs. (4) and (5).

When Eq. (14) is used to judge the initial failure of material, attention should be paid to the characteristic area. If the grid is continuous, the point on the face is not only semi-spring, but also interpolation point (Fig. 14). There are two parallel springs at the same place, which means the contact force between two blocks is undertaken by 8 springs. Thus for each spring, the real spring force would be half, and the real characteristic area for each semi-spring and semi-edge should be half too.

Fig.16.Failure angle test model.

6. Precision test

6.1.Failurecriteriatest

To test the accuracy of the maximum tensile criterion and Mohr–Coulomb criterion, four types of numerical cases are designed (Fig. 15), and the accuracy of tensile failure, shear failure and toppling failure could be verif i ed. Results (Table 1) show that the numerical and theoretical critical failure values are almost the same, suggesting that the main idea of semi-spring and semi-edge combined contact model is reasonable and correct.

6.2.Semi-springcontactmodeltest

To prove the accuracy of semi-spring contact model, the failure angle test of granite blocks is executed. The experiment has been conducted by Li et al. (2005, 2007) and the model is shown in Fig. 16. The model consists of 74 granite blocks, with the size 10 cm × 10 cm × 20 cm, 10 cm × 10 cm × 10 cm, and 5 cm × 5 cm × 5 cm, respectively. The cohesion and tensile strengths of interface between blocks are 0, and the friction angle is 26°. During the experiment, uplift the platform gradually until the model failure, and then record the failure angle θ. Numerical result of failure angle is 21°, which is consistent with the experimental result (19.5°–21.8°). The failure process of granite blocks is shown in Fig. 17, which also agrees with the experimental result. Besides, when the granite blocks topple, face–face (with face vertex in face), face–vertex, vertex–vertex and vertex–edge will happen. This numerical test also conf i rms that the semi-spring contact model could detect the four types of traditional contact (mentioned in Fig. 12) well.

Fig.15.Numerical cases for verif i cation of the failure criteria.

Table1Statistics of numerical cases.

Fig.17.Failure process of granite blocks.

Fig.18.Movement of block under gravity.

6.3.Semi-edgecontactmodel

To check the reliability and accuracy of semi-edge contact model, a numerical model containing three types of contact (shown in Fig. 18) has been set up. The model contains 4 blocks. Three blue blocks are fi xed, and the purple block falls down under gravity. The movement of top block (Fig. 18) reveals that semi-edge contact model could detect edge-face, face–face (without face vertex in face) and edge-edge contact well. Besides, the numerical sliding acceleration (g= 2.927 m/s2) is almost the same as the analytical solution (g= 2.928 m/s2), which could be calculated by

7.SimulationofJiweishanlandslide

7.1.Background

Jiweishan landslide (Figs. 19 and 20) is a typical bedrock landslide in Wulong, Chongqing, China, which occurred on June 5, 2009, with a total volume about 7 × 106m3. The dimension of the landslide mass is about 780 m (length) × 260 m (width) × 60 m (thickness). According to fi eld investigation and geological analysis, the principal failure mechanism is the virtual dip slipping, with the real dip direction 345°∠21°(Fig. 21). When landslide occurred, the key body crumbled fi rst, and then the sliding body slided along the sliding bed; after mass center moved out of sliding bed, toppling happened.

7.2.Realmodelsimulation

For simulating the sliding process and studying the relationship between friction angle and sliding distance, the real model based on the contour line is set up (Fig. 22), where groups 1–3 are bedrock, key body, and sliding body, respectively.

Fig.19.Photo of Jiweishan landslide.

Fig.20.Sketch of sliding body and key body (unit of the length: m).

Fig.21.Sketch of virtual dip slipping.

By combined contact model mentioned above, the fi nal failure state is achieved (Fig. 23). In Fig. 23, sliding faces mean interfaces between different groups, while other faces mean the interfaces between different elements in sliding mass. In Fig. 23, with the increase of inner friction angle, the sliding distance and fragmentation degree of landslide mass decreases gradually.

7.3.Simplemodelsimulation

Based on simple model, typical failure mode (Fig. 25) is obtained, and the relationship between real dip direction (represented by 2 components) and the critical failure friction angle is also obtained. The results are shown in Figs. 26–28, where DACW means dip anglecomponent in west, DACN means dip angle component in north, and CFFA presents critical failure friction angle.

Fig.22.Real model.

Fig.23.Final failure states in different inner friction angles.

In Figs. 26–28, when DACW keeps constant, the CFFA increases linearly with the increase of DACN, with the slope 0.90. When DACN keeps constant, CFFA decreases linearly with the increase of DACW, with the slope of–0.19. If the CFFA keeps constant, the relationship between DACN and DACW will be linear approximately.

Fig.24.Simple model.

Fig.25.Typical failure mode.

Fig.26.Relationship between DACW and DACN associated with CFFA.

Fig.27.Relationship between DACN and CFFA associated with DACW.

Fig.28.Relationship between DACW and CFFA associated with DACN.

8. Conclusions

Semi-spring and semi-edge combined contact model is used to calculate the contact forces between FEM elements, and it is also suitable for arbitrary convex polyhedrons. According to abovementioned model, six contact types could be reduced to two, which are semi-spring and target face contact and semi-edge and target edge contact, so it will save some computing time to some extent. Some simple numerical cases are presented to show the validity and accuracy of this model.

Inner friction angle is an important mechanical parameter for the slide of Jiweishan. With the change of the friction angle, stability state, failure mode and sliding distance also alter. Besides, there is a good linear relationship between dip angle component (DACN and DACW) and critical failure friction angle (CFFA).

Conf l ict of interest

The authors do not have any possible conf l icts of interest.

Acknowledgments

The authors would like to thank the supports from the National Basic Research Program of the Ministry of Science and Technology of China (Grant No. 2010CB731506), the National Key Technology Research and Development Program of the Ministry of Science and Technology of China (Grant No. 2012BAK10B01) and the Youth Science Fund of National Natural Science Foundation of China (Grant No. 11302230).

Beyabanaki SAR, Mikola RG, Hatami K. Three-dimensional discontinuous deformation analysis (3D DDA) using a new contact resolution algorithm. Computers and Geotechnics 2008;35(3):346–56.

Cheng YM, Chen WS, Zhang YH. New approach to determine three-dimensional contacts in blocks system: penetration edges method. International Journal of Geomechanics 2006;6(5):303–10.

Cundall PA. Formulation of a three-dimensional distinct element model—part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1988;25(3):107–16.

Hart R, Cundall PA, Lemos J. Formulation of a three-dimensional distinct element method model—part II: mechanical calculations for motion and interaction of a system composed many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1988;25(3):117–25.

Itasca Consulting Group Inc. 3DEC-3D distinct element code. Minneapolis, USA: Itasca Consulting Group Inc; 1987.

Jing LR, Stephansson O. Fundamentals of discrete element methods for rock engineering theory and applications. Amsterdam: Elsevier; 2007.

Keneti AR, Jafari A, Wu JH. A new algorithm to identify contact patterns between convex blocks for three-dimensional discontinuous deformation analysis. Computers and Geotechnics 2008;35(5):746–59.

Li SH, Zhao MH, Wang YN, Rao Y. A new numerical method for DEM-block and particle model. International Journal of Rock Mechanics and Mining Sciences 2004;41(3):436.

Li SH, Lian ZZ, Wang JG. Effect of rock mass structure and block size on the slope stability-physical modeling and discrete element simulation. Science in China Series E (Engineering and Materials Science) 2005;48(Supp.):1–17.

Li SH, Wang JG, Liu BS, Dong DP. Analysis of critical excavation depth for a jointed rock slope by face-to-face discrete element method. Rock Mechanics and Rock Engineering 2007;40(4):331–48.

Liu J, Kong XJ, Lin G. Formulations of the three-dimensional discontinuous deformation analysis method. Acta Mechanica Sinica 2004;20(3):270–82.

Ma ZS, Feng C, Liu TP, Li SH. A GPU accelerated continuous-based discrete element method for elastodynamics analysis. Advanced Materials Research 2011;320:329–34.

Nezami EG, Hashash YMA, Zhao DW, Ghaboussi J. A fast contact detection algorithm for 3D discrete element method. Computers and Geotechnics 2004;31(7):575–87.

Wang JQ, Lin G, Liu J. Static and dynamic stability analysis using 3D-DDA with incision body scheme. Earthquake Engineering and Engineering Vibration 2006;5(2):273–83.

Wang YN, Zhao MH, Shihai L, Wang JG. Stochastic structural model of rock and soil aggregates by continuum-based discrete element method. Science in China Series E (Engineering and Materials Science) 2005;48(Supp. 1):95–106.

Wu JH, Juang CH, Lin HM. Vertex-to-face contact searching algorithm for threedimensional frictionless contact problems. International Journal for Numerical Methods in Engineering 2005;63(6):876–97.

Wu JH. New edge-to-edge contact calculating algorithm in three-dimensional discrete numerical analysis. Advances in Engineering Software 2008;39(1):15–24.

Yeung MR, Jiang QH, Sun N. A model of edge-to-edge contact for threedimensional discontinuous deformation analysis. Computers and Geotechnics 2007;34(3):175–86.

*Corresponding author. Tel.: +86 10 82544148.

E-mailaddress:fengchun@imech.ac.cn (C. Feng).

Peer review under responsibility of Institute of Rock and Soil Mechanics, Chinese Academy of Sciences.