A Study of Contact Detection between Noncircular Particles in Discrete Element Method

One of the key points of modeling noncircular particles in the discrete element method is the contact detection of particles. In this study, a general contact detection algorithm of two-dimensional particles with analytic shape functions is provided. e contact detection of particles with strictly convex shape function, such as ellipses and superellipses, is solved by Newton–Raphsonmethod, and a grid method is provided to deal with heart-shaped particles. e grid method can be generalized into a particle system, in which the shape function is not convex.e accuracy and stability of the algorithm are veried by a series of tests. For the collision of a pair of ellipses with an aspect ratio of α a/b 1000, the eciency is not worse than the Newton–Raphson method. For random packing under gravity, no residual kinetic energy is observed, and the force that acts on the bottom is equal to the gravity, that is, the system reaches a mechanical equilibrium state. After equilibrium, the process of hopper discharge is also simulated.e present method is suitable for arbitrarily shaped particles with analytic shape functions in two-dimensional cases. In addition to superellipses, the method can also be applied to other particle systems as long as the shape functions are given.


Introduction
e granular matter is ubiquitous in nature. Various amorphous structures can be found in aggregates of nonspherical or noncircular particles [1]. ese aggregate models are usually used to study the structures of many materials, such as, glass [2], nanoparticles [3], and colloids [4][5][6]. e discrete element method (DEM) can easily obtain the microscopic properties that are di cult to obtain in experiments, and has been an e ective method to study the process of the particle system. Scienti c attention has broadened the study of spherical aggregates, which are the simplest shapes in the past few decades [7][8][9]; moreover, the packings of nonspherical or noncircular particles have attracted considerable attention [10][11][12].
In DEM, contact detection is the most time-consuming part. Moreover, the computational time is notably enhanced when determining the interaction of nonspherical or noncircular particles whose mathematically rigorous treatment is notably nontrivial [13,14]. e simplest noncircular and nonspherical particles are ellipses and ellipsoids, whose surfaces are described by a quadratic equation. If two variables are found in the polynomial, then we have a conic (curve) and a quadric surface for three variables [15].
Several contact detection algorithms, such as intersection strategies [16], curvature simpli cations [17], Perram and Wertheim (PW) overlap potential [18][19][20][21][22], and characteristic polynomial methods [13][14][15]23], have been introduced and developed [24]. e main idea in PW overlap potential algorithm is described as follows: Two particles are xed at their centroids. en, one or both of them are rescaled so that they have only one common point. Next, the scaling factor and the contact point are obtained. e two particles are in external tangency if the scaling factor is equal to 1. e particles are separated if the scaling factor is higher than 1 and overlapped if the factor is smaller than 1. A requirement of particle shape being smooth and convex is important in this algorithm because it makes the normal vector n(r) be a one-to-one function. A unique normal vector exists at any point on the particle surface, and one can determine the point on the surface that corresponds to normal vector as long as it is provided. e overlap potential can be determined analytically for ellipsoids, and the maximum of the overlap potential can be easily found numerically [25].
In characteristic polynomial methods, the relative position of two ellipsoids is related to the root pattern of their characteristic equation. Specifically, two ellipsoids are separated if and only if their characteristic equation has two distinct positive roots. Furthermore, they touched each other externally if and only if one positive double root exists. is result allows us to compute the configuration of two ellipsoids by simply determining the number of the positive roots of a polynomial, which can be derived by the classical Descartes' rule and the modified sign variation number of the signed subresultant sequence. e algebraic conditions given by Wang et al. [23] laid the theoretical foundation for the follow-up practical applications in the collision detection for ellipsoids.
e solution of the quartic polynomial equation for two-dimensional (2D) elliptical particles is complicated and unstable because of considerable round-off errors caused by the limited computer accuracy. Details of the corresponding analysis have been described in Reference [26]. Similarly, the solution of the sixth-order polynomial equation for three-dimensional (3D) ellipsoidal particles is more complicated and unreliable because of several roundoff errors [27]. To avoid directly solving a quartic polynomial equation or a sixth-order polynomial equation, a parametric equation of ellipsoid is utilized to transform the sixth-order polynomial equation into a quartic polynomial equation by two approximations [28]. Furthermore, a golden section search algorithm is used instead of a sixth-order polynomial equation to detect the overlapping of ellipsoids [29]. However, the algorithm is restricted to a narrow range of aspect ratios.
is study attempts to present a universal method to determine the contact point of 2D particles that have analytical shape functions for particles with strictly convex shape functions, such as ellipses and superellipses. e Newton-Raphson (NR) method is used to obtain the scaling factor and the contact point, and a grid method is used to deal with heart-shaped particles with concave shape function. e grid method is also suitable for ellipses and superellipses and is verified by simulations of random packing under gravity and hopper discharge. e proposed method never diverges without an initial iteration point as needed in the NR method. e present method is suitable for arbitrarily shaped particles with analytic shape functions in 2D cases. In addition to superellipses, the method can also be applied to other particle systems as long as the shape functions are given.

Contact Detection of Superellipses
Compared with ellipses, studies on superellipse are infrequent [21]. For example, a superellipse is given by the shape function (in the body-fixed coordinate frame) where ζ ≥ 1 is an exponent, a and b are the semiaxis lengths, and α is the value of equipotential surface. An exponent 1/ζ is added above to properly normalize the convex function and define the particle shape, even though it is not strictly necessary. It is a simple ellipse when ζ � 1 and f(x) � 1. It is a rectangle with sides 2a and 2b when ζ ⟶ ∞ and f(x) � 1. Another definition of superellipse is based on where n and m are positive exponents, and f(x) � C defines different equipotential surfaces. e value of C determines the size of a superellipse. Figure 1 shows several superellipses with various equipotential surfaces at the same centroid and orientation (θ � π/4), where the shape parameters are n � 6, m � 8, a � 2, and b � 4. e equipotential surface value of C controls the size of a particle. at is, a larger value of C corresponds to a larger particle when the shape function is the same. e main idea of the PW overlap potential algorithm is that two particles can be rescaled to contact each other if they are separated or overlapped and to obtain the scaling factor (or the value of C) and the contact point position. Over many repetitions, two particles are in external tangency if the scaling factor is equal to 1.
Here, a specific example is provided to explain the idea. e parameters of three particles, namely, i, j, and k, are listed in Table 1, in which θ is the orientation, x c and y c are the position coordinates of the centroid, v x and v y are the translational velocities, and ω is rotational velocity.
Particles i and j are separated, as shown in Figure 2(a). A scaling factor is used to enlarge the particles simultaneously. If particles are overlapped, as shown in Figure 2(b), they can be reduced to be in external tangency by using the scaling factor. erefore, whatever particles are separated or overlapped, only one scaling factor guarantees that they are in external tangency. Two particles are in external tangency if the scaling factor is equal to 1. ey are separated if the scaling factor is higher than 1 and overlapped if the scaling factor is smaller than 1.
us, the problem of contact detection is translated into the solution of the scaling factor. e solving process is presented in detail as follows. e shape functions of particles i and j in the body-fixed coordinate frame are eir shape functions in the space-fixed coordinate frame can be transformed to where X β is the centroid of particle β in the space-fixed coordinate frame, and R β is the rotational matrix between the body-fixed coordinate frame and the space-fixed coordinate frame, β � i or j. If the orientation of particle β is θ β , then the expression of the rotational matrix is

Advances in Materials Science and Engineering
If the two particles are tangent, then only one point ) because this point is on both particles i and j. At this point, the gradients of both surfaces are parallel (with opposite direction), namely, ∇ X F i (X) + ∇ X F j (X) � 0. If the two particles are separated, as shown in Figure 2(a), a scaling factor λ (or C) that is used to rescale the particles and guarantee one unique contact point between particles must be available [19,30]. Given λ ∈ (0, 1),λF i (X) � 1 means that particle i is enlarged to the equipotential surface of C i � 1/λ, and (1 − λ)F j (X) � 1 means that particle j is enlarged. If the rescaled particles are in external tangency, then the following equations exist: Introducing α 2 � (1 − λ)/λ and bringing the same equations obtained by Cleary et al. [31] Usually, these equations can be solved by the NR method.
e NR method provides a very efficient way of converging to a multi-dimensional root. Let X denote the entire vector of values X n and F denote the entire vector of functions F n , n � 1, ..., N. In the neighborhood of X, each of the functions F n can be expanded in Taylor series, e matrix of partial derivatives that appear in equation (8) is the Jacobian matrix J, which is described as en, equation (8) can be rewritten as follows: By neglecting terms of order δX 2 and higher and setting F(X + δX) � 0, we obtain a set of linear equations for the corrections δX that move each function closer to zero simultaneously, namely, Equation (10) can be solved by the LU decomposition of the matrix [32]. e corrections are then added to the solution vector, and the process is iterated to convergence. e aforementioned solution method is substituted into equation (7), After simple variable substitutions, the transition formulas of gradients between the space-fixed coordinate frame and body-fixed coordinate frame are expressed as follows: x where R T β , X βc are the transposed matrix of the rotational matrix R β and centroid of particle β, β � i or j. e Jacobian matrix of particle β reads where [33,34]. e Jacobian matrix of equation (11) is Substitute equations (13)-(17) into equations (11) and (12), iteration stops if either the sum of the magnitudes of F i (X) is less than a tolerance, or the sum of the absolute values of the corrections δX is less than a tolerance. e particles in Figure 2 are rescaled by using the NR method to find the external tangent point. Obviously, particles i and j need to be enlarged, and particles i and k should be reduced. e value of the equipotential surface is scaled from 1 to 6.2893 for particles i and j and to 0.3789 for particles i and k. e tangent points are located at X T1 � (14.6715, 3.7305) and X T2 � (− 2.8060, − 11.7405) with a tolerance of 10 − 8 , respectively (Figure 3). e average computational time is 19.4 ms.
If particles move and rotate with the translational and angular velocities listed in Table 1, the collision process is shown in Figure 4. e collision time is t � 0.0197 s, and the collision point is at X T � (13.844, 6.3381). It takes 0.265 s (440 time-steps) to calculate the process. e blue curves illustrate the initial positions of particles, the green ones illustrate the trajectories, and the red ones illustrate the positions at the collision moment.
NR method is a universal method used to determine the contact point of particles with strictly convex shape functions. e solution converges steadily to the unique contact point through a scaling factor if the initial iteration point is provided properly. If the shape function is not strictly convex and the solution is not unique, then it is likely to cause the divergence of the NR method. For example, it will       fail to determine the contact point between heart-shaped particles.

Contact Detection of Heart-Shaped Particles
e shape functions of heart-shaped particles are defined in different ways, such as equations (18) and (19) e size of heart-shaped particles is controlled by a subtracted constant (here, the subtracted constant is "1") in equation (18) instead of the value of the shape function.
us, the generalized shape function can be written as Particles with different values of α are shown in Figure 5, and the shape functions are not strictly convex. e normal vector n(x) of each point on the particle is not a one-to-one function, which may be the reason for the failure of using the NR method.

Advances in Materials Science and Engineering
A new algorithm, the so-called grid method, is introduced to address the contact detection of heart-shaped particles.
e main idea of the grid method is simple. Similarly, a specific example is provided to illustrate this method. ree heart-shaped particles are denoted by i h , j h , and k h with the parameters as listed in Table 2, and their relative positions and orientations are shown in Figure 6. In Table 2, θ is the orientation, x c and y c are the position coordinates of the centroid, v x and v y are the translational velocities, and ω is rotational velocity.
Particles i h and j h are separated, but i h and k h are overlapped in Figure 6. erefore, two scaling factors must be provided to guarantee that particle i h is in external tangency with particles j h and k h . e key problem is to determine the smallest value of the scaling factor that is used to rescale a particle to touch with another one. Generally, particle j h is fixed and particle i h is rescaled to find the smallest α and the tangent point. First, one should calculate the points on the surface of particle j h with the centroid position X j , orientation R j and α j , that is, the following equation should be solved.
en, search the solutions of equation (21) to find the smallest value F min of us, the corresponding point X min is the contact point when particle i h is rescaled by α min , and α min satisfies For convenience, the reason to determine X min and α min in the body-fixed coordinate frame is only a rotation in the space-fixed coordinate frame with the relative position unchanged.
Given a function in the body-fixed coordinate frame of particle i h , denote A � x 2 + y 2 ,B � x 2 y 3 , and substitute equation (23) into the shape function  Advances in Materials Science and Engineering 3 , which means that the minimal g corresponds to the minimal α. A point x min is in a series of points on the surface of particle j h which makes the value of g in equation (22) be the smallest g min . In other words, when particle i h is separated from particle j h , particle i h is enlarged slowly, and the first common point x min is on the equipotential surface of α min , thereby meeting the condition Another way is to directly solve the equation f(R i (X − X i ), α) � 0 for the points on the surface of particle j h and finds x min and α min . However, computational efficiency and time have not been significantly improved in this way. e aforementioned process is used to deal with the example as shown in Figure 7. For the case of separation, particle i h is enlarged to the surface of α � 2.5532 and the contact point is located at X � (3.3989, 6.6597). For the case of overlapping, particle i h is reduced to the surface of α � 0.3674, and the contact point is located at X � (2.5675, 5.9296). e average computational time (approximately 19.2 ms) to detect the tangent point and the scaling factor is almost the same as that of superellipses.
If particles move and rotate with the translational and angular velocities listed in Table 2, then the collision process is shown in Figure 8. e collision time is t = 0.0088 s, and  When the grid method is used to determine the contact point of superellipses, the computational time is 18 ms and slightly lower than that of using the NR method. Nevertheless, the grid method is suitable for ellipses with a very large aspect ratio of 1000 and a size dispersion of 40 times, in which case the NR method will fail. e parameters of ellipses are listed in Table 3. e contact detection is shown in Figure 9.
Figures 9(a) and 9(c) illustrate the contact detection between particles i and j by enlarging particle i, and Figures 9(b) and 9(d) illustrate the contact detection between particles i and k by reducing particle i. If particles i and j move and rotate with the translational and angular velocities listed in Table 3, then the collision process is shown in Figures 10. Figures 10(a) and 10(b) are the initial and final positions of the two particles, respectively. In Figure 10(c), the green curves illustrate the initial positions of particles, the blue ones illustrate the trajectories, and the red ones illustrate the positions at the collision moment. A local zoom of the collision process is shown in Figure 10(d).

Contact Detection Used in DEM
e grid method is used in the following superellipse system to further verify its validity. e contact forces and torques of simulation and parameters are listed in Tables 4 and 5, more details can be found in previous studies [7,35]. e particle system contains 8100 superellipses with shape function f 1 (x) ≡ |x/1.0| 10 + |y/1.0| 10 � 1. e initial body-fixed coordinate frame coincides with the space-fixed coordinate frame. In the beginning, particles are packed cubically, and no overlapping occurs between them. e rotational and translational energies are zero. en, the particles fall under gravity onto the bottom, which consists of particles with shape function f 2 (x) ≡ |x/1.5| 16 + |y/1.5| 16 � 1. e periodic boundary condition is applied in the x direction. e initial and final states are shown in Figure 11.   Figure 12. At the final equilibrium state, the kinetic energy is almost 0, and the force that acts on the bottom is equal to the gravity of the system. us, the final state is a stable equilibrium state.
After reaching the equilibrium state, a hole is opened at the bottom, and its width is 24 mm, which is 12 times the width of the particle. Figure 13 shows the snapshots of particle flow at different times. After 0.11 s, 67 particles drop through the hole, a semicircular structure of force arch comes into being, and the particle flow stops. e stable arch is shown in Figure 6, and the results are consistent with that of a previous study [36]. e contact detection proposed in this work can also be used to study the granular flow in a hopper similar to Reference [37]. e granular flows of the elliptical particle system shown in Figure 14 are simulated with the grid method under different outlet sizes. e relationships between the discharge rate and the outlet size are calculated and according to Beverloo's law [38], as shown in Figure 15. e two-dimensional hopper (Figure 14) consists of a rectangular channel of height (H � 800 mm) and width (W � 80 mm). e outlet is situated at the center of the bottom and has a length (D), which is varied at will during the simulation process. e major principal axis and minor principal axis of the elliptical particles are d maj � 2 mm and d min � 0.66 mm. e particles are initially randomly placed in the closed hopper, and then they settle under the gravity to minimize the gravitational potential and are allowed to leave the hopper from the outlet centered at the bottom. During the steady discharge process, the number of elliptical particles as a function of time is recorded every 200 time-steps, from which we measure globally the instantaneous flow rate as ΔN/Δt, where ΔN is the number of particles passing through the exit during the time interval Δt. e mean flow rate Q is obtained by temporal averaging the instantaneous flow rate during the steady flow state. Based on the above calculation, the dependence of discharge rate on the outlet size is examined, as shown in Figure 15.

Discussion and Conclusion
e iterative convergence direction and step length are difficult to control when solving the contact point between ellipses with large eccentricity by using the NR method, and singular matrices often appear when solving the inverse Jacobian matrix, which makes the solution fail. erefore, it is difficult to use the NR method in discrete element simulation of the particle system. However, solving the inverse matrix is not needed in the grid method which is not affected by the major and minor axes of particles. For instance, particles i and j are in Figures 16(a) and 16(b) have the same shape function f(x) � x 2 /0.1 2 + y 2 /100 2 � 1. eir centroid coordinates are X i � (1.2, 3.4) and X j � (2.1, − 19.5) and angular coordinates are θ i � π/5.6 and θ j � π/4.5, respectively. When using the NR method to calculate the contact point, singular matrices often appear when solving the inverse Jacobian matrix. is fails correct solution because the eccentricity is too large. But by using the grid method and enlarging particle j to satisfy f j (x) � x 2 /0.1 2 + y 2 /100 2 � 13.6177 in the contact detection, the contact point is determined at (54.4032, − 81.2724) as shown in Figures 16(c) and 16(d).
In the simulation, the time-, length-, and mass-unit values are ms, mm, and mg, respectively. If the time-step is set as 1 × 10 − 7 s, the collision contact point will be located at (− 16.5993, − 20.3341, − 14.2559) after 1924 time-steps which means that there are 1924 contact detections of arbitrary positions and angular coordinates. e calculation process of the grid method takes 0.986460 s, and the time of the NR method is 1.865421 s. Under the same solution accuracy ( < 10 − 3 ), the time consumption of grid method is less than that of NR method.
In addition to dealing with superellipses, the use of grid method in detecting the contact of heart-shaped particles is a universal approach and is suitable for all particles with analytical shape functions. e advantage of the method is that the inverse of the Jacobian matrix and the initial iteration point are unnecessary. e grid method will never diverge, and the problem wherein singular matrices appear in the NR method is avoided.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.