A fast contact detection algorithm for 3-D discrete element method
Erfan G.Nezami,Yousf M.A.Hashash *,Dawei Zhao,Jamshid Ghaboussi
2230Newmark Civil Engineering Laboratory,Department of Civil and Environmental Engineering,University of Illinois at
Urbana-Champaign,205North Mathews Avenue,Urbana,IL 61801,USA
Received 8December 2003;received in revid form 21July 2004;accepted 11August 2004
Abstract
In the discrete element method,determining the contact points between interacting particles and the associated contact normals at each time step is a critically important and time consuming calculation.Common-plane (CP)algorithm is one of the more effec-tive methods for contact detection when dealing with two-dimensional polygonal or three-dimensional polyhedral particles.A new approach,called fast common plane (FCP)method,is propod to find the common plane between polygonal particles.FCP approach recognizes that a common plane has identifying characteristics,which dramatically reduce the arch space for the com-mon plane.In two-dimensions,t
he CP is found by checking only 5possible candidate planes.In three-dimensions,the candidate planes fall within 4types related to the geometry of the particles and their relative positions.Numerical experiments reveal that in three dimensions FCP algorithm can be up to 40times faster than available arch methods for finding the common-plane.Ó2004Elvier Ltd.All rights rerved.
1.Introduction
The idea of using discrete particles in numerical simulations was introduced in the work of Born and Huang [1]and Maradudin [2],where atoms are de-scribed via the concentrated mass and contact force of interacting atoms.Einstadt [3]models bonds be-tween atoms by springs.Cundall [4]introduces the dis-crete element method (DEM)to simulate large deformations in jointed rock formations.Cundall and Strack [5]extend DEM to analyze asmblies of ideal-ized granular particles compod of circular disks and spheres.
During the past two decades,DEM has proved to be a reliable tool to study the behavior of granular materi-als in both micro-and macro-scale.DEM simulations have been ud in large scale geophysical applications such as landslides [6,7]and ice flows [8],as well as many
industrial and mining applications including dragline excavation [9],mixing in tumblers,[10,11],and silo fill-ing [11].
At the micro-level,DEM is ud to investigate distri-bution and evolution of micro parameters (such as con-tact normals and contact forces),and their relation to macro-parameters (such as stress and strain)[12,13].The numerical results are in good agreement with micro-mechanical obrvations from experimental tests on nat-ural sands [14,15],idealized photo-elastic nsitive rods [16,17],and spherical glass [18].
The complexity of the behavior of granular materials and substantial computational time and effort required for DEM computations,however,had limited most of DEM simulations prior to 1990to asmblies of circular disks or spheres.Examples are the DEM codes BALL [19]and TRUBAL [20].Lin and NG [21]provide an extensive list of available circular or spherical DEM codes.In most soil-mechanics related applications,the assumption of spherical particles fails to capture esn-tial aspects of mechanical behavior of the particulate material.
0266-352X/$-e front matter Ó2004Elvier Ltd.All rights rerved.doi:10.pgeo.2004.08.002
*
Corresponding author.Tel.:+12173336986;fax:+12173339464.
E-mail address:hashash@uiuc.edu (Y.M.A.Hashash).
/locate/compgeo
Computers and Geotechnics 31(2004)
沥青铀矿
575–587
As a result,during the last decade,many efforts have been made to incorporate non-spherical particles in DEM implementation.This includes ellips[22,23], connected circular gments[24],super-quadratics[25], polygons[26],ellipsoids[21],polyhedrons[26,27].
2.Contact detection schemes in DEM
Regardless of rapid advances in computer hardware and parallel computation techniques,the huge computa-tional time and effort required to calculate and update contact forces are still a major hindering factor in large scale DEM simulations.For complex particle geome-tries,such as three-dimensional polyhedrons,contact detection subroutines can easily take up to80%of the total analysis time.Applicability of a DEM code is di-rectly related to the efficiency of the employed contact detection scheme.
Contact detection in DEM is usually performed in two independent stages.Thefirst stage,referred to as neighbor arch,is merely a rough arch that aims to provide a list of all possible particles in contact.Among available algorithms for neighbor arching,the most re-cent ones include No Binary Search(NBS)contact detection algorithm[28]and DESS algorithm[29].A re-view of neighbor arch m
ethods is available in[30].
In the cond stage,called geometric resolution,pairs of contacting particles obtained from thefirst stage are examined in more detail tofind the contact points(or contact area if distributed contact forces are considered) and calculate the contact forces.Geometric resolution algorithms strongly depend on complexity of the geo-metric reprentation of particles.For example,if the boundaries of the particles are implicitly reprented by a single function f(x,y,z)=0,then a clod form solu-tion is likely to be available(for example e[5]for contacts between disks and spheres[22],for two-dimen-sional ellips,and[21]for three-dimensional ellipsoids). Efficiency of the contact detection schemes are mostly controlled by the simplicity of the resulting equations.
Where the boundary cannot be reprented by a sin-gle function f(x,y,z)=0,such as in polygons or polyhe-drons,the contact detection can be quite cumbersome. Barbosa[26]introduces a simple algorithm for contact detection between polyhedrons that requires comparing all the vertices of one particle to all faces of the other one and vice versa.The algorithm has a high computa-tional complexity of order O(N2),with N being the num-ber of vertices.Williams and OÕConnor[30]introduce Discrete Function Reprentation algorithm,DFR, which achieves a computational complexity of order O(N).In DFR,the contact between particles is calcu-lated by consi
dering the interaction between the bound-ing boxes of particles.Krishnasamy and Jakiela[31]and later Feng and Owen[32]introduce energy-bad meth-ods forfinding the contact forces,in which a potential energy function is defined for each contact as a function of the overlap area.Cundall[33]introduces the well-known class of‘‘Common-Plane’’(CP)methods:‘‘A common plane is a plane that,in some n,bicts the space between the two contacting particles’’.If the two particles are in contact,then both will interct the CP,and if they are not in contact,then neither inter-cts the CP.As a result of using CP,the expensive par-ticle-to-particle contact detection problem reduces to a much faster plane-to-particle contact detection problem. Once the CP is established between two particles,the normal to the CP defines the direction of the contact normal,which in turn defines the direction of the normal contact force between the two particles.This is espe-cially advantageous for vertex-to-vertex or edge-to-vertex contacts,where the definition of the contact normal is a non-trivial problem.The method has a com-plexity of order O(N)and has been successfully imple-mented in three-dimensional DEM code3DEC[27].
The following ctions introduce a new approach to obtain the CP for polygonal(2-D)and polyhedral(3-D)particles.Determination of contact points and con-tact forces are not discusd.Particles are assumed to be convex,rigid or deformable,while concave particles can be modeled as a combination of veral convex par-ticles attached to each other.
3.Definition of the common plane
The CP is identified by its unit vector normal,n[h1], and any point V0on it as shown in Fig.1.For any point V in the space,the‘‘distance’’d V of that point to any arbitrary plane in the space is defined as
d V¼nÁðV0ÀVÞ;ð1Þwhereby,n is th
e unit vector normal to the plane and V0 is any point on that plane.Both V and V0are described
576 E.G.Nezami et al./Computers and Geotechnics31(2004)575–587
美国队长 彩蛋
in a global Cartesian coordinate system.Eq.(1)divides the space into positive and negative half-spaces,with points in positive half-space have positive distances and points in negative half-space have negative distances to the plane.For any polygonal or polyhedral particle A the‘‘distance’’d A of the particle to any plane in the space is defined as
d A¼
maxðd V
A
Þif d C
A
<0
double date
minðd V
A
Þif d C
A
>0
()
;ð2Þ
where d V
A½h2 is the distance of a vertex V on the particle
to the plane(Eq.(1)),and min{Æ}and max{Æ}reprent minimum and maximum values,respectively,taken
vlaiover all vertices of the particle.d C
A is the distance of
the centroid of the particle to the plane.If a face of the particle is parallel to the plane then more than ver-
tex can define the distance d AÁd C
A ¼O½h3 is of no prac-
tical interest as the CP will never pass through centroid of a particle.Subscripts and superscripts in all equa-tions denote particles and vertices(points),respectively. The vertex(or each of the vertices)that define the dis-tance in(2)is called the‘‘clost vertex’’of that particle to the plane.
For any two particles A and B,a CP is the plane which meets the following three conditions:
Condition1.
Centroids of particles A and B are located on opposite sides of the CP.In this paper it is assumed that the centroid of particle A is located on the negative side and that of particle B in the positive side of CP, Fig.2.
Condition2.
The gap,defined as d BÀd A,is a maximum. Condition3.
d A=Àd B,
d A and d B ar
e the distances o
f particles A and B.
Condition1guarantees that the CP is between the particles.The gap d BÀd A is only a function of direction n of the CP and is independent of the location of the plane in space.Conquently,Condition2identifies the direction n by maximizing the gap.Condition3spec-ifies the location of the CP by tting d A=Àd B.For p-arated particles the gap is always positive(d A<0and d B>0),while for particles in contact the gap is always negative(d A>0and d B<0).
Whenever d BÀd A>TOL,where TOL is a small pos-itive ur-defined tolerance,then the particles are recog-nized as not in contact,no CP is developed.A‘‘potential contact’’is a contact for which0<d BÀd A<TOL the particles are likely to develop new contacts in the next few time steps.A real contact is a contact for which d BÀd A<0.Contact points and forces are found only for real contacts.
4.Conventional algorithm forfinding the CP
Cundall[33]suggests a two-stage procedure forfind-ing the CP:thefirst stage specifies one point on the CP (referred to as the reference point,point M in Fig.3). The cond stage is an iterative process,in which a nor-mal vector n,corresponding to the maximum gap,is found by rotating the CP around the reference point. In two-dimensions,the CP is a line and the rotation is performed around the reference point M.In three-dimensions,two arbitrary orthogonal axes are chon in the CP with their origin at the reference point.The CP is then perturbed around each of them in both neg-ative and positive n.If any perturbation produces a gap larger than that of the current CP,the new CP re-places the current one.In this ca,the clost vertices and the reference point is updated bad on the newly found CP[34].If all the perturbations produce smaller
guardiumE.G.Nezami et al./Computers and Geotechnics31(2004)575–587577
gaps than that of the current CP,the next iteration starts with a smaller perturbation.The iteration process starts from an initial guess(either the CP from the previous time step or the perpendicular bictor of the line that connects the centroids of the particles),and continues until the direction of the CP is found with reasonable accuracy.At any stage of iteration,if the gap exceeds a positive tolerance TOL then the iterative process halts and the contact is deleted.A gap larger than TOL indi-cates that the particles are too far from each other to make a contact.The total number of iterations depends on the accuracy of the initial guess of the CP.In general, the algorithm requires a large number of iteration steps. The number of iteration steps is especially high for the first-time formation of the CP,where the initial guess and the actual CP are very different.
5.Fast identification of common plane candidates
When two particles are not in contact,the definition of the CP can be utilized to limit the number of candi-date common planes and thus significantly reduce the computational cost of common plane lection.
5.1.CP identification in2-D
Statement:In two-dimensions,the CP can be found by checking only5possible candidate planes.
The following provides a stepwi proof of the state-ment above leading to identification of the5possible candidate planes.
Proof.Let A and B be the clost vertices for two not-in-contact particles A and B,respectively.
(i)The CP pass through the midpoint M of gment AB.
Let h measure the angle between the CP and the perpendicular bictor(PB)of the gment AB,as shown in Fig.4.Then
j d A j¼j MA j cos h and j d B j¼j MA j cos h½h4 : The Condition3(Section3)of CP definition, (d A=Àd B),implies that j d A j=j d B j or
j MA j cos h¼j MB j cos h)j MA j¼j MB j:
[CP should pass through the midpoint M of gment AB.
(ii)CP is completely located within the space Sedifice
Space S is the area formed by rays Mm1,Mm2,Mm3 and Mm4,drawn from the midpoint M,parallel to edges AA1,BB1,AA2,and BB2,respectively(the shaded area in Fig.5)(A ray is the part of a straight line beginning at a given point and extending limitlessly in one direction).
Assume that line L,portion of which is located outside space S,is a candidate common plane(Fig.5). Then,vertex B1is clor to this line than vertex B.This implies that AB1and not AB are the clost vertices. This contradicts AB being the clost vertices and the geometric arrangement of the particles.Therefore,line L cannot be a candidate common plane.Similarly all lines located partially or completely outside space S cannot be candidate common planes.
[The common plane is completely located within space S.
(iii)The CP should produce the smallest angle with the PB of the line AB
From Fig.4,d BÀd A=j AB jÆcos h
From Condition2,d BÀd A is maximum
)cos h is maximum
绝望的主妇[angle h is minimum.
(iv)The CP is one offive candidate planes
The CP is the line that
is completely located in space S from proof(ii),and
578 E.G.Nezami et al./Computers and Geotechnics31(2004)575–587
makes the smallest possible angle with the PB of AB from proof (iii).If the PB of the gment AB,is completely located in space S (Fig.6(a)),then:
From proof (iii)the line that makes the smallest possible angle h with PB is the PB itlf.The PB also satisfies proof (ii),
[The common plane is the PB (type a below).
If the PB is not located completely inside the space S (Fig.6(b)),then:
美国艺术留学The PB is not the common plane.
The common plane is the line with the smallest possi-ble angle h to the PB (proof iii).[The common plane is one of the boundary rays Mm 1,Mm 2,Mm 3or Mm 4,(Mm 1in the figure,type b below).Any line inside space S other than on the boundaries will make a larger angle h .
missile
Note that in Fig.6(a)candidate common planes superimpod to rays Mm 1and Mm 2will partly be outside space S and can be eliminated as candidate planes as per proof (ii).However,from an algorithmic/practical numerical implementation point of view it is easier to maintain the as candidate planes and eliminate them using check Conditions (1)&(3)as defined in Section 3
[the CP is one of the following candidates:
Type a :The perpendicular bictor of gment AB .Type b :The lines passing through the mid-point of gment AB and parallel to edges AA 1,AA 2of particle A ,or parallel to edges BB 1and BB 2of particle B .
The number of candidate planes is limited to five.h
5.2.CP identification in 3-D
Statement :In three-dimensions,the candidate planes fall within 4types related to the geometry of the parti-cles and their relative positions.
Proof.Let A and B be the clost vertices for two not-in-contact particles A and B ,respectively.
(i)The CP pass through the midpoint M of gment AB
Similar to the 2-D ca,the CP should pass through the midpoint M .Note that in 3-D,PB and CP are both planes,rather than lines,and the angle h measures the dihedral angle between PB and CP (Fig.7).The dihedral angle is the angle made by two perpendiculars to the interction line of two planes,one in each plane.
[The CP pass through the midpoint M of gment AB
一辈子的英文
.
E.G.Nezami et al./Computers and Geotechnics 31(2004)575–587579
