A general method for modelling deformable structures in DEM

744 Molecular Dynamics (MD) or Discrete (Distinct) Element Method (DEM)A general method for modelling deformable structures in DEM IV International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2015 E. O˜nate, M. Bischoff, D.R.J. Owen, P. Wriggers & T. Zohdi (Eds) A GENERAL METHOD FOR MODELLING DEFORMABLE STRUCTURES IN DEM A. EFFEINDZOUROU1, K. THOENI1, B. CHAREYRE2AND A. GIACOMINI1 1Centre for Geotechnical and Materials Modelling The University of Newcastle Callaghan, NSW 2308, Australia http://www.newcastle.edu.au 2Laboratoire 3SR Domaine Universitaire BP53, 38041 Grenoble Cedex 9, France http://www.3sr-grenoble.fr Key words: DEM, Cylinder, Grid, Membrane, Arbitrary shape, Minkowsky sum Abstract. The discrete element method was originally designed to handle rigid spherical particles. In this work, the method is extended to model arbitrary deformable structures, such as membrane and container structures. The paper focuses on the implementation of a deformable discrete element called PFacet (particle facet) in the open-source framework YADE. A PFacet element is geometrically constructed by the Minkowsky sum of a triangular facet and a sphere. Contact detection is based on the three primitives: sphere, cylinder and facet. It allows the representation of arbitrary deformable structures by using the standard contact models for spheres. Finally, the model is validated and an example is presented illustrating the capabilities of the model to simulate arbitrary deformable structures. 1 INTRODUCTION The discrete element method (DEM) is broadly used to simulate granular materials. The method originally considers the interaction between rigid spheres. The modelling of rounded non-spherical particles allows to take into account the shape of the particles. Indeed, grains in nature and many important phenomena cannot be reproduced just using spherical particles. The simplest method for modelling non-spherical particles is to either bond or clump spheres together [1]. In the first case spheres belonging to the same aggregate interact via a bond where attractive forces are allowed and bonds can brake. In the case of a clump, the aggregate is rigid and the forces on all the spheres making up the clump are added together to obtain the total force. Another way of modelling rounded non-spherical particles is by superquadratic ellipsoids [2] or rounded polyhedra 1 745 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini and spheropolytopes [3]. The latter are discussed more in detail in the following section since a similar concept will be used in the present work. Spheropolytopes are generated from the Minkowski addition of a polytope by a sphere, which is nothing more than the object resulting from sweeping a sphere around a polytope.The idea appeared in [4] for flexible objects and in [5] for convex non-spherical particles. In the latest, a linear triangulation-based contact detection method based on a weighted Delaunay triangulation was used in order to improve the performance. Numerical experiments indicated that the overall complexity of the contact detection is proportional to the number of particles. An extension of the method was published in [6] for simulating non-convex non-spherical particles in 2D. The particle shape was represented by the classical concept of a Minkowski sum, which permits the representation of complex shapes without defining the object as a sum of spherical or convex particles. A well-defined conservative and frictional interaction between these bodies was derived. The method was then extended to 3D by [3]. Starting from the modelling of deformable chained cylinders in 3D [7], this paper presents a general method for modelling deformable non-convex structures with the DEM as published in [8] and implemented into the open-source framework YADE [9]. Firstly, the three primitives used in the model and their implementation are discussed. Then the model is validated and some examples are presented. 2 GENERAL FORMULATION The contact model of two interacting spheres is presented herein since all possible interactions considered in this work are based on it. The constitutive model used for the sphere–sphere interaction links the relative displacement and relative rotation to the contact force #» Fand contact moment # » Macting at the contact between two interacting bodies. The normal contact force #» Fnand the incremental shear force d #» Fsare calculated as: #» Fn=kn #» un(1) d #» Fs=ks˙ #» usdt(2) where knand ksare the contact stiffnesses associated to the normal and shear force, #» unis the normal distance (or overlap) between the two spheres, ˙ #» usis the relative shear velocity, and dtis the time step. Considering two spheres S1and S2, the normal and tangential stiffness of the sphere– sphere contact are defined as follows: kn=2E1R1E2R2 E1R1+E2R2 (3) ks=2E1R1ν1E2R2ν2 E1R1ν1+E2R2ν2 (4) where Eiis the Young’s modulus, Rithe radius and νithe Poisson’s ratio associated to sphere iwith i∈{1,2}. 2 746 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini The relative rotation of the spheres is used to define the twisting moment # » Mtand the bending moment # » Mb. By adopting a vector representation #» Ω12 of the relative rotation, the moments can be calculated as: # » Mt=kt #» Ωt 12 (5) # » Mb=kb #» Ωb 12 (6) where ktand kbare the contact stiffness associated to the twisting and bending moment respectively, and #» Ωt 12 and #» Ωb 12 are the twisting and bending components of the relative rotation associated to the two spheres. Plastic deformation is taken into account via the following elastic limits: � #» Fn�≤σl nA � #» Fs�≤Fntan ϕ+σl sA Mt≤σl nIb R Mb≤σl sIt R (7) where σl nand σl sare the tensile and shear strengths, ϕ= min(ϕ1,ϕ 2) is the friction angle, A=πR2is the reference surface area, It=πR4/4 and Ib=πR4/8 are the reference polar and bending moments of inertia respectively, and R= min(R1,R 2) is the reference radius of the contact. 3 CYLINDER ELEMENTS The cylinder element was introduced in [7]. It is the basic entity for the deformable element described in this paper and, hence, it is summarised in this section. A cylinder is composed of two nodes and a tube. Geometrically it corresponds to the Minkowski sum of a sphere and a polyline. The deformation of a cylinder is characterised by the positions and orientations of its two nodes. The mass of the cylinder element is lumped into its two nodes. Three types of interactions are considered: sphere–sphere, sphere–cylinder and cylinder– cylinder. The sphere–sphere contact model is used to calculate the sphere–cylinder and cylinder–cylinder interactions. In order to compute the sphere–cylinder interaction, a virtual sphere is introduced. The virtual sphere is centred at the position of the projection of the contact point between the sphere and the cylinder on the segment connecting its nodes. The cylinder and virtual sphere have the same radius. The displacements at the surface of the cylinder are assumed to vary linearly between its nodes. The translational and rotational velocities of the virtual sphere are linearly interpolated between the cylinder’s nodes. Then, the forces and moments are distributed on the two nodes. In [7] chained cylinders (or interconnected cylinders) were used. This allows the modelling of elastic perfectly plastic beams that can withstand normal, shearing, bending 3 747 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini and twisting loadings. For a cylinder element with two nodes, N1and N2, the contact stiffnesses are defined as: kn=EtA L kt=GtIt L ks=12EbIb L3 kb=EbIb L (8) where Etand Ebare the tensile and bending moduli of the material respectively, Gtis the shear modulus associated with the twisting moment and L=�N1−N2�is the length of the beam element. 4 GRIDS A grid can be generated using a series of nodes connected to an arbitrary number of cylinders. The modelling of a grid involves the implementation of the interaction between non-adjacent cylinders (no common nodes). In order to handle cylinder–cylinder interactions, the contact forces between two cylinders C1and C2are determined. The distance between two cylinders C1and C2is calculated as the distance between two segments (Fig. 1), bearing in mind that each cylinder is geometrically represented by the Minkowski sum of a sphere and a segment. Fig. 1. Distance between two cylinders. The segments associated to the cylinder C1and C2are A=[Nc1 0Nc1 1] and B=[Nc2 0Nc2 1] respectively, where Nci jcorrespond to the position of the nodes jof cylinder Ciwith 4 748 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini i∈{1,2}and j∈{1,2}. The parametric equations of the segments Aand Bare: A(s)=Nc1 0+s#» a∀s∈[0,1] (9) B(t)=Nc2 0+t #» b∀t∈[0,1] (10) with #» aand #» bbeing the direction vectors of the segments Aand B. The vector wbetween Aand Bcan then be calculated as: w(s, t)=A(s)−B(t)∀s, t ∈[0,1] ×[0,1] (11) The objective is to find (sC,tC)∈[0,1] ×[0,1] such that: #» wC=w(sC,t C)=A(sC)−B(tC) = min (s,t)∈[0,1]×[0,1] w(s, t) (12) The segment [ACBC] is uniquely defined as simultaneously perpendicular to both lines if the segments Aand Bdo not intersect or are not parallel. This implies that the vector #» wC=w(sC,t C) is uniquely perpendicular to the segment vector of the directions of #» a and #» band the conditions #» a·#» wC= 0 and #» b·#» wC= 0 are satisfied. Hence, the following system of equations can be written:      #» wC=A(sC)−B(tC) #» a·#» wC=0 #» b·#» wC=0 (13) The solution of Eq. (13) is given by:            sC=(#» a· #» b)( #» b·#» w0)−� #» b�(#» a·#» w0) �#» a�� #» b�−(#» a· #» b)2 tC=�#» a�( #» b·#» w0)−(#» a· #» b)(#» a·#» w0) �#» a�� #» b�−(#» a· #» b)2 (14) where #» w0=NC1 0−NC2 0. For a cylinder–cylinder contact each cylinder is associated to a virtual sphere at the contact point. This approach is similar to the sphere–cylinder contact. The virtual spheres SC1and SC2are centred respectively in A(sc) and B(tc). Contact laws can then be inherited, and the force and moment linearly distributed on each cylinder node. 5 PFACET ELEMENTS The grid model described in Section 4 is extended to model arbitrary deformable objects. To achieve this objective a new element is introduced, the so-called PFacet (particle 5 749 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini facet). The PFacet element is geometrically constructed by the Minkowsky sum of a triangular facet and a sphere. It is composed of 3 nodes (spheres) and 3 connections (cylinders). The deformation is defined by the set of positions and orientations of its nodes. The constitutive behaviour of the PFacet element is simulated via an ordinary interaction law acting between the 3 nodes. The mass of the PFacet is equally lumped into the 3 nodes. The introduction of the PFacet element involves the consideration of the following three interactions: sphere–PFacet, cylinder–PFacet and PFacet–PFacet.        (a)       (b)       (c) Fig. 2. (a) Triangular skeleton with the nodes, (b) PFacet geometry where the blue zone corresponds to the inside of the PFacet element and (c) P′projection of the contact point Pon the plane that contains the PFacet nodes. (a) (b) (c) Fig. 3. (a) Sphere–PFacet interaction, (b) Cylinder–PFacet interaction, and (c) PFacet– PFacet interaction. For a sphere–PFacet interaction (Fig. 3a), the coordinates of the contact point Pare characterised by the position of the 3 PFacet nodes. The pointP′is the projection of the contact point P(Fig. 2c) on the plane containing the nodes of the PFacet element. The barycentric coordinates (p1,p2,p3) of the point P′(Fig. 2a) are defined as: 6 750 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini p1=� # » N1N3�( # » N1N2· # » N1P′)−( # » N1N2· # » N1N3)( # » N1N3· # » N1P′) � # » N1N2�� # » N1N3�−( # » N1N2· # » N1N3)( # » N1N2· # » N1N3) p2=� # » N1N2�( # » N1N3· # » N1P′)−( # » N1N2· # » N1N3)( # » N1N2· # » N1P′) � # » N1N2�� # » N1N3�−( # » N1N2· # » N1N3)( # » N1N2· # » N1N3) p3=1−p1−p2 (15) The barycentric coordinates of P′are used to distinguish if the contact is inside the PFacet element (the blue zone in Fig. 2b) or not. If p1≥0 and p2≥0 and p3<1, the contact is characterised as inside the PFacet element otherwise a sphere–cylinder interaction is considered. The virtual sphere s′ Pis positioned in P′. Both, the virtual sphere and the PFacet element have the same radius. The PFacet element radius is defined as the radius of the PFacet node located at the ends of each cylinder. The virtual sphere’s displacement and rotation are linearly interpolated between the 3 nodes of the PFacet N1,N2and N3. The translational (#» v′ P) and rotational (#» ω′ P) velocities are then calculated from the translational (#» vN1,#» vN2and #» vN3) and rotational (#» ωN1,#» ωN2and #» ωN3) velocities of the three nodes: #» v′ P=p1 #» vN1+p2 #» vN2+p3 #» vN3 #» ω′ P=p1 #» ωN1+p2 #» ωN2+p3 #» ωN3 (16) Introducing a virtual sphere s′ Pallows the use of the contact model for a classical sphere–sphere interaction. The contact force, #» F, is calculated for the virtual sphere and then the force is distributed over the nodes as follows: #» FNi=pi #» F(17) where #» FNiis the force applied at the node Niwith i∈{1,2,3}. A cylinder–PFacet contact (Fig. 3b) is decomposed into a node–PFacet ( i.e. sphere– PFacet) contact between the nodes of the cylinder and the PFacet, and a cylinder–cylinder contact between the cylinder and PFacet cylinders. A PFacet–PFacet interaction (Fig. 3c), is a combination of node–PFacet and cylinder– cylinder interactions. All these types of interactions have been described above. 6 VALIDATION A tube of length l(Fig. 4a) is subjected to a bending moment. The cylinder is fixed at one end, precisely the nodes initially positioned at x= 0 are only fixed along the x-axis. The bending moment Mbis applied at the other end (x=l). A distribution of forces {f1...f Nl}such that k=Nl k=1 fk= 0 and k=Nl k=1 xkfk=Mbis applied on each nodes which 7 751 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini are initially positioned at x=l. The forces fkare linear in yin the initial configuration. The forces are updated constantly during the deformation to reflect the positions xkin the deformed state. Mbis increased progressively in order to plot the load–displacement relation. The changes of Mbare sufficiently slow to maintain a quasi-static regime of deformation. The hollow cylinder is potentially able to buckle in bending. For the boundary conditions presented above, the relevant buckling limit is defined in [10]. It corresponds to an ovalisation of the cross section. The maximum bending moment is calculated as: Mmax b=2√2Kn 9(1 −ν2)d(18) where Knis the stiffness of the beam in tension, νis the Poisson’s ratio of the wall, and dis the thickness of the wall. This expression is used for comparison in Fig. 4e. In the calculation of Mmax b, the values of Knand νare obtained from a simulated tension test on the same tube. When a non-linear displacement is detected in the simulation, the applied moment is kept constant. The buckling is identified when, for a given moment Mb, the deformation of the cylinder increase with time. During the buckling phenomenon, the cylinder enters a dynamic mode of deformation. Finally, the tube reaches an equilibrium via a winding mechanism (Fig. 4d). In the simulation, buckling is reached for a moment value close to the theoretical limit. One can observe that the final state involves self-interaction in the tube (Fig. 4d). It shows that the method developed allows to investigate complex post-buckling behaviours. In addition, the capability of solving problems involving self interacting structures and large deformations is as well demonstrated. 7 EXAMPLE In the following section a potential application of the developed discrete model is presented. The ability to simulate deformable objects is investigated by simulating a impact of a sphere on a tube. The tube is deformable and composed of 640 PFacet elements with a radius of 2 mm. The tube has a radius of R=0.1 m and a length of 0.8 m. The material parameters used for the tube are: Young’s modulus Et= 1 MPa and Poisson’s ratio νt=0.3. The impacting sphere is first modelled by a solid sphere (Fig. 5a) and then by a hollow deformable sphere represented by 128 PFacet elements with a radius of 2 mm (Fig. 5a). The material parameters of the sphere are: Young’s modulus Es= 100 MPa and Poisson’s ratio νs=0.3 . The sphere has a radius R=0.1 m and an impact velocity of v=2m/s. To focus on the deformation of the impact gravity was not applied during the simulation. Fig. 5b and Fig. 5c illustrate the final state of the tube after impact of the solid sphere. The tube is deformed and the sphere is stopped with the creation of a crater. The deformation of the tube almost corresponds to its radius. Fig. 6b and Fig. 6c show 8 752 A. Effeindzourou, K. Thoeni, B. Chareyre and A. Giacomini (a) (b) (c) (d) 0 100 200 300 400 500 600 700 800 Torque (N/m) 0.000 0.001 0.002 0.003 0.004 0.005 Vertical deformation (m) Simulation Brazier [10] Linear solution (e) Fig. 4. Buckling of a tube subjected to an end bending moment: (a)–(d) deformation and (e) comparison with theoretical limit. the final state of the tube after impact of the deformable sphere. It can be seen that the tube is not as deformed as in the previous case and the impacting sphere deforms itself. 8 CONCLUSIONS The paper presents a method for the simulation of arbitrary deformable objects. The model is an extension of the flexible cylinder element described in [7]. A new element is introduced into the open-source framework YADE [8]. The so-called PFacet (particle facet) is an element geometrically constructed by the Minkowsky sum of a triangular facet and a sphere. Contact detection is based on the three primitives: sphere, cylinder and facet. It allows the representation of arbitrary deformable structures by using the standard contact models for spheres. The ability of the model to reproduce the behaviour of a structure under bending is validated. A tube subjected to an end bending moment is investigated and an excellent agreement with the theoretical buckling limit is observed. 9