Efficient implementation of superquadric particles in Discrete Element Method within an open-source framework

Particle shape representation is a fundamental problem in the Discrete Element Method (DEM). Spherical particles with well known contact force models remain popular in DEM due to their relative simplicity in terms of ease of implementation and low computational cost. However, in real applications particles are mostly non-spherical, and more sophisticated particle shape models, like superquadric shape, must be introduced in DEM. The superquadric shape can be considered as an extension of spherical or ellipsoidal particles and can be used for modeling of spheres, ellipsoids, cylinder-like and box(dice)-like particles just varying five shape parameters. In this study we present an efficient C++ implementation of superquadric particles within the open-source and parallel DEM package LIGGGHTS. To reduce computational time several ideas are employed. In the particle–particle contact detection routine we use the minimum bounding spheres and the oriented bounding boxes to reduce the number of potential contact pairs. For the particle–wall contact an accurate analytical solution was found. We present all necessary mathematics for the contact detection and contact force calculation. The superquadric DEM code implementation was verified on test cases such as angle of repose and hopper/silo discharge. The simulation results are in good agreement with experimental data and are presented in this paper. We show adequacy of the superquadric shape model and robustness of the implemented superquadric DEM code.


Introduction
In many engineering applications different types of particles have to be stored, transported, mixed, or segregated. Despite that, knowledge of the static and dynamic behavior of particulate solids is still not fully understood. Such knowledge is of major importance for a proper design of processing units of silos, rotating drums, and others [23]. The Discrete Element Method (DEM), developed by Cundall and Strack [14], has proven to be an efficient tool for modeling granular materials. In DEM granular material is treated as a system of distinct interacting particles. Each particle has own mass, velocity, position, and contact properties; it obeys Newton's second law and is tracked individually. Together with the rapidly increasing computational power available, DEM becomes more and more popular among engineers and researchers. A comprehensive overview of major DEM applications can be found in [56].
Many DEM codes still employ disks (in 2D) and spheres (in 3D) to represent particle shapes due to their implementation simplicity and efficiency in speed of contact detection, which results in faster code development and lower computational time. The rolling friction correction can be theoretically linked to various physical effects to model particle non-sphericity using spherical particles, as emphasized by Ai et al. [2]. Moreover, the contact force models that include normal and tangential forces for a pair of interact-ing spheres are already established. An overview of the most popular contact forces is given in [55]. However, particles in granular and powder materials in nature and industry are mostly non-spherical. Moreover, spherical particles behave differently than complex-shaped particles, not only on the single grain level but also as an assembly. As summarized by Cleary [10], non-spherical particles differ from the spherical ones in the following ways: compacity of packed heap, resistance to shear and roll, and, as a result, ability to block the flow. Therefore, the physical validity of results obtained from simulations using spherical particles is usually questionable [31,48].
Many approaches have been suggested in the literature to handle particle non-sphericity. Previously, Lu et al. [34] have summarized the main theoretical developments in nonspherical DEM and reviewed its applications. The most popular approach in the DEM community, according to Lu et al. [34], is the multisphere (MS) approach [1,19,30,31,37,50]. In this method simple spheres are allowed to overlap and glued together to represent complex shapes. The method has the advantages that any shape can be represented by a set of glued (or prime) spheres and contact detection together with force calculation is based on that for spheres. One of the disadvantages is the fact that high accuracy of the shape approximation requires a significant number of prime spheres. Markauskas et al. [37] showed that approximation of ellipsoidal particles by 25 prime spheres increases CPU time by factor of 17. Marigo and Stitt [36] studied the influence of particle shape representation by the MS approach for a system of cylindrical pellets and found that about 160 primary spheres are required to be in agreement with experimental data. Another disadvantage of the MS method is the occurrence of multiple contact points [31] since an approximation of a convex particle (e.g., ellipsoid) by MS-particle is always non-convex if the number of primary spheres is more than one. The number of interparticle contacts increases with the increase of the number of the prime spheres [37]. Thus, a reasonable number of prime spheres should be chosen.
Polygonal (in 2D) and Polyhedral (in 3D) particles as introduced by Cundall [13] have been widely used in DEM to model granular materials. In this approach the particle surface is approximated by line segments (in 2D) or by triangles (in 3D), thereby providing a high level of versatility in particle shape representation. Different algorithms for collision detection were developed by Cundall [13], Chang and Chen [8], Boon et al. [5], and Nezami and Hashash [41]. The major drawback of this method is that the question of how contact forces between two colliding polyhedral bodies are calculated is still not completely answered.
The superquadric shape is an extension of spheres and ellipsoids. This shape was first introduced in mathematics by Barr [3], used in DEM by Williams and Pentland [52] in 2D DEM and by Cleary in 3D DEM simulations [9,12] and later used by Lu et al. [35]. The superquadric equation given by Barr is as follows: where a, b, c are the half-lengths of the particles along its principal axes, and n 1 and n 2 are blockiness parameters. Parameters n 1 = n 2 = 2 give an ellipsoid, and a cylinder is obtained if n 1 = 2 and n 2 2 and a box-like particle if n 1 2 and n 2 2 ( Fig. 1). Superquadrics give an excellent trade-off between model complexity and shape flexibility by simply changing 5 shape parameters (a, b, c, n 1 , n 2 ) in formula (1). However, the use of superquadric particles is limited in that sense that only ellipsoidal, box-like and cylinder-like particles can be modeled. Another disadvantage of the superquadric approach is that the contact detection procedure can be implemented, possibly, only by using typically computationally expensive iterative methods (Newton's method), convergence properties of which decrease with increase of blockiness parameters n 1 and n 2 .
Unfortunately, there is a lack of literature that could provide detailed descriptions of algorithms necessary for implementation of superquadric particles in DEM. In this study we will present the non-spherical DEM approach providing all necessary mathematical tools for an efficient implementation of superquadric particles in the DEM based on open-source DEM package LIGGGHTS [28], such that the reader can understand the underlying algorithms and analytical expressions for particle-wall contact and minimum bounding sphere. We show good versatility of the approach for the practical range of blockiness parameters. Validation work along with several application examples is presented.
2 Numerical model 2.1 Motion of an arbitrarily shaped particle In DEM each particle i obeys Newton's second law and is tracked individually by solving explicitly their trajectories: where m i and X Ci are the mass and the position of a particle center, F i is the total force acting on a particle i that is the sum of normal particle-particle, tangential particle-particle and external non-particle forces like gravity. The contact force and the contact point between two non-spherical particles depend on particles' orientation, and hence accurate determination and integration of orientation becomes critical. Accurate determination of a particle's orientation is also critical for the determination of the angular velocity of a particle. Orientation of a non-spherical particle is usually described as a rotation of the coordinate vectors {e 1 , e 2 , e 3 } that define the global observer-fixed reference frame to the coordinate vectors {ê 1 ,ê 2 ,ê 3 } that define the local body-fixed reference frame. This rotation can be tracked by rotation vectors [7] or by quaternions [22] that are singularity-free in contradistinction to the methods based on Euler angles. The quaternion of rotation q can easily be constructed from the unit axis of rotation e and the angle of rotation α around this axis: The rotation matrix A = A(q) is constructed from the quaternion components: By definition of the rotation matrix: A · e 1 =ê 1 , A · e 2 =ê 2 , A · e 3 =ê 3 , A −1 = A T . The orientation of a particle can be updated every time step using the following expression [32]: where sign "•" denotes quaternion multiplication [22], ω i = (ω i x , ω iy , ω iz ) T = A −1 i Ω i is the angular velocity in the particle-based (canonical) coordinate system, and Ω i is the angular velocity in the observer-fixed coordinate system. Rotational motion of a particle is tracked by the following equation in the observer-fixed coordinate system: where L i = I i · Ω i is the angular momentum of the particle i, I i is the tensor of inertia, and T i is the total torque acting on a particle i with respect to the particle center. Note that the normal force can also produce torque and must be taken into account while calculating T i . For a spherical particle Eq. (6) above can easily be resolved with regard to Ω since its tensor of inertia is always constant and only possesses nonzero elements in the diagonal: where E is the identity tensor. In general case the t A ·Î · A −1 , wherê I is the principal tensor of inertia, i.e., the tensor of inertia in particle-based coordinate system which contains non-zero entries only in the diagonal:Î x ,Î y andÎ z . Moving to the body-fixed coordinate system yields the following general form of Eq. (6): where t i = (t i x , t iy , t iz ) T = A −1 i T i is the total torque acting on a particle in the particle-based coordinate system. Equation (7) in conjunction with Eq. (5) can be solved by various methods, such as described by Miller et al. [38], Walton and Braun [51] and Omelyan [42][43][44][45]. Corresponding analytical expressions for volume and principal moments of inertia can be found in works by Jaklič et al. [25,26].

Neighbor search
In this section we describe the broad phase of the contact detection algorithm. To perform simulations of large-scale systems, it is essential to optimize the computational strategy. The number of potential contact pairs can be minimized by employing a neighbor list to exclude particle pairs that are a priori too far from each other to be in contact [28]. Different techniques for constructing neighbor lists have been proposed in the literature. These include the Verlet-Neighbour List [49], Linked Cell method [46], NBS algorithm [39], and MR linear contact detection algorithm [40]. The number of potential contact pairs in a neighbor list can be reduced with the help of bounding volumes. A bounding volume is a simple volume that encapsulates a more complex body. According to Ericson [18] doing cheap bounding volume intersection tests before performing more complex geometric tests results in a significant performance gain since the amount of work needed to determine a collision is reduced and computational time is saved by rejecting contact pairs whose bounding volumes do not intersect. Such bounding volumes include bounding spheres and oriented bounding boxes (OBB) (Fig. 2).
The bounding sphere is the most memory-efficient bounding volume. The spheres' intersection check consists in evaluating if the distance between the centers of the spheres is less or equal than the sum of their radii. If the answer is yes, than the narrow phase of the contact detection must be used. In this section we present an accurate analytical solution for the minimum bounding sphere, i.e., the bounding sphere of minimum volume, for a superquadric particle using its implicit shape equation. We seek the point (x, y, z) on the particle's surface which has the largest distance to the particle center. This distance is the minimum bounding radius, which gives the minimum volume. This condition in terms of an optimization problem can be expressed as follows: where f (x, y, z) is the superquadric equation in the bodybased coordinate system (Eq. (1)). Without loss of generality we require x > 0, y > 0, z > 0. Using Lagrange multipliers this optimization problem can be rewritten as a system of non-linear equations: x a n 2 −1 x a Substituting x = ax, y = αbx, z = βcx gives: Doing one more substitution γ = (1 + α n 2 ) n 1 /n 2 −1 yields: The solution of this system provides the following: The minimum bounding radius can now be easily obtained: r = x 2 + y 2 + z 2 . For spherical and ellipsoidal particles (n 1 = n 2 = 2), the system of equations degenerates, which has four solutions: The minimum bounding radius becomes r = max(a, b, c). For the case of the cylinder (n 1 > 2, n 2 = 2), without loss of generality, let a > b. This will give α = 0. Other unknowns can be found using the equations above.
While the bounding sphere is an orientation invariant approximation of a particle, the oriented bounding box (OBB) can capture orientation and aspect ratio of the particle. The minimum oriented bounding box is a rectangular block with semi-axes a, b, c with the center located at the center of the particle in question and oriented as the particle. The intersection check methods between two OBBs are usually based on the concept of the separation axis and can be found in [17,18,47].

A contact detection algorithm
Equation (1) defines a superquadric surface in its local (canonical) coordinate system. We will refer to the function f (x) ≡ x a n 2 + y b n 2 n 1 /n 2 + z c n 1 − 1 as the shape function. If for a certain point (x, y, z) T the value f < 0 then the point is located inside the particle, if f > 0, then (x, y, z) T is outside the particle. If f = 0, then (x, y, z) T lies on the particle surface.
For practical use it is necessary to be able to define a superquadric particle with respect to a global coordinate system which is usually observer-fixed. This is done by applying the usual translation and rotation operations. The shape function F of a superquadric particle with center X C and quaternion q in a global frame is given by the following expression: The points X and X C are defined in the global frame. For a contact detection algorithm it is also necessary to be able to calculate 1st (gradient) and 2nd (Hessian matrix) derivatives of the shape function. The gradient of the shape function calculated at a point x = (x, y, z) T in the local frame is the following: The second derivatives are given by x a where ν = x a . The corresponding gradient vector and Hessian matrix read The first and second derivatives of the shape function at a point X in the global coordinate system now can easily be established by calculating them in the local frame and applying the transition formulas: Now we are ready to formulate the contact detection problem in terms of an optimization problem. Consider two superquadric particles with two times continuously differentiable (automatically fulfilled if n 1 ≥ 2 and n 2 ≥ 2) shape functions F 1 (X) and F 2 (X) defined in global frame. Following Houlsby [24] and extending the algorithm for 3D case we seek a "midway" point P between the particles and "closest" to both. In other words, solve the following optimization problem: Applying the Lagrange multipliers approach gives the Lagrange function in the following form: The equation ∇ X,λ L(X, λ) = 0 gives the condition for the stationary point: Regrouping terms yields Introducing μ 2 = (1 − λ)/(1 + λ) brings us to the following system introduced by Cleary et al. [12]: This system of 4 equations with 4 unknowns has to be solved at each DEM time step for each pair of particles. This system gives a mathematical condition for interparticle contact detection. If for point X 0 conditions F 1 (X 0 ) < 0 and F 2 (X 0 ) < 0 are fulfilled (Fig. 3), the contact between two particles takes place with the contact point X 0 and the overlap direction n 12 = ∇ F 1 /||∇ F 1 || or n 12 = −∇ F 2 /||∇ F 2 || calculated at the contact point. Newton's method for this system can be expressed as is the Jacobian matrix of the right-hand side term Φ: For stability reasons it is necessary to find a scalar parameter at every iteration to ensure convergence of the algorithm. There are several methods to obtain such a scalar parameter. One of them is first to check if α = 1 satisfies ||Φ(Z n+1 || < ||Φ(Z n )||. If not, let α := α/2 and repeat or use the Golden section algorithm [27] with termination if any α satisfying ||Φ(Z n+1 )|| < ||Φ(Z n )|| is found. The solution from previous DEM time step can be used as a starting point (initial guess). Usually a few Newton iterations (1 − 3) are required to converge, depending on a user defined tolerance ε 1 for termination criterion ||Φ(Z n+1 )|| < ε.
If there is no information on the contact point from the previous step, let (a 1 0 , b 1 0 , c 1 0 , n 1 10 , n 1 20 ) and (a 2 0 , b 2 0 , c 2 0 , n 2 10 , n 2 20 ) be the shape and blockiness parameters of particle 1 and 2. The following approach is suggested: (1) Find the contact point for a pair of volume equivalent spheres with radii r 1 and r 2 and centers located at the same points as for given particles, X 1 C and X 2 C . These spheres defined as superquadrics have the following shape and blockiness parameters: The analytical solution for the sphere-sphere contact point is Use this point as a starting point.
(2) Choose number of steps N and calculate 3. Modify shape and blockiness parameters.
(4) Calculate the contact point for particles with shape parameters (a 1 , b 1 , c 1 , n 1 1 , n 1 2 ) and (a 2 , b 2 , c 2 , n 2 1 , n 2 2 ) using the iterative algorithm described above and the last computed starting point. Use the found contact point as a new starting point for the next step. Steps 1-5 of the iterative procedure listed above define the step-wise linear transition of the spherical shape parameters to the shape parameters of the given particles 1 and 2 where k is the iteration number. After k = N steps the shape and blockiness parameters will be the same as initial ones by construction of the procedure: As a result, this procedure will give the contact point for the given pair of particles. This step-wise procedure ensures the convergence of the method and does not affect computational time significantly, since it must be called only once per pair of particles. We found that for the exponents n 1 ,n 2 ≤ 8 convergence is guaranteed. In addition, Cleary and Sawley [11] studied influence of the blockiness parameters on the mass flow rate from a hopper and found that the values n 1 , n 2 > 8 fail to make any difference to the nature of the flow. Thus, we can safely recommend to use blockiness parameters in the range between 2 and 8.

Particle-wall contact
For industrial applications of DEM it is necessary to be able to resolve the contact between a particle and a flat surface. A flat surface (the wall) is usually defined by any point x w on it and the unit normal vector n defined in the particle-based coordinate system. This yields the equation of the wall: To establish the contact point we first seek a point x = (x, y, z) T on the particle surface that has the minimum/ maximum distance to the wall (Fig. 4). The mathematical condition in terms of an optimization problem is expressed by maximize n · x subject to f (x, y, z) = 0.
The optimization problem has two solutions with different signs for (x, y, z, λ). Without loss of generality the normal vector n is directed outwards with respect to the particle and its components n x , n y , and n z are positive. Hence, we seek a point x, y, z with positive signs. Applying the Lagrange multipliers approach we can rewrite the optimization problem as follows: Performing the same variable change as in the previous section gives the following system: The solution of this system provides the following expressions: For the outer normal n = (n x , n y , n z ) T with components of any signs the solution can be generalized: The normal overlap vector δ between the particle and the wall now can be easily established by calculating the projection x * of the contact point onto the wall: To calculate the overlap vector and the contact point in the global coordinate system, transition formulas (see Eq. (18)) can be applied. A corresponding algorithm for interaction between superquadric particles and walls of arbitrary shape has been developed in LIGGGHTS. This algorithm employs the solution for the contact between a particle and a flat wall presented above. However, the description of this algorithm is beyond the scope of this paper.

Contact force calculation
In the spherical Discrete Element Method the two following approaches are common: the hard-sphere and the soft-sphere approach. In the hard-sphere approach (eventdriven), particles are assumed as rigid bodies, a sequence of collisions is processed, one collision at a time without the contact forces being explicitly considered. In the soft-sphere approach (time-driven) particles are allowed to deform slightly(overlap), and the contact forces are calculated as functions of the overlap [55]. This overlap is not real but intends to model the deformation of the interacting particles at a contact point in an indirect way.
Di Renzo and Di Maio [15] suggested using the linear spring Hertzian model without cohesion for the particleparticle and particle-wall contacts. This model employs the following formula for the interparticle contact force, acting from a spherical particle i with radius R i and center at point X Ci on a spherical particle j with radius R j and center at point X C j : where • F n,i j = k n,i j δ n,i j + γ n,i j v n,i j is the normal force component, • F t,i j = k t,i j δ t,i j + γ t,i j v t,i j is the tangential force component, • δ n,i j = δ n,i j n i j is the normal overlap vector, i j dτ is the tangential overlap [28].
Both normal and tangential forces contain a spring force and a damping force. The coefficients k n , k t , γ n , γ t are calculated from the material properties (density, coefficient of restitution, Poisson ratio, Young's modulus, shear modulus), overlap, and radii of the particles. The corresponding expressions can be found in [4].
However, the models above are only applicable for spherical particles. Feng and Owen [20] proposed theoretical framework for developing energy-conserving normal contact models for arbitrarily shaped particles. It has been established that the normal force must be a potential field vector associated with a potential function φ that is a function of the overlap volume [20]. However, the accurate calculation of the overlap volume may be computationally expensive and thus become not applicable for a case with millions of particles.
Previously, Zheng et al. [54] modified the Hertzian model taking into consideration two principal radii of curvature and applied it to ellipsoidal particles. They showed good agreement with the results calculated by means of the finite element method (FEM). However, the extension of the model to superquadric particles becomes quite complex. In this study we propose the following approach. The contact point X 0 and the contact direction n i j define the contact line. We seek points X i and X j as the nearest (with respect to X 0 ) intersection points between the contact line and the particles' surfaces. In other words, we solve the following non-linear algebraic equations separately with respect to the scalars α i > 0 and α j < 0 at every DEM time step for each pair of overlapping particles: Then the normal overlap vector is These equations with respect to α i and α j are easier to solve if moved to their own local reference frames: Note that in both reference frames scalars α i and α j are the same for each particle. The equations can be easily solved by Newton's iterations: Calculation of the coefficients k n , k t , γ n , γ t in(35) requires knowledge of the equivalent radius R eq [15]: where R i and R j are the radii of particles i and j. We will use R = 1/K as the particle radius, where K = K mean = 1 2 (κ 1 + κ 2 ) is the mean local curvature coefficient [21] calculated at X i and X j for each particle correspondingly, and κ 1 and κ 2 are the principal curvature coefficients: Alternatively, one can use the Gaussian curvature coefficient K = K Gauss = √ κ 1 κ 2 [21]: Obviously, K mean = K Gauss = 1/R for a spherical particle of radius R.
One of the disadvantages of the proposed methodology (and hence of the possible extension of the method by Zheng et al. [54] to superquadrics) is that the curvature coefficients may become zero leading to infinite curvature radii if superquadric exponents n 1 and n 2 are more than 2, especially in face-to-face contact. The mean curvature radius becomes infinite if both principal curvature coefficients are zero. This occurs at 6 points on the particle surface: x = y = 0, z = ±c, y = z = 0, x = ±a and x = z = 0, y = ±b (in the local reference frame). The Gaussian curvature radius becomes infinite if any of the principal curvature coefficients is zero. This occurs if x = 0, y = 0 or z = 0 (in the local reference frame). For this reason, we limit the curvature radius: R curvature = min(R curvature , q R vol ), where R vol is radius of the volume equivalent sphere, q is the limiting coefficient that must be chosen in advance. In the current implementation of LIGGGHTS q = 10 is used. The influence of the choice of the curvature radius and the limiting coefficient q on the simulation results is to be studied in the future publications.

Contact force between two ellipsoidal particles
Here we bring in contact two ellipsoidal particles (Table 1) oriented parallel to each other (Fig. 5). Particles are consid-  [54] ered to be elastic and frictionless. The overlap distance was varied in the range between 0 and 5μm. Three radius models were tested: the mean curvature radius, the Gaussian curvature radius, and radius of the volume equivalent sphere. The normal contact force between overlapping ellipsoidal particles is plotted as a function of the overlap and presented in Fig. 6. Calculated results are compared with the FEM analysis carried out by Zheng et al. [54] and added to Fig. 6. It can be observed that for elastic contacts the normal force can be well described by the Hertz theory using the Gaussian curvature radius for the particle radii R i and R j in Eq. (38) since it gives minimal discrepancy with respect to the FEM results. Despite this, only the mean curvature radius is used in the further test cases in this paper since it becomes infinite only at 6 points of the non-ellipsoidal (superquadric exponents n 1 , n 2 > 2) particle surface rather than on the infinite number of points if the Gaussian curvature is used, as described in the previous section. A more comprehensive comparison of different curvature radii using different particle shapes at different orientations is to be discussed in a future paper.

Settling of particles under gravity and simulation speed
In the following test case we compare particles (Table 2) with different blockiness parameters in terms of computational time. The code was compiled with g++ (5.2) compiler and run in serial mode on an Intel Core i7-4790 processor-based desktop computer. A total of 1000 particles were distributed randomly (Fig. 7, left) in the simulation box 0.1 × 0.1 × 0.25 m and allowed to settle under gravity along Z-direction and form a static packed bed (Fig. 7, right). Periodic boundary  However, in order to make an accurate comparison between different methods, the program codes must be run on the same hardware.

Particle-wall impact
In this test, as described by Kodam et al. [29], a cylindrical particle oriented at a specified angle, impacts a wall with a specified translational speed normal to the wall and zero angular velocity. The contact is assumed to be frictionless and without gravity. The post-impact angular and translational velocities, ω + y and V + z correspondingly, according to Kodam et al. [29], are given by where m is the mass of the cylinder, ε is the coefficient of restitution at the point of contact, V − z is the pre-impact translational velocity of the cylinder, α is the angle between the cylinder's face and the line joining the contact point and the center of the particle, θ is the angle between the cylinder's face and the wall, and I yy is the moment of inertia around the y-axis. The parameter r is the distance between the cylinder's center and the corner point C (Fig. 9), which is assumed to be fixed. Particle parameters are listed in Table 3.
The post-impact angular and translational velocities were calculated for various orientation angles θ for the DEM simulations and compared with analytical expressions in Fig. 10. The wall was removed immediately after collision to prevent the secondary contact that occurs in reality at low and high impact angles and is not taken into account in Eqs. (41) and (42).
For angles 5 • θ 86 • the superquadric DEM gives relative good agreement with analytics; however, for other angles there is small disagreement mainly due to the error in shape approximation, and, as a result, because of the corner point C, that is non-static with respect to the particlebased coordinate system. This corner point is always static at impact angles θ = 0 • and θ = 90 • for true cylinders. According to [29], multisphere simulations (with 54 prime spheres) show significant errors over most of the orientation range.

Piling of particles
For the second validation test case superquadric particles with the following shape parameters were used: a = 2.0 mm, b = 2.0 mm, c = 1 mm, n 1 = n 2 = 4. Particle parameters along with simulation setup data are listed in Table 4. They were compared with volume equivalent spherical particles of radius R = 1.836 mm. Domain boundaries are represented by rigid walls of the same material as the particles. In both  cases the heap was formed by continuously dropping particles from a small area located above the center of the heap. As can be seen from Fig. 11 the heap shape for spherical particles (angle of repose 31 • ) differs from the heap shape for superquadric particles (angle of repose 40 • ). The heap becomes stable 1s after the dropping is stopped with almost zero maximum angular and translational velocity, which testifies the stability of the algorithms. These simulations show importance of using non-spherical particle shapes in the Discrete Element Method. Having the same material properties the non-spherical particles can demonstrate different behavior in comparison to spherical ones just by changing the shape of the particles.

Angle of repose
In this test case the following particles were used: sugar cubes, "M&M's" and chewing gum. The particles are randomly distributed in a cylindrical volume with random orientation. They are allowed to settle under gravity in the cylindrical volume. After the settling is completed, the vertical wall of the cylinder moves upwards such that some particles fall and leave the computational domain, while other particles remain on the plate and form a heap as a result. For this test case a set of experiments was conducted and compared with the simulations. The corresponding shape parameters for each sort of the particles were found and are listed in Table 5. Corresponding superquadric approximations of the particles in question are presented in Fig. 12. Material properties were chosen the same for each particle shape and are listed in Table 6.   Young's modulus, particle and wall (Pa) 1.0 × 10 6 Poisson's ratio, particle and wall 0.3 Coefficient of friction, particle-particle 0.6 Coefficient of friction, particle-wall 0.4 Coefficient of restitution, particle and wall 0.2 Fig. 11 Piling problem, simulation snapshots at the final time step, superquadrics (left) and volume equivalent spheres (right) without rolling friction, 4680 particles From the pictures (Figs. 13, 14 and 15) it can be seen that there is only a qualitative agreement between the experiments and the simulations mainly due to relatively big size of the particles used in comparison to the cylinder diameter.

Static packing of cylinders
In this validation test we simulate a static packing of cylinders, defined as superquadric particles. The particles were dropped into a cylindrical container, and compared to the experimental data provided by Kodam et al. [29]. The DEM material properties along with particle properties are listed in Table 7. Particle size parameters, a, b, and c, were chosen such that particles in the simulation have the same diameter and length as in the experiment: a = b = D/2,   (Table 7). Volume and principal moments of inertia of the cylinders defined as superquadrics have values smaller than those for true cylinders due to the rounded edges. The difference between a true cylinder and its superquadric approximation decreases with the increase of the blockiness/roundness superquadric shape parameter n 1 , however, leading to less stability of the method. Hence, a compromise value must be chosen. The superquadric cylinder in this simulation was set to have the same mass as the true cylinder by increasing density by 7 % and setting blockiness parameters to n 1 = 6.0, n 2 = 2.0. DEM time step was set t = 10 −5 s. The image of the final state from the simulation is shown in Fig. 16. The final experimental fill height according to Kodam et al. [29] is 53.3±2.0 mm, while superquadric DEM simulation gives the fill height of roughly 52.0 ± 3.0 mm which is in good agreement with the experiment. Kodam et al. [29] simulated packing of the cylinders with the multisphere approach and found that 9-sphere particles underpredict the bed height by 21 %, while 54-sphere particles underpredict the fill height by 8 %.

Hopper discharge
We also use the method in simulating the discharge of Smar-ties®chocolate candies from a flat bottom hopper. These candies have ellipsoidal shape with the following shape parameters: 2a = 13.56 mm, 2b = 2c = 7.19 mm, n 1 = n 2 = 2.0, and density 1338 kg/m 3 . The hopper is 290 mm along the X-direction, 55mm in Y-direction, and has an orifice of 54 mm in X-direction and 55mm in Y-direction. Gravity is oriented along the Z-direction. The particle-particle Fig. 18 Ratio of the particles remaining in the hopper as a function of time. Experimental data (Liu et al. [33]) vs. simulation and particle-wall coefficient of friction was set 0.4. Young's modulus, Poisson's ratio, and coefficient of restitution were chosen to be 10 GPa, 0.29, and 0.5 correspondingly. The hopper geometry and particle properties used in the DEM simulation are the same as those used by Liu et al. [33] and Dong et al. [16]. The DEM time step size was set t = 2 × 10 −5 s. A total of 5500 particles are dropped in the hopper to form a bed of 0.4 m height and remain motionless until the orifice is opened. After the settling is completed, the orifice located at the center of the bottom is opened and particles discharge from the hopper by gravity exhibiting Vshaped flow pattern. The simulation results are shown in Fig.  17 at different time steps during the discharge along with the snapshots observed in the experiment by Liu et al. [33]. Good agreement between the simulation and the experiment can also be found in terms of discharge rate, i.e., ratio of the remaining particles in the hopper at different times ( Fig.  18), which proves validity of the shape model. However, one can see that experimental results are consistently higher than the numerical results. The discrepancy between experimental and simulation values may be due to the coefficient of friction that must be calibrated to fit experimental data.

Hopper discharge-influence of the aspect ratio
In the next validation test case we apply the superquadric DEM code to modeling of the discharge of cylindrical particle from a flat bottom hopper (Fig. 19). The simulated particles have density ρ = 2500 kg/m 3 , fixed bottom diameter D = 1 mmf, while the aspect ratio α = h/D, i.e., the ratio between the height and the bottom diameter, was varied. Young's modulus, Poisson's ratio, and coefficient of restitution are the same as in the previous section. The hopper is 11D along X-direction, 4D along Y-direction. Periodic boundary conditions are applied in the Y-direction. The orifice is 3.6D along X-direction and 4D along Y-direction. A total of 2100 particles were simulated with aspect ratio α = 0.33, 1400 particles with α = 0.5, 700 particles with α = 1, and 467 particles with α = 1.5 correspondingly. As a consequence, the total volume of particle before the discharge was 514 mm 3 for all shapes. If the coefficient of friction is small (μ = 0.1), it can be seen (Fig. 20) that the volume of discharged particles as function of time is comparable for different shapes and is almost linear until the discharging is about to end [16], while the number of discharged particles, obviously, varies for different shapes. This phenomenon occurs in rheology of granular particles and can be a test case for validation of a non-spherical DEM code [32].

Conclusion
The superquadric shape model was implemented as a separate surface model in the open-source DEM package LIGGGHTS® [28] which is an extension of the open-source package LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [46]. Both are massively parallel and written in C++. The program codes are available for public download.
The superquadric DEM has shown promising results along with qualitative and quantitative agreement with experimental data reported in the literature. This paper shows versatility and applicability of the superquadric DEM. The methods for contact detection, which can easily take up to 80 % of computational time [41], and contact force calculation between superquadric particles have been described in detail by using their implicit equations. The formulation employs the "contact point" which is midway and closest to both particles. The corresponding algorithms for particle-mesh interaction have also been developed but are not presented in this paper. The superquadric DEM code has been applied to various DEM problems which prove robustness and efficiency of the implemented algorithms. The methods have shown to be fast, and can be further optimized.
The superquadric particles are expected to give more accurate results than multisphere approximations for more reasonable computational time. Detailed comparison of superquadrics and multispheres in DEM is to be done in the future. The proposed methodology has the potential to be further extended for any other type of particles defined by a potential/shape function. The code is expected to be available for public download in 2017.