Numerical Methods and Implementation in Geotechnical Engineering – Part 2

By Y.M. Cheng, J. H. Wang and Li Liang

Numerical Methods and Implementation in Geotechnical Engineering explains several numerical methods that are used in geotechnical engineering. The second part of this reference set includes more information on the distinct element method, geotechnical optimization analysis and reliability analysis. Information about relevant additional numerical methods is also provided in each chapter with problems where applicable. The authors have also presented different computer programs associated with the materials in this book set which will be useful to students learning how to apply the models explained in the text into practical situations when designing structures in locations with specific soil and rock settings.

This reference book set is a suitable textbook primer for civil engineering students as it provides a basic introduction to different numerical methods (classical and modern) in comprehensive readable volumes.

• Language: English
• Publisher: Bentham Science Publishers
• Release date: Apr 20, 2020
• ISBN: 9789811437427

Book preview

PREFACE

For most of the geotechnical problems, particularly those related to real life problems, analytical solutions are usually not available. For both research and practical applications, numerical methods and computer programs are required for many cases. In the recent forty years, many numerical methods have evolved for various kinds of engineering problems. Engineers are now well adapted to the uses of different computer programs for the solution of engineering problems. There is however a major drawback in the current engineering practice in that most of the engineers are not familiar with the basics of the numerical methods, the methods of implementations and the limitations of the numerical methods/programs. In fact, to a certain extent, the methods of implementations and the limitations of the numerical methods are related. In many internal studies using different commercial numerical programs, the authors sometimes found noticeable or even completely different results with different programs or the same program with different default setting for a given problem, and this situation is not uncommon. For a problem with unknown solution, how an engineer assess the acceptability of the computer results is a difficult issue that needs serious attention. In several technical meetings in the Hong King Institution of Engineers, the authors have discussed with some engineers about the appreciation of the limitations of the daily-used engineering programs. If two computer programs can produce significantly different results, how an engineer determine the acceptability of the results actually require deeper knowledge about the basics of the numerical methods and implementations. Interestingly, the authors like to ask the students a question Different answers can be obtained from different commercial programs. Which results should be accepted, and why should those results be accepted?. In general, the authors challenge the students (undergraduate and graduate students) every year for this question, and virtually this question is never answered properly. The problems in the assessment of the numerical results will also be discussed in this book, which is seldom addressed in other books or research papers.

The authors have participated in different types of geotechnical research and consultancy works in different countries, and has written a book Frontier in Civil Engineering, Vol.1, Stability Analysis of Geotechnical Structures, which is well-favored by many students, engineers and researchers. Most of the books on numerical methods seldom address the actual procedures in numerical implementations, but many postgraduates actually need to develop computer programs to consider special constitutive models, loadings, numerical methods, boundary conditions and other effects. In view of the limitations of most of the books at present, the authors would like to write a new book on numerical methods and the implementations based on their previous works, and this new book should be useful for senior undergraduates, postgraduates, engineers as well as researchers.

In this book, finite element method, optimization method, plasticity based slip line method, limit analysis method, distinct element method, Smoothed-Particle Hydrodynamics Method, Spectral Element Method and Material Point Method will be introduced. The present book will not cover dynamic problems which is a big topic, and hopefully this will be covered later by the authors in another book. The authors will also try to explain the methods of implementation for some of these methods through sample computer programs. Sample programs are given and discussed to assist students in developing programs for their own uses. These programs are not meant to be efficient or up-to-date, but will help the students in learning about the implementation of some numerical methods. This book should not be taken as a classical textbook, as the authors do not intend it to be. There are many new contributions to numerical methods in geotechnical engineering over the last 30 years, and many topics can be covered by individual books for detailed discussion. There is also no way for the authors to cover all numerical methods in details in this book. This book is a basic introduction to some more commonly used numerical methods in geotechnical engineering which have been used by the authors for teaching and research, with the discussion of some common commercial program problems, programming techniques and applications.

ACKNOWLEDGEMENTS

The authors would like to acknowledge the support from the Hong Kong Research Grant Council through the project PolyU 5128/13E, the National Natural Science Foundation of China (Grant No. 51778313), the Cooperative Innovation Center of Engineering Construction and Safety in Shangdong Blue Economic Zone. The present book is also partly supported by CityU Strategic Research Grant for unfunded GRF/ECS (SRG-Fd): Enhancement of Building Information Modelling (BIM) in Construction Safety by GIS: 4D Modelling, Geo-spatial Analysis and Topography Modelling (Preliminary Study) (CityU Project No. 7004631).

CONSENT FOR PUBLICATION

Not applicable.

CONFLICT OF INTEREST

The authors confirm that this chapter contents have no conflict of interest.

Y.M. Cheng

School of Civil Engineering, Qingdao University of Technology,

Qingdao,

China

Department of Civil and Environmental Engineering, Hong Kong,

Polytechnic University,

Hong Kong

Distinct Element Method and Other Numerical Methods

Y.M. Cheng, J.H. Wang, Li Liang, W.H. Fung Ivan

Abstract

In this chapter, the basic formulation of the distinct element method is introduced, which will be followed by more discussions on the blocky and particle formulation of the DEM. This will be followed by some applications of the DEM to several types of problems. As the alternatives to the DEM, DDA and manifold for two-dimensional and three-dimensional problems will be introduced with case studies. There are also many new numerical methods that have evolved in the last twenty years. For example, some recent continuity based methods as such as the SEM, SPH and MPM will also be introduced in this chapter, with applications to some slope failure problems in Hong Kong.

Keywords: DDA, Distinct element, Manifold, Meshless, MPM, Particle flow, SEM, SPH.

1.1. INTRODUCTION TO DISTINCT ELEMENT METHOD

While the finite element is the most popular method adopted by the engineers for practical works, there are some other numerical methods which are under active research by various researchers. These methods are generally not mature enough or difficult to be used for general engineering problems, nevertheless, they can be applied in many geotechnical problems governed by large displacement, discontinuity or even separation of materials. The authors do not aim at introducing all of these special numerical methods, but will discuss some of these methods with applications.

The most commonly used numerical methods for continuous systems are the finite difference method (FDM), the FEM and the boundary element method (BEM). The basic assumption adopted in these numerical methods is that the materials concerned are continuous throughout the physical processes. This assumption of continuity requires that at all points in a problem domain, the material cannot be torn open or broken into pieces. The neighbourhood of all material points will remain unchanged throughout the whole physical process. Some special algorithms have been developed to deal with material fractures in continuum

mechanics based methods, such as the special joint elements by Goodman (1976) and the displacement discontinuity technique in BEM by Crouch and Starfield (1983). However, these methods can only be applied with limitations (Jing et al. 1993):

(1) large-scale slip and opening of fracture elements are prevented in order to maintain the macroscopic material continuity;

(2) the amount of fracture elements must be kept to relatively small so that the global stiffness matrix can be maintained well-posed, without causing severe numerical instabilities; and

(3) complete detachment and rotation of elements or groups of elements as a consequence of deformation are either not allowed or treated with special algorithms.

The authors have adopted the zero thickness interface element by Goodman (1976) in developing a soil-structure interaction program. One of the practical problems in using this joint element is the difficulty in mesh generation. Currently, most of the engineers will specify a thin layer of material to simulate the interface instead using the joint element directly in finite element analysis. Rock masses are always dissected by joints, faults, cracks or other discontinuities which control the failure and sliding of the masses. The complicated structure of rock masses and the discontinuities lead to a complicated non-linear mechanical response. Research work in the deformation in jointed rock mass is very important in modeling complicated important geotechnical problems. For such system, the use of the classical continuity based numerical methods will face great difficulty if the deformation is large. The use of the discontinuity based approach will be more appropriate under such condition.

Before a slope starts to collapse, the factor of safety serves as an important index in both the LEM and SRM to assess the stability of the slope. The movement and growth after failure have launched which are also important and in many cases cannot be simulated on the continuum model. This should be analyzed by the discontinuous mechanics. The commonly used discontinuous numerical methods include: Finite Element Method (FEM) with Joints, Discrete Element Method (DEM), Rigid-Block Spring Method (RBSM), Contact-Spring Model, Discontinuous Deformation Analysis (DDA) and Numerical Manifold Method (NMM). Joint element method with zero thickness was proposed by Goodman, Taylor and Brekke in 1968. It is very effective for problems with only few discontinuities and small deformation. The conditions of no penetration and no tension at the interface however cannot be strictly satisfied. Sometimes, it may lead to great errors between computation results and actual measured results. There was still a major progress in modelling discontinuous problem at that time. The extension of the Goodman joint element to higher order element for improved accuracy has been achieved by many researchers. Desai (1984) also proposed joint element where a small thickness is allowed. Desai’s element is easier for programming purpose, but the formulation of the stiffness matrix needs some special treatments as the thickness of the element is very small. A major limitation of the joint element method is that if there are many discontinuities presented in a rock mass, the use of joint element is practically very difficult to be adopted. The separation of elements also cannot be modeled by the use of joint element, as the method of analysis is still based on the continuum finite element method. Besides the blocky structure of rock, masonry and similar structures, discontinuous medium in form of assembly of grains is also commonly found. There are hence, two major discontinuity based numerical approach: blocky and granular. Actually, the granular approach is a special case of the blocky type problem, but due to the shape of the grains, contact detection can be a much easier task, and special techniques are available for the solution which is not possible for the blocky type problem.

In continuum description of soil material, the well-established macro constitutive models and the parameters can be measured experimentally. On the other hand, a discrete element approach will consider the material as distinct grains or particles that interact with each other. The micro-parameters are difficult to be measured directly. The commonly used distinct element method is an explicit method based on the finite difference principles which is originated in the early 1970s by a landmark work on the progressive movements of rock masses as 2D rigid block assemblages (Cundall, 1971). Later, the works by Cundall are developed to the early versions of the UDEC and 3DEC codes (Cundall, 1980; Cundall and Hart, 1985, 1992). The method has also been developed for simulating the mechanical behavior of granular materials (Cundall and Strack, 1979), with an early code BALL (Cundall, 1978) and Trubal (Strack and Cundall 1984) which later evolved into the codes of the PFC group for 2D and 3D problems of particle systems (Itasca, 1995). Lemos (1983) developed a coupled DEM/boundary element formulation. Later, Lorig (1984) developed the pre and post-processor for DEM analysis, and the code is later modified to HYDEBE (Hybrid discrete element boundary element) which is written in Fortran IV. Through continuous developments and extensive applications over the last three decades, a great body of knowledge and a rich field of literature about the distinct element method have been accumulated. The main trend in the development and application of the method in rock engineering is represented by the history and results of the code groups UDEC/3DEC.

Based on Trubal, Thornton and Randall (1988) further developed the contact model which lead to the GRANULE code. Ng and Dobry (1991) further developed the Trubal code to CONBAL with applications to quartz sand. There are other codes like DIBS by Walton (1982), DISC by Ting et al. (1989), GLUE by Bathurst and Rothenburg (1989) and JP2 by O’Sullivan (2002). Currently, there are many open source particle flow codes (Oval, LIGGGHTS, ESyS, Yade, ppohDEM, Lammps) as well as commercial programs. In general, DEM is still limited to basic research instead of practical applications, as there are many limitations which include: (1) difficult to define and determine the micro-parameters; (2) there are still many drawbacks in the use of matching with the macro response to determine the micro-parameters; (3) not easy to set up a computer model; (4) not easy to include structural element or water pressure; (5) extremely time-consuming to perform an analysis; (6) pre and post-processing are not easy or trivial. It should also be noted that DEM can also be formulated by an energy based implicit integration scheme, however, this approach is not commonly adopted due to the requirement for the solution of a large matrix. The Discontinuous Deformation Analysis (DDA) method is similar in many respects to the force based DEM, and a large matrix will also be obtained from the DDA analysis.

In DEM, the models can either be soft or hard model, which depicts that penetration is either allowed or not during the analysis. This classification should not be mixed up with the rigid or deformable block/particle formulation where the block or particle is rigid or deformable. In general, rigid block/particle formulation is simple in the formulation and analysis, and this approach is suitable for many granular medium or where the stress state is not high. For deep buried tunnels/mines where the insitu stresses are high, the deformation of the blocks or particles can be appreciable and should not be neglected. On the other hand, for the penetration, the authors find that if penetration is completely eliminated in the analysis, the time for analysis can be extremely long, and more importantly, oscillation is commonly encountered and divergence may sometimes be found.

Irrespective of DEM or DDA, block or granular formulation, the use of damping is critical in the calculation. To dissipate the kinetic energy so as to achieve equilibrium, the use of artificial damping is necessary, or else the calculation will keep on fluctuating. There are three major types of damping that can be used: viscous damping, self-adaptive damping and Coulomb damping. For viscous damping, the damping force is proportional to the velocity of the block/particle. For maximum efficiency of computation, the damping coefficient should be chosen to be slightly less than the critical damping value. For low frequency vibration, the use of mass damping can be more effective, while for high frequency vibration, the use of stiffness damping can be more effective. It is also possible to combine these two damping in the actual computation. For a complicated system with different sizes of blocks, it is sometimes difficult to choose a suitable damping before the analysis. In some cases, the introduction of damping introduces body forces and generates incorrect failure modes. The critical damping value usually depends on the eigenvalue of the system, which is however difficult to be assessed for a non-linear system. The use of a single damping values to different sizes of blocks or particles can also introduce some local stability problems. To overcome this problem, Cundall proposed the use self-adaptive damping which is commonly used in many commercial programs. There are two major types of self-adaptive damping, namely viscous type and local self-adaptive type. Self-adaptive damping varies at different locations in a system. It can overcome many of the problems in the viscous damping, but usually the system will be over-damped. Coulomb type is friction type damping for which damping is always pointing opposite to the direction of movement. A detailed treatment about the various computational techniques in DEM has been provided by Jing (2007).

For the time-step in the analysis, it must be controlled so as to achieve convergence in the analysis. A rough guideline is to use eq.(1.1) for the time step

Since ωn=2π/T, hence

Where m and k are the mass and stiffness of the contact spring respectively. In practice, Δt will be chosen as Tmin/10 or similar to ensure convergence in the computation. Towards this, some users may choose to increase the inertia mass, the time-step and efficiency of the computation. This is an approach commonly used for static problem, but this approach cannot be applied for dynamic problems.

In particle DEM, the packing of granular material can be defined from the statistical distributions of the grain size and porosity, and the particles are assigned normal, shear stiffness and friction coefficients in the contact relation. Two types of bonds can be represented either individually or simultaneously; contact and parallel bonds, respectively (Itasca, 1995). Although the individual particles are solid, these particles are only partially connected at the contact points which will change at different time step. Under low normal stresses, the strength of the tangential bonds of most granular materials will be weak and the material may flow like a fluid under very small shear stresses. Therefore, the behaviour of granular material in motion can be studied as a fluid-mechanical phenomenon of particle flow where individual particles may be treated as ‘molecules’ of the flowing granular material. In many particle models for geological materials in practice, the number of particles contained in a typical domain of interest will be very large, similar to the large numbers of molecules.

One of the primary objectives of the particle model is the establishment of the relations between microscopic and macroscopic variables/parameters of the particle systems, mainly through micromechanical constitutive relations at the contacts. Compared with a continuum, particles have an additional degree of freedom of rotation which enables them to transmit couple stresses, besides forces through their translational degrees of freedom. At certain moment, the positions and velocities of the particles can be obtained by translational and rotational movement equations and any special physical phenomenon can be traced back from every single particle interaction. Therefore, it is possible for DEM to analyze large deformation problems and flow process which will occur after slope failure has initiated. The main limitation of DEM is that there are great difficulties in relating the microscopic and macroscopic variables/parameters, hence DEM is mainly tailored towards qualitative instead of quantitative analysis up to the present.

DEM runs according to a time-difference scheme in which calculation includes the repeated application of the law of motion to each particle, a force-displacement law to each contact, and a contact updating scheme. Generally, there are two types of contact in the program which are the ball-wall contact and the ball-ball contact. In each cycle, the set of contacts is updated from the known particles and known wall positions. Force-displacement law is firstly applied on each contact, and new contact force is then calculated according to the relative motion and constitutive relation. Law of motion is then applied to each particle to update the velocity and the direction of travel based on the resultant force, and the moment and contact acting on the particles. Although every particle is assumed as rigid material, the behavior of the contacts is characterized using soft contact approach in which finite normal stiffness is taken to represent the stiffness which exists at the contact. The soft contact approach allows small overlap between the particles which can be easily observed. Stress on particles is then determined from this overlapping through the particle interface.

There have been major developments in DEM over the last twenty years, and many books and research papers have appeared which cover different aspects of the topic. Among the, the works by Konietzky (2004), Jing and Stephansson (2007), Matuttis et al. (2014), Wu (2012), Li et al. (2017), Jebahi et al. (2015), Sarhosis et al. (2016), Norouzi et al. (2016), Wriggers and Panagiotopoulos (1999), Thornton (2015), Chen et al. (2013), Zhao et al. (2011), Mohammadi (2003), Munjiza (2004), O'Sullivan (2011), Zhang and Sanderson (2002), Munjiza et al. (2012), Oñate and Owen (2011a), Oñate and Owen (2011b), Halsey and Mehta (2002) are suggested for further readings to the readers.

1.2. GENERAL FORMULATION OF DEM

DEM runs according to a time-difference scheme in which calculation includes the repeated application of the law of motion to each particle, a force-displacement law to each contact, and a contact updating scheme. The blocks or grains can be deformable or nondeformable in DEM, but deformable formulation is more tedious in the computations. Generally, deformable assumption is used for blocks while nondeformable assumption is used for grains. This is, in general, in compliance with the fact that block formulation which is more commonly used for tunnel or mining problems under high confining stresses, and the deformation of the blocks are generally not negligible. On the other hand, for grains, this formulation is more suitable for surface type problems or problems under low confining stresses. In each cycle, the set of contacts is updated from the known positions. Force-displacement law is first applied on each contact. New contact force is calculated and replaces the old contact force. The force calculations are based on pre-set parameters like normal stiffness, density and friction. Next, law of motion is applied to each particle to update its velocity, direction of travel based on the resultant force, moment and contact acting on the blocks/particles. Force-displacement law is then applied to continue the circulation.

1.2.1. The Force-Displacement Law

The force-displacement law is described for both the ball-ball and ball-wall or for the similar block contacts. The contact arises from contact occurring at a point. For the ball-ball contact, the normal vector is directed along the line between the ball centers. For the ball-wall contact, the normal vector is directed along the line defining the shortest distance between the ball center and the wall. The contact force vector Fi is composed of normal and shear component in single plane surface

The force acting on particle i in contact with particle j at time t is given by

where rj and ri stand for particle i and particle j radii, lij(t) is the vector joining both centres of the particles and kn represent the normal stiffness at the contact. The shear force acting on particle i during a contact with particle j is determined by:

where f is the particle friction coefficient, ks represent the tangent shear stiffness at the contact. The new shear contact force is found by summing the old shear force (min Fij(t)) with the shear elastic force. Δsij stands for the shear contact displacement-increment occurring over a timestep Δt.

where Vijs is the shear component of the relative velocity at contact between particles i and j over the timestep Δt. The choice of timestep can have a major influence on the efficiency and accuracy of the computation. A large timestep can lead to poor accuracy or even divergence, while a small timestep will lead to major increase in the computation time. The displacement and rotation of each block or particle are updated at every time-step according to the following

1.2.2. Law of Motion

The motion of the particle is determined by the resultant force and moment acting on it. The motion induced by resultant force is called translational motion. The motion induced by resulting moment is rotational motion. The equations of motion are written in vector form as follows:

- (Translational motion)

- (Rotational motion)

where xi and θi stand for the translational acceleration and rotational acceleration of the particle or block i. Ir stands for moment of inertia. Fid and Mid stand for the damping force and damping moment. Unlike finite element formulation, there are now three degree of freedom for 2D problem and six degree of freedom for 3D problems. In Cundall and Strack’s explicit integration distinct element approach, solution of the global system of equations is avoided by considering the dynamic equilibrium of the individual particles rather than solving the entire system simultaneously. That means, Newton’s law of motion is applied directly. This approach also avoids the generation and storage of the large global stiffness matrix that will appear in finite element analysis. On the other hand, the implicit DDA approach that will be discussed later will generate a global stiffness matrix which is even larger than that in the finite element analysis. Actually, an implicit DEM formulation for DEM is possible, but this will also result in a large matrix so that many DEM codes do not adopt such formulation.

In a typical DEM simulation, if there is no yield by contact separation or frictional sliding, the particles will vibrate constantly and the equilibrium is difficult to be achieved. To avoid this phenomenon which is physical incorrect, numerical or artificial damping is usually adopted in many DEM codes, and two most common approaches to damping are the mass damping and non-viscous damping. For mass damping, the amount of damping that each particle feels is proportional to its mass, and the proportionality constant depends on the eigenvalues of the stiffness matrix. This damping is usually applied equally to all the nodes. As this form of damping introduces body forces which may not be appropriate in flowing regions, it may influence the mode of failure. Alternatively, Cundall (1987) proposed an alternative method where the damping force at each node is proportional to the magnitude of the out-of-balance-force, with a sign to ensure that the vibrational modes are damped rather than the steady motion. This form of damping has the advantage that only accelerating motion is damped and no erroneous damping forces will arise from steady-state motion. The damping constant is also non-dimensional and the damping is frequency independent. As suggested by Itasca (2004), an advantage of this approach is that it is similar to the hysteresis damping, as the energy loss per cycle is independent of the rate at which the cycle is executed. While damping is one way to overcome the non-physical nature of the contact constitutive models in DEM simulations, it is quite difficult to select an appropriate and physically meaningful value for the damping. For many DEM simulations, particles are moving around each other and the dominant form of energy dissipation is for frictional sliding and contact breakages. The choice of damping may affect the results of computations. Currently, most of the DEM codes allow the use of automatic damping, and manually prescribed damping may be applied under some special circumstances.

To capture the inherent non-linearity behavior of the problem (with generation and removal of contacts, non-linear contact response and stress-strain behavior and others) the displacement and contact forces in a given time-step must be small enough so that in a single time step, the disturbances cannot propagate from a particle further than its nearest neighbours. For most of the DEM programs, this can be achieved automatically, and the default setting is usually good enough for normal cases. It is however sometimes necessary to manually adjust the time-step in some special cases when the input parameters are unreasonably high or low. Most of the DEM codes use the central difference time integration algorithm which is a second order scheme in timestep.

1.2.3. Measuring Logic for Particle Formulation

For particle based DEM, it is found that there are large fluctuations in the results of analysis between particles, which reflect the influence of the local phenomenon. Such results are not surprising, as the results are highly sensitive to the interaction between particles and the timestep under which the results are monitored. Such local results can be meaningless unless the results are monitored over a long time and larger region. A number of quantities in a DEM model are hence defined with respect to a specified measurement circle. These quantities include coordinate number, porosity, sliding fraction, stress and strain rate. The coordination number and stress are defined as the average number of contacts per particle. Only particles with centroids that are contained within the measurement circle are considered in the evaluation of these values. In order to account for the additional area of particles that is being neglected, a corrector factor based on the porosity is applied to the computed value of stress. Using such smoothing techniques, the influence of the local phenomenon can be smoothed for global analysis.

Since measurement circle is used, the stress in particle is described as the two in-plane force acting on each particle per volume of the particle. Average stress is defined as the total stresses in the particles divided by the volume of measurement circle. Thus, the shape of particle is not important, as the reported stress is easily scaled by volume unit. The average stress tensor σij of the material is defined by

Where N is the number of particles contained within V, σijp is the average stress tensor in particle (p), and Vp is the volume of the particles. Homogenization is the transition of the constitutive relations from the microscopic to the macroscopic level, which is required for the proper interpretation of the results and to avoid the fluctuation in the local results. Without such homogenization, the engineers will face vast amount of results which are impossible to be interpreted. Moreover, the engineers will also be frustrated by the major fluctuation in the results so that the results will become meaningless. Homogenization is not just a simple average process. All the basic physical laws should not be violated in the averaging process. Cheng et al. (2009) pointed out the importance of averaging and smoothing of results for the proper interpretation of the wall pressure on silo. In eq.(1.13), a volume has to be defined in the averaging, but there is no universal guideline for this representative elementary volume. For granular problem, this volume may not be a critical issue, and the basic requirement is that it cannot too small to encompass only very few particles. For blocky system, such volume may not exist. The use of averaging is hence not commonly adopted for blocky discontinuous system.

1.2.4. Contact Constitutive Models

The constitutive model acting at a particular contact consists of three parts: a stiffness model, a slip model and a bonding model. The stiffness model provides an elastic relation between the contact force and relative displacement. The slip model enforces a relation between the shear and normal contact forces such that the two contacting balls may slip relative to one another. The bonding model serves to limit the total normal and shear forces that the contact can carry by enforcing the bond-strength limits. For ball to ball contact, the ideal force displacement relation is given by Hertz (1881) which relates the normal stiffness to one-third order of the normal force and two-third order of the shear stiffness of the grain material, while the shear stiffness is a function of the normal stiffness which is controlled by the Poisson ratio. Such relation is not easy to be used in computation, and normal stiffness is commonly assumed to be constant for many problems for simplicity.

To determine the contact forces, the contact detection and contact resolution must be evaluated which are actually computationally intensive. An elegant contact detection algorithm is required to keep track of the neighbor list which can keep on changing with time. Contact resolution requires the accurate calculation of the contact geometry and kinematics, overlap depth/separation, relative tangential motion, overlap area or volume. A contact constitutive model is then used to calculate the contact forces from this information. In most of the DEM codes, the contact is simplified to be a single point. The strains at the contacting particles and the non-uniform stress distributions within the particles are usually not explicitly considered in the DEM simulation. The overlap between the particles is considered to represent completely the effect of the deformation.

The basic DEM formulation provides no resistance to rotation at the contact points while the non-convex and rough surfaces for real soil particles contacts provide resistance to rotation at the contact points. While contact constitutive models have been proposed to account for the rolling resistance at the contact points, the energy dissipation during spin or resistance to spin motion is rarely considered in the current DEM models. The parallel bond model which is available in some commercial programs considers these components of motion by a cement between the particles. There are also various rotation resistance models which have been proposed by various researchers (Iwashita and Oda 1998, Jiang et al. 2005) with success.

1.2.5. Model Generation

For a blocky DEM model, there are three major methods for block system generation: constructive solid geometry (CSG), successive space division (SSD) and boundary representation (BR). In the CSG method, basic topological shapes are generated and combined to form more complex geometry. This approach is commonly adopted for granular materials by forming spheres, ellipses or other shapes to form more complex particles. For the SSD method, a master block which is usually the domains of interest is successively divided, and the program as given in the appendix of this chapter relies on this method of model generation. This method is commonly used for blocky type DEM codes, and the subdivision depends on the definition of discontinuities in form of joints and fractures. In general, a smooth surface division is generated, and the effect of the rough surface is simulated by means of appropriate friction factor and joint model. It should be noted that free blocks can be generated from this scheme, and free blocks may eventually depart from the master block. An example will be provided later for the illustration of this phenomenon. There is a major limitation to the SSD method: incapable to generate non-convex blocks. Luckily, non-convex blocks are not common for practical purposes. The BR method uses closed surfaces and polyhedra to represent the boundaries of blocks. This method requires the user to define vertices, edges, faces and bodies to define fractures, shape, orientation and location as well as the connectivity relations between the bodies This method is more difficult to understand, but is more realistic to model real blocky system with arbitrary shapes, however, the users need to define huge amount of information to form a DEM model, hence this method is not commonly adopted for practical purposes. More details about the mathematical formulation for blocks generation and definitions can be found in Jing and Stephansson (2007).

To generate a DEM granular model which is later subjected to loading, there are four major groups of methods which are shown in Fig. (1.1). In the section expand method, balls were generated and expanded within small box without the top wall, and the process repeats until the total area is filled. In rain method, the balls are dropped into the section until the falling rain has filled the section. Simple expand and rain method appears to be more promising in that it provides a balance between performance and model generation time. It should be noted that for generation of a large model, it can be extremely time-consuming in the computation. The generation of a suitable model is however vital to the study of granular materials. For example, the void ratio, particle size distribution, combinations of several particles to form irregular or nonconvex particles are extremely important towards the understanding of the more realistic behavior of granular materials. More details can be found in Sullivan (2011).

For blocky system, the joint/fracture and the spacing are usually specified for a solution domain. More than 1 set of discontinuity can be defined for a system, and generally, a system will be less stable when the number of discontinuity set increases. As shown in Table 1.1, if there is only one joint set in the problem, the tunnel is safe for different type of shapes. On the other hand, 2 joint sets can lead to failure and falling of blocks easily. The authors have developed simple block generation scheme based on the joint set information as given by the users, which can be found in the appendix of this chapter. More advanced options for block generation are available in some commercial programs. With reference to Fig. (1.2), two set of joints are used to define the blocks for a system by Cheng (2008), and some triangular blocks will be formed due to the formation of the horseshoe tunnel. Due to the excavation, the self-weight of the blocks gradually induce failure of the system, and the progressive failure and falling of blocks are seen in Fig. (1.3). The formation of a realistic blocky system can be very important for the proper analysis of many underground construction works.

Fig. (1.1))

Generation of a DEM particle model.

Fig. (1.2))

Two set Joint 30 &120 (deg) with 1m (x-direction) spacing for a horseshoe tunnel.

Fig. (1.3))

Generation of blocks based on 2 discontinuity sets for a horsehoe tunnel, with the progressive failure and free falling of blocks.

1.2.6. Modelling of Joints/Contacts

A contact surface constructed among two block edges represents a rock joint. The data elements are formed to demonstrate point contacts for each pair of blocks that are in touch. This situation includes blocks touch along their common edge, corner meet a corner or corner contacts with another corner are valued. Blocks may hung-up by the sharp corner, and stress concentration would also be induced by sharp corner. However, it is able to achieve a realistic representation via rounding the corners in order that blocks can slide past one another smoothly during such situation. The blocks/particles can either be deformable or undeformable, and this will affect the model in generating the contact model. If the block/particle is deformable, overlap is allowed, and it is advisable to define the contact in terms of the normal and shear contact relations. If the blocks/particles are undeformable, a directional spring may be sufficiently good to model the contact behavior.

Table 1.1 The DEM analysis result under different number sets of joints and tunnel shape.

There are 2 major types of contacts, namely corner-to-corner contacts and edge-to-corner contacts. The former is more critical because of a rock joint closed along its whole length. It is assumed that the joints extending the joint among the two contacts and separating it into half with each mid-length supporting its own contact stress. In each associated length and point contact, it is required to compute the incremental normal and shear displacements.

For the normal direction, it is assumed that the stress-displacement relation is linear and is controlled by the normal stiffness kn

where

Δσn = increment of effective normal stress

Δun = increment of normal displacement

Considering a joint, if the tensile strength T is exceeded or equivalently less the effective normal stress, effective normal stress is equal to zero. Similarly, a constant shear stiffness (ks) controls the shear strength. The normal and shear stiffness are difficult to be measured directly, which is well-known for everyone working in DEM. The authors adopt the parameters matching method, where the micro-parameters are tuned so that the macro-behaviour can be reflected in the computation. Such technique has been used by the authors and many other researchers. The shear stress (τs) is governed by both cohesive strength (C) and frictional angle (φ).

or, if

Where

= the incremental shear displacement in elastic component

= the total incremental shear displacement.

It is possible that Blocks could be rigid or deformable in the distinct element method. In rigid blocks, the basic formulation demonstrates a set of distinct blocks that do not vary their geometry due to applied loading. It occurs in the material with low deformability and high strength and/or the low-stress environments. However, the individual blocks cannot be assumed to be rigid for many applications. The deformable blocks are partitioned of the finite-difference triangular elements. In a model consisting of a system of continuous and discontinuous joints, the vertices of the triangular elements are treated as grid-points. For each grid-point, eq.(1.18) formulates the equations of motion.

Where

s = the surface enclosing the mass m lumped at the grid-point,

nj = the unit normal to s,

Fi = the resultant of all external forces applied to the grid-point (from block contacts or otherwise), and

gi = the gravitational acceleration.

Contacts between two blocks are usually determined by the Common-Plane or other similar methods. For particles, due to the special nature of the geometry, more efficient contact detection can be achieved. For blocks, the nature of contact can be no contact, vertex to vertex, vertex to edge, vertex to face, edge to edge, edge to face and face to face. To determine the correct type of contact with the minimum time is extremely important in DEM, as contact detection consume a major proportion of the computer time.

1.3. DISTINCT ELEMENT PARTICLE ANALYSIS OF 3D SLOPE

The authors have applied the two-dimensional particle DEM method to model laboratory and field silo and debris flow tests results. In general, the results of analysis are qualitatively acceptable, but usually not quantitatively. There are also limited success in obtaining reasonable results quantitatively, but in general this is not easy. For finite element/difference based slope stability analysis, the assessment of the factor of safety and critical slip surface can be determined with ease. After the initiation of the failure, the displacement will be very large which is difficult to be assessed by the use finite element method. More critically, there will be separation and loss of continuity which is difficult to be assessed by the conventional numerical method, yet the assessment of the spread of failure is important in many practical applications. Towards this issue, the use of distinct element method will be more suitable. There are however very few applications of DEM in slope stability analysis, and this section is devoted towards this application.

1.3.1. DEM Analysis of 3D Slope with Curvature

Rassam and Williams (1999) conducted a survey on the curvature effect on fill slope stability with concave and convex faces by using SRM in FLAC2D. At present, there are only limited works devoted to the effect of curvature on the factor of safety using LEM and SRM. Practically, there is no previous studies about the post-failure conditions and progressive failure mode of 3D slope with curvature with DEM analysis.

An axisymmetric model is generated for the analysis of a homogeneous concave slopes as shown in Fig. (1.4). The radius of curvature of the concave slope is given by R0/H0 equals 0.857, where R0 and H0 are the radius of curvature and height of slope respectively. The slope angle in the numerical model is 45° with c’=0.5kPa and ϕ’=26.5°. Wall to ball friction in the five faces of boundary in PFC3D model is set to 0.5, which coincides with the soil friction as ball to ball friction in particle flow code. The model is simplified into cube in this numerical analysis so as to reduce the tremendous computer time required for analysis, however, in practical condition, the concave slope section is surrounded by continuous soil mass. That is to say, the boundary wall in numerical modeling is used as other soil particles around in reality. So the boundary condition is the crucial factor in a realistic analysis of the problem. Then the in-situ state of stress in the ground before any excavation or construction is generated by setting the initial conditions which is controlled by gravity in a particle flow model. The in-situ stress state is also an important factor because it will influence the subsequent behavior of the model. Initial stresses, however, are not prescribed directly in PFC3D. The initial stress state must be derived from the initial conditions specified for the compaction of the ball assembly, and it is generated from the gravitational loading in particle flow analysis. Based on the procedures as mentioned above, the initial state of a cut concave slope is generated in Fig. (1.4).

The failure process of a concave slope with small cohesion under gravity is analyzed in Fig. (1.5). In the beginning, the soils at the two corner side start to collapse while the soil at the middle of the slope crest remains basically stationary. With increasing settlement at the two corners, the soil in the middle part