A Consistent Time Level Implementation Preserving Second-Order Time Accuracy via a Framework of Unified Time Integrators in the Discrete Element Approach

1  Introduction

In order to numerically describe the behavior of granular materials (not continuous body) represented as a collection of large number of distinct particles, the discrete/distinct element method (DEM), proposed in [1], is being widely used for its ease of implementation and versatility. In the DEM technique, each particle element can undergo large displacement and rotation, and the deformations of individual particle elements are assumed to be small, compared with their displacements; hence we can treat the elements as rigid bodies and the deformation is usually described in terms of the relative normal displacement and the relative tangential displacement between the particles in contact. Consequently, the evaluation of those two displacements has a significant influence on both the dynamic response of particles and the numerical accuracy of the simulation.

One of the common applications for the DEM is particle packing. Various packing methods using DEM have been studied previously such as, Formation of a pile, i.e., sand piling process [2–6], pouring or depositing particles under gravity [7–12]. The DEM particle packing is also used in pharmaceutical field for powder tableting [13,14]. Another field of application using DEM is simulating particle flow; and a commonly studied example is ball mill mixer simulation; some of which include, horizontally rotating drums [15–18], tumbling ball mill [19–21]. Drum and mill are commonly used in industries for mixing, drying or coating subjects; thus, there is extensive interest among researchers in this study. Besides, the DEM is still studied and utilized in recent studies such as solid-liquid flow simulation [22], impact disruption of cohesive soft sphere simulation [23], granular shear flows of dry flexible fibers [24] and non-spherical complex granular simulation [25].

Furthermore, the DEM not only simulates granular flow, but it can also be used for solid body analysis and assessing fracture as studied in [26]. Recently, the DEM has been expanded to be able to solve multi-physics applications where granular particles interact with a continuous body, be it structural or thermodynamics. An Extended Discrete Element Method (XDEM) has been proposed in [27] which couples continuous body analysis method and DEM together without losing small scale information from commonly used volume averaging concept. The XDEM allows those microscale information from particles to accurately influence a macroscopic body. The XDEM opens up a wide variety of possible physical applications in many fields of study.

In transient dynamic simulation of the DEM type problems, contact forces and particle motion need to be evaluated and updated at every time step. The widely used time integrations in the DEM numerical simulations include the explicit Euler method, fourth-order Runge-Kutta (RK-4) method, and central difference method/velocity verlet method, etc. It is straightforward to implement these methods in the DEM simulation; however, these explicit time stepping techniques possess several numerical disadvantages in addition to the conditionally stable features; for example, the Euler methods are first-order accurate in time, RK-4 method requires additional computations of high-order derivatives of differential variables [28,29], and the central difference method is a non-dissipative scheme, i.e., it cannot introduce numerical dissipation if one needs to. Also, past practices poorly model the tangential physics of particle interactions. Recently, a unified time integration framework, under the umbrella of the so-called generalized single-step single-solve (GSSSS) family of algorithms, including not only implicit, but also explicit single-step single-solve algorithms, which can also be written in a linear three-step form, has been proposed in [30]. This single simulation tool kit, which has been originally developed in the field of computational structural dynamics, encompasses not only most existing algorithms that are second-order time accurate and are commonly used over the past 50 years or so, but it additionally inherits a wide variety of new implicit and explicit time integration schemes also with second-order time accuracy, including the explicit central difference method as a special case. To implement these time integrations in DEM and preserve the numerical accuracy of the algorithms, the contact force is required to be evaluated in a consistent way with the selected time integrations in both normal and tangential directions. According to our knowledge, most of the existing literature have been embedding the first-order accurate updating evaluation of tangential displacement into the first/second-order accurate time integrations. The specifics are shown in Section 2, leading to reduced accuracy.

In this work, the main objective is to show the consistent implementation of the GSSSS family of algorithms to the simulation of granular assemblies with DEM while retaining all the advantages of numerous methods within the present framework of algorithms. It is worthy to mention that one major drawback of the DEM technique, as well as the other particle based methods or meshless methods [31–35], is contact detection or neighboring particle search. To circumvent this issue, the linked list is employed to detect the contacting neighbor particles. In addition, Open Multi-Processor (Open MP) is implemented to reduce the simulation computational time. The outline of the rest of the paper is as follows. Section 2 presents the basic theory of DEM formulations and both the linear and the nonlinear model of contact forces and torques are presented in detail. The GSSSS unified time integrators are highlighted in Subsection 2.2 and the specifics of the proposed consistent tangential displacement evaluation is derived and compared with the existing tangential displacement evaluation technique carefully. The effectiveness of the proposed consistent tangential displacement evaluation is demonstrated and illustrated using simple numerical examples in Section 3. Finally, the conclusions are drawn and presented in Section 4. This study shows that the DEM with explicit schemes within the GSSSS family of algorithms achieves and consistently preserves second-order time accuracy; the consistent time level implementation and preserving the second-order time accuracy for the presented results for a wide variety of time integrators in the present unified framework include features such as no overshoot in displacement or no overshoot in velocity or no overshoot in both displacement and velocity for a given set of initial conditions of the particles and differs from algorithm to algorithm.

2  Discrete Element Method with the Explicit GSSSS Family of Algorithms

2.1 Basic DEM Formulations

DEM is a numerical simulation technique for a system of small “distinct” particles/elements, such as powders and granular materials. For the sake of simplicity, we assume the shape of every particle/element in a system is a disk in the two-dimensional Euclidean space E2 or a sphere in the three dimensional Euclidean space E3, although DEM can treat any shapes of elements in general, see [36]. In DEM, we assume that each particle is rigid, i.e., no deformation; however, particles are allowed to overlap slightly one another at the contact points.

Let qi(t):[0,∞)→Rndim denote the position component vector of the ith particle (i=1,2,⋯,nele) with respect to the Cartesian coordinate frame B, where nele and ndim denote the number of particles/elements in the system and number of dimension, respectively. Suppose that the ith and jth particles are in contact at time t as shown in Fig. 1a. The overlap (penetration) between particle i and j, which is denoted by nδij∈R, is defined as

nδij=(ri+rj)−dij(1)

where ri and rj denote the radii of the ith and jth particles, respectively; and dij is the center distance between the particles i and j and dij:=‖qj−qi‖=(qj−qi)⋅(qj−qi) is the center distance between the particles i and j. We assume the overlap is sufficiently small relatively to the particle sizes. If the inequality dij<ri+rj is met, the overlap nδij>0 occurs; and the contact forces between patricles i and j are taken into consideration. Define the unit normal vector between these particles as nij:=qj−qi‖qj−qi‖=qj−qidij,‖nij‖=1. The translational velocity component vectors of particles i and j (at the center points Ci∈Endim and Cj∈Endim, respectively) with respect to B are given as the (total) time derivatives of the position components vectors, i.e., vi=q˙i and vj:=q˙j, respectively (◻˙≡d◻/dt denotes the total time derivative of ◻). Let ωi and αi denote the angular velocity and angular acceleration vectors of particle i, respectively (for i=1,2,…,nele). Therefore, the velocity at point Pi∈Endim, as shown in Fig. 1a, vPi, can be written as vPi=vi+ωi×ri=vi+riωi×nij, where ri=rinij. Similarly, vPj, the velocity of point Pj, can be written as vPj=vj+ωj×rj=vj−rjωj×nij, where rj=−rjnij. Hence, the relative velocity vector at point Pi with respect to point Pj is obtained as

vPi/Pj=vPi−vPj=vi/j+(riωi+rjωj)×nij(2)

where vi/j:=vi−vj denotes the relative translational velocities of particle i with respect to j. Since the projection of vPi/Pj parallel to nij is given as

nvPi/Pj=(vPi/Pj⋅nij)nij=(vi/j⋅nij)nij,(3)

the projection of vPi/Pj perpendicular to nij, which is obtained by subtracting nvPi/Pj from vPi/Pj, can be written as

svPi/Pj=vPi/Pj−(vPi/Pj⋅nij)nij(4)

Figure 1: (a) Pictorial illustration of particles i and j in contact; and (b) the Kelvin-Voigt model for particles i and j in contact. O∈E3 is an origin of a cartesian coordinate frame B for E3 with a right-handed orthonormal basis for the associated vector space

Note that nvPi/Pj and svPi/Pj are the relative velocity vector of point Pi with respect to point Pj in the normal (nij) and tangential (tij) directions, respectively. Substituting Eqs. (2) into (4) yields:

svPi/Pj=vi/j−(vi/j⋅nij)nij+(riωi+rjωj)×nij (5)

Hence, tij=svPi/Pj||svPi/Pj||,‖tij‖=1.

Similarly, we can derive the relative acceleration vectors of point Pi with respect to point Pj in the normal (nij) and tangential (tij) directions. Since the accelerations at Pi and Pj are given as aPi=ai+ωi×(ωi×ri)+αi×ri=ai+riωi×(ωi×nij)+riαi×nij and aPj=aj+ωj×(ωj×rj)+αj×rj=aj−rjωj×(ωj×nij)−rjαj×nij, where ai∈Rndim and aj∈Rndim are the linear acceleration vectors of particles i and j, respectively, the relative acceleration vector at point Pi with respect to point Pj may be obtained as

aPi/Pj=aPi−aPj=ai/j+(riαi+rjαj)×nij+ri[(ωi⊗ωi)−(ωi⋅ωi)Indim]nij+rj[(ωj⊗ωj)−(ωj⋅ωj)Indim]nij(6)

where ai/j:=ai−aj denotes the relative translational accelerations of particle i with respect to j, and Indim∈Rndim×ndim denotes the identity tensor. Hence, the relative acceleration vector of point Pi with respect to point Pj in the normal (nij) and tangential (tij) directions are obtained as

naPi/Pj=(aPi/Pj⋅nij)nijandsaPi/Pj=(aPi/Pj⋅tij)tij(7)

respectively.

Contact Forces and Torques

The resultant contact force and torque acting on particle i are defined as the summations of the contact forces and torques, respectively, acting on particle i from all particles j(≠i) that are in contact with particle i. The contact forces due to particle(s) j can be modeled as follows.

Linear Model: The contact force acting on particle i from particle j can be uniquely decomposed into the contact forces in the normal direction (nij) and tangential direction (tij) as fij=nfij+sfij∈Rndim. In DEM, the contact forces are usually modeled by the Kelvin-Voigt model, as shown in Fig. 1b; hence, both the normal and tangential contact forces on particle i due to particle j have the stiffness/spring and damping terms, and the friction on the surface is also considered. For the normal contact force, nfij, using the stiffness and damping constants in the normal direction, kn and ηn, respectively, we have

nfij=kn(−nδij)+ηn(−nvPi/Pj)=−[kn(nδij)+ηn(vi/j⋅nij)]nij(8)

where nδij is the relative displacement vector between particle i and j in the normal direction, i.e., nδij∈Rndim, may be defined as the displacement while the particles i and j are in contact for the time interval, [tc1,tc2]⊂[0.∞), i.e.,

nδij=∫tc1tc2nvPi/Pjdt=nδijnij=[(ri+rj)−dij]nij(9)

The viscous damping coefficient ηn in Eq. (8) is usually selected as ηn=−2ln⁡CRm^ijknπ2+(ln⁡CR)2, where CR∈[0,1]⊂R+ denotes the coefficient of restitution, and m^ij, which is sometimes called the reduced mass, is defined as

m^ij:=[mi−1+mj−1]−1=mimjmi+mj(10)

Similarly, the tangential contact force, sfij, can be written as

sfij=ks(−sδij)+ηs(−svPi/Pj)(11)

where ks and ηs are the stiffness and damping coefficients in the tangential direction, respectively. The tangential relative velocity, svPi/Pj, in the equation above is given by Eq. (5). The relative displacement vector between particles i and j in the tangential direction, i.e., sδij∈Rndim, may be defined as the displacement while the particles i and j are in contact for the time interval, [tc1,tc2]⊂[0.∞), i.e.,

sδij=∫tc1tc2svPi/Pjdt(12)

We usually take ηs≃ηn unless otherwise the values of the damping coefficients are critical in the simulation. If the condition ||sfij||>μ||nfij||, where μ>0 denotes the coefficient of friction, is met, the “slip” between particles i and j on the surface at the contact point occurs, and therefore, the tangential contact force, given in Eq. (11), is replaced with the frictional force, sfij=−μ||nfij||tij with nfij, given in Eq. (8). Once we get all the contact forces, fij, acting on particle i, we can finally obtain the resultant contact force on particle i as fi=∑j≠incontifij, where nconti denotes the number of the particles that are in contact with particle i at time t. On the other hand, the resultant torque on particle i can be written in the form hi=∑j≠incontiri×sfij.

Nonlinear Model: There exist several nonlinear models developed for DEM in the literature to obtain more accurate, realistic simulation results. A common nonlinear model of DEM, proposed in [37], is based on the classical Hertzian contact theory and Mindlin theory [38,39]. The normal contact force, nfij, based on the Hertizian contact theory for two spheres, is still defined in the same form as Eq. (8),

nfij=−kn(nδij)−ηn(nvPi/Pj)(13)

However, as opposed to the linear model, nδij and stiffness kn are defined as

nδij=(nδij3/2)nij(14)

kn=43r^ijE^ij(15)

with nδij, defined as in Eq. (1), r^ij:=[ri−1+rj−1]−1=rirjri+rj, and E^ij:=[1−νi2Ei+1−νj2Ej]−1=EiEj(1−νi2)Ej+(1−νj2)Ei where Ei and νi are the Young modulus and Poisson ratio of particle i, respectively. For two spherical particles with the same radius (r=ri=rj) and the same material properties (E=Ei=Ej and ν=νi=νj), the stiffness defined in Eq. (15) may be written in the form, kn=E2r3(1−ν2). In the case of a spherical particle i with a wall, we assume rj→+∞, i.e., r^→ri; and thereby, kn=4ri3EiEwall(1−νi2)Ewall+(1−νwall2)Ei, where Ewall and νwall are the Young modulus and Poisson ratio of the wall. According to [37], the normal damping coefficient ηn in Eq. (13) is suggested to be ηn=κm^ijkn(nδij1/4), in which m^ij is as defined in Eq. (10). For κ>0, which can be regarded as dependent on the coefficient of restitution between two sphere particles or between a sphere particle and a wall, it is to be determined empirically by analysts (Usually, α∈[0.01,1]). Note that the nvPi/Pj in Eq. (13) is still defined as given in Eq. (3).

On the other hand, the tangential contact force, sfij, for a nonlinear model is defined, based on the Hertzian contact theory and Mindlin theory [38,39] with the assumption that “no-slip” is allowed between two sphere particles i and j, as

sfij=−ks(sδij)−ηs(svPi/Pj)(16)

with

ks=[18bij2−νiGi+2−νjGj]−1,bij:=r^ijnδij(17)

where Gi:=Ei/[2(1+νi)] is called the shear modulus of particle i, and bij is the radius of a contact area between sphere particles i and j. For two sphere particles with the same radius (r=ri=rj) and the same material properties (E=Ei=Ej and ν=νi=νj, i.e., G:=Gi=Gj), we have ks=2G2−ν2r(nδij). In the case of the contact between sphere particle i and a wall, we assume the wall (rj→+∞) is rigid, i.e., Ej→+∞, since the elastic tangential displacement of the wall due to the contact with particle i is usually negligible; hence, Eq. (17) leads to the following expression: ks=8Gi2−νiri(nδij). For most applications, we take the same value for the normal and tangential damping coefficients: ηs≃ηn.

2.2 Transient Analysis: A Consistent Time Level Implementation of Unified Time Integrators

Applying the Newton-Euler equation to particle i for i=1,2,…,nele, we get

Miai =fiappl(18)

Jiαi=hiappl(19)

where Mi=miIndim∈Rndim×ndim, with the identity tensor Indim∈Rndim×ndim, and ai∈Rndim denote the mass tensor and the acceleration vector of particle i, respectively; and Ji∈Rndim×ndim and αi∈Rndim denote the inertia tensor (with respect to the center of mass) and the angular acceleration vector of particle i, respectively. The right-hand side of Eq. (18), i.e., fiappl∈Rndim, is the applied force on particle i, and it includes any force acting on the particle, including conservative and nonconservative forces, as well as the contact forces described in the previous subsection. Likewise, hiappl∈Rndim in Eq. (19) is the applied torque vector of particle i.

In a unified format for all particles (i=1,2,…,nele) in the system, Eqs. (18) and (19) can be written in the form Ma=fappl and Jα=happl, respectively, where M=diag(M1,M2,…,Mnele)∈RN×N, J=diag(J1,J2,…,Jnele)∈RN×N, a=(a1T,a2T,…,aneleT)T∈RN, α=(α1T,α2T,…,αneleT)T∈RN, fappl=(f1applT,f2applT,…,fneleapplT)T∈RN, h=(h1applT,h2applT,…,hneleapplT)T∈RN with N:=ndimnele. Or simply,

[M00J]⏟∗M(aα)⏟∗a=(fapplh)⏟∗f∀t∈[0,∞)(20)

Time Discretization

Consider the temporal discretization of the governing equations for a time interval I=[t0,tf]⊂[0,∞) split into subintervals, i.e., I=[t0,tf]=⋃n=0f−1[tn,tn+1], in which t0 and tf denote the initial and final time levels of simulation, respectively. The time step size is defined as Δt:=tn+1−tn>0, and it is assumed to be constant for simplicity. Temporally discretizing Eq. (20) by way of the normalized time weighted residual methodology [40] results in the following algorithmic framework:

∗M∗a~n+1=∗f(∗q~n+1,∗v~n+1,tn+W1)(21)

where the algorithmic acceleration, velocity, and configuration vectors are defined as

∗a~n+1=∗an+W1Λ6Δ∗an+1∗v~n+1=∗vn+ΔtW1∗an,∗q~n+1=∗qn+ΔtW1∗vn+Δt22W2∗an(22)

respectively, for n=0,1,…,f−1. The associated update equations to advance from t=tn to t=tn+1 are

∗qn+1=∗qn+Δt∗vn+Δt22∗an+λ3Δ∗an+1Δt2∗vn+1=∗vn+Δt∗an+Δtλ5Δ∗an+1,∗an+1=∗an+Δ∗an+1(23)

The algorithmic parameters, also known as the DNA [discrete numerically assigned] markers, are defined with the input scalar algorithmic parameters (ρ∞min,ρ∞max,ρ∞s), denoting the minimum principal root, the maximum principal root, and the spurious root [41], respectively, together with η, as

W1=W2=11+ρ∞s,W1Λ6=2+ρ∞min+ρ∞max+ρ∞s−ρ∞minρ∞maxρ∞s(1+ρ∞min)(1+ρ∞max)(1+ρ∞s)λ3=η(1+ρ∞min)(1+ρ∞max),λ5=3+ρ∞min+ρ∞max−ρ∞minρ∞max2(1+ρ∞min)(1+ρ∞max)(24)

for the U0(ρ∞min,ρ∞max,ρ∞s) family, and

W1=3+ρ∞min+ρ∞max−ρ∞minρ∞max2(1+ρ∞min)(1+ρ∞max),W2=2(1+ρ∞min)(1+ρ∞max)W1Λ6=2+ρ∞min+ρ∞max+ρ∞s−ρ∞minρ∞maxρ∞s(1+ρ∞min)(1+ρ∞max)(1+ρ∞s),λ3=η2(1+ρ∞s),λ5=11+ρ∞s(25)

for the V0(ρ∞min,ρ∞max,ρ∞s) family. The U0 and V0 represent zero-order overshooting behaviors in the configuration and the velocity fields, respectively. Note that no two schemes within the algorithmic framework can have the same DNA markers that completely distinguish one scheme from another. Any scheme within this framework is guaranteed to be second-order accurate by approximating ∗qn≈∗q(tn), ∗vn≈∗q.⁡(tn), and ∗an≈∗q..⁡(tn−ϕ) at the time step n where the acceleration time level shifting parameter is defined by ϕ:=W1Λ6−W1; see [42] for the details. At time t=t0, the acceleration need to be computed from ∗a0=∗M−1∗f(∗q0,∗v0,t0) where ∗q~0=∗q(t0) and ∗v~0=∗q.⁡(t0) are given initial conditions. In DEM simulation, we update from n to n+1 via Eqs. (21)–(23) after evaluating the force vectors, as shown in the previous subsection, at time tn.

Algorithmic Evaluation of sδij:

Suppose particles i and j are in contact during a time interval Ic=[tc1,tc2]⊂I and split it into subintervals, Ic=[tc1,tc2]=⋃n=c1c¯[tn,tn+1] with tc¯+1≡tc2. Following the temporal discretization as given in Eq. (21), the tangential contact force, sf~ij, is evaluated at q=q~n+1, v=v~n+1, and ω=ω~n+1. That is, the algorithmic tangential force may be written as

sfij|q~n+1,∗v~n+1=:sf~ijn+1=ks(−sδ~ijn+1)+ηs(−sv~Pi/Pjn+1)(26)

where sv~Pi/Pjn+1 and sδ~ijn+1 are the algorithmic relative tangential velocity vector and the algorithmic relative tangential displacement vector defined as follows: sv~Pi/Pjn+1:=v~i/jn+1−(v~i/jn+1⋅n~ijn+1)n~ijn+1+(riω~in+1+rjω~jn+1)×n~ijn+1, where v~i/jn+1=vi/jn+ΔtW1Λ4ai/jn, n~ijn+1=q~j−q~i‖q~j−q~i‖, and

sδ~ijn+1=‖s⧫δ~ijn+1‖t~ijn+1(27)

where

s⧫δ~ijn+1:=sδijn+ΔtW1Λ1svPi/Pjn+Δt2W2Λ2saPi/Pjn(28)

t~ijn+1=sv~Pi/Pjn+1||sv~Pi/Pjn+1||(29)

Note that sδijn in Eq. (28) is given by

sδijn=‖s⧫δijn‖tijn(30)

where

s⧫δijn:=sδijn−1+ΔtsvPi/Pjn−1+Δt22saPi/Pjn−1+Δt2λ3ΔsaPi/Pjn,ΔsaPi/Pjn:=saPi/Pjn−saPi/Pjn−1(31)

tijn:=svPi/Pjn||svPi/Pjn|| (32)

Remark

1.    If the contact starts at time tn, i.e., t=tn=tc1, we have sδijn=0 in Eq. (28).

2.    When sv~Pi/Pjn+1=0, the algorithmic tangential unit vector, t~ijn+1, defined by Eq. (29), becomes singular. In this case, we assume that particle j is sliding on the surface of particle i, and therefore use t~ijn+1=s⧫δ~ijn+1‖s⧫δ~ijn+1‖ instead; and we have sδ~ijn+1= s⧫δ~ijn+1. Similarly, if svPi/Pjn=0, use tijn=s⧫δijn‖s⧫δijn‖ instead of Eq. (32).

3.    When slip occurs at the algorithmic time level tn+W1∈Ic, i.e., ‖sf~ijn+1‖>μ‖nf~ijn+1‖, in which nf~ijn+1 is the algorithmic normal contact force defined as nf~ijn+1:=−kn[(ri+rj)−‖q~jn+1−q~in‖]n~ijn+1−ηn(v~i/jn+1⋅n~ijn+1)n~ijn+1, and the algorithmic tangential force sf~ijn+1 defined in Eq. (26) is replaced with the following algorithmic frictional force, sf~ijn+1=ks(sδ~ijn+1)=−μ‖nf~ijn+1‖t~ijn+1. Therefore, the algorithmic tangential displacement vector, sδ~ijn+1, defined in Eq. (27), should be replaced with sδ~ijn+1=−μks‖nf~ijn+1‖t~ijn+1. Similarly, in the case of ‖sfijn‖>μ‖nfijn‖, employ sδijn=−μks‖nfijn‖tijn instead of Eq. (30).