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

Authors: Chengbo Li; Qianqian Gao; Xianchang Li
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

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. The contact detection of particles with strictly convex shape function, such as ellipses and superellipses, is solved by Newton–Raphson method, and a grid method is provided to deal with heart-shaped particles. The grid method can be generalized into a particle system, in which the shape function is not convex. The accuracy and stability of the algorithm are verified by a series of tests. For the collision of a pair of ellipses with an aspect ratio ofα = a/b = 1000, the efficiency 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. The 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.


---

## Body

## 1. Introduction

The granular matter is ubiquitous in nature. Various amorphous structures can be found in aggregates of nonspherical or noncircular particles [1]. These aggregate models are usually used to study the structures of many materials, such as, glass [2], nanoparticles [3], and colloids [4–6]. The discrete element method (DEM) can easily obtain the microscopic properties that are difficult to obtain in experiments, and has been an effective method to study the process of the particle system. Scientific attention has broadened the study of spherical aggregates, which are the simplest shapes in the past few decades [7–9]; moreover, the packings of nonspherical or noncircular particles have attracted considerable attention [10–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]. The 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–22], and characteristic polynomial methods [13–15, 23], have been introduced and developed [24].The main idea in PW overlap potential algorithm is described as follows: Two particles are fixed at their centroids. Then, 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. The two particles are in external tangency if the scaling factor is equal to 1. The 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 vectornr 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. The 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. Speciﬁcally, 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. This result allows us to compute the conﬁguration 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 modiﬁed sign variation number of the signed subresultant sequence. The algebraic conditions given by Wang et al. [23] laid the theoretical foundation for the follow-up practical applications in the collision detection for ellipsoids. The 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 round-off 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.This 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. The 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. The grid method is also suitable for ellipses and superellipses and is verified by simulations of random packing under gravity and hopper discharge. The proposed method never diverges without an initial iteration point as needed in the NR method. The 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.


## 2. 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)(1)fx=xa2ζ+yb2ζ1/ζ=α,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 deﬁne the particle shape, even though it is not strictly necessary. It is a simple ellipse when ζ=1 and fx=1. It is a rectangle with sides 2a and 2b when ζ⟶∞ and fx=1. Another definition of superellipse is based on(2)fx=xan+ybm,where n and m are positive exponents, and fx=C defines different equipotential surfaces. The 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.Figure 1
Superellipses of various equipotential surfaces.The equipotential surface value ofC controls the size of a particle. That is, a larger value of C corresponds to a larger particle when the shape function is the same. The 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. The parameters of three particles, namely,i, j, and k, are listed in Table 1, in which θ is the orientation, xc and yc are the position coordinates of the centroid, vx and vy are the translational velocities, and ω is rotational velocity.Table 1
Parameters of superellipses.
Particlea/mmb/mmnmθ/radxc/mmyc/mmvx/mm·s−1vy/mm·s−1ω/rad·s−1i2.520262π/35−913020j31046π/471520−501k31046π/4−8−5———Particlesi 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. Therefore, 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. They are separated if the scaling factor is higher than 1 and overlapped if the scaling factor is smaller than 1. Thus, the problem of contact detection is translated into the solution of the scaling factor. The solving process is presented in detail as follows.Figure 2
Two kinds of relative positions of superellipses. (a) Separated superellipses. (b) Overlapped superellipses.
(a)(b)The shape functions of particlesi and j in the body-fixed coordinate frame are(3)fix=xaini+ybimi,fjx=xajnj+ybjmj.Their shape functions in the space-fixed coordinate frame can be transformed to(4)FiX=fiRiX−Xi,FjX=fjRjX−Xj,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, β=iorj. If the orientation of particle β is θβ, then the expression of the rotational matrix is(5)Rβ=cosθβsinθβ−sinθβcosθβ.If the two particles are tangent, then only one pointXT satisfies fiRXT−Xi=fjRXT−Xj because this point is on both particles i and j. At this point, the gradients of both surfaces are parallel (with opposite direction), namely, ∇XFiX+∇XFjX=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,λFiX=1 means that particle i is enlarged to the equipotential surface of Ci=1/λ, and 1−λFjX=1 means that particle j is enlarged. If the rescaled particles are in external tangency, then the following equations exist:(6)λ∇XFiX+1−λ∇XFjX=0,a,FiX=FjX,b,Introducingα2=1−λ/λ and bringing the same equations obtained by Cleary et al. [31](7)∇XFiX+α2∇XFjX=0,a,FiX=FjX,b.Usually, these equations can be solved by the NR method.The NR method provides a very efficient way of converging to a multi-dimensional root. LetX denote the entire vector of values Xn and F denote the entire vector of functions Fn, n=1,...,N. In the neighborhood of X, each of the functions Fn can be expanded in Taylor series,(8)FnX+δX=FnX+∑m=1N∂Fn∂XmδXm+OδX2.The matrix of partial derivatives that appear in equation (8) is the Jacobian matrix J, which is described as(9)Jnm=∂Fn∂Xm.Then, equation (8) can be rewritten as follows:(10)FX+δX=FX+J⋅δX+OδX2.By neglecting terms of orderδX2 and higher and setting FX+δX=0, we obtain a set of linear equations for the corrections δX that move each function closer to zero simultaneously, namely,(11)J⋅δX=−FX.Equation (10) can be solved by the LU decomposition of the matrix [32]. The corrections are then added to the solution vector,(12)Xnew=Xold+δX,and the process is iterated to convergence.The aforementioned solution method is substituted into equation (7), FX=F1X,F2X,F3X and variables X=X,Y,α are obtained, where(13)F1X=∂FiX∂X+α2∂FjX∂X,F2X=∂FiX∂Y+α2∂FjX∂Y,F3X=FiX−FjX.After simple variable substitutions, the transition formulas of gradients between the space-fixed coordinate frame and body-fixed coordinate frame are expressed as follows:(14)∇XFβX=RβT∇xfβx,(15)x=RβX−Xβc,where RβT, Xβc are the transposed matrix of the rotational matrix Rβ and centroid of particle β, β=iorj.The Jacobian matrix of particleβ reads(16)Jβ=RβTHxfβxRβ,where Hxfβx=∂2fβ/∂x2∂2fβ/∂x∂y∂2fβ/∂y∂x∂2fβ/∂y2 is the Hessian matrix of the shape function fβx [33, 34].The Jacobian matrix of equation (11) is(17)J=J11i+α2J11jJ12i+α2J12j2α∂Fj∂XJ21i+α2J21jJ22i+α2J22j2α∂Fj∂Y∂Fi∂X−∂Fj∂X∂Fi∂Y−∂Fj∂Y0.Substitute equations (13)–(17) into equations (11) and (12), iteration stops if either the sum of the magnitudes of FiX is less than a tolerance, or the sum of the absolute values of the corrections δX is less than a tolerance.The particles in Figure2 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. The 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. The tangent points are located at XT1=14.6715,3.7305 and XT2=−2.8060,−11.7405 with a tolerance of 10−8, respectively (Figure 3). The average computational time is 19.4 ms.Figure 3
Contact detection of separated and overlapped particles. (a) Enlarged superellipses. (b) Reduced superellipses.
(a)(b)If particles move and rotate with the translational and angular velocities listed in Table1, the collision process is shown in Figure 4. The collision time is t=0.0197s, and the collision point is at XT=13.844,6.3381. It takes 0.265 s (440 time-steps) to calculate the process. The 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.Figure 4
Collision process of superellipsesi and j.NR method is a universal method used to determine the contact point of particles with strictly convex shape functions. The 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.


## 3. Contact Detection of Heart-Shaped Particles

The shape functions of heart-shaped particles are defined in different ways, such as equations (18) and (19)(18)fx=x2+y2−13−x2y3,(19)fx=x2+235x23−y2−1.The 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. Thus, the generalized shape function can be written as(20)fx,α=x2+y2−α3−x2y3=0.Particles with different values ofα are shown in Figure 5, and the shape functions are not strictly convex. The normal vector nx 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.Figure 5
Varied sizes of heart-shaped particles.A new algorithm, the so-called grid method, is introduced to address the contact detection of heart-shaped particles. The main idea of the grid method is simple. Similarly, a specific example is provided to illustrate this method. Three heart-shaped particles are denoted byih, jh, and kh 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, xc and yc are the position coordinates of the centroid, vx and vy are the translational velocities, and ω is rotational velocity.Table 2
Parameters of heart-shaped particles.
Particlesαθ/radxc/mmyc/mmvx/mm·s−1vy/mm·s−1ω/rad·s−1ih12π/92.555402jh2π/657−20−50−3kh2π/637———Figure 6
Two kinds of relative positions of heart-shaped particles. (a) Separated heart-shaped particles. (b) Overlapped heart-shaped particles.
(a)(b)Particlesih and jh are separated, but ih and kh are overlapped in Figure 6. Therefore, two scaling factors must be provided to guarantee that particle ih is in external tangency with particles jh and kh. The 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 jh is fixed and particle ih is rescaled to find the smallest α and the tangent point. First, one should calculate the points on the surface of particle jh with the centroid position Xj, orientation Rj and αj, that is, the following equation should be solved.(21)FjX,αj=fRjX−Xj,αj=0.Then, search the solutions of equation (21) to find the smallest value Fmin of(22)FiX,0=fRiX−Xi,0.Thus, the corresponding pointXmin is the contact point when particle ih is rescaled by αmin, and αmin satisfies fRiXmin−Xi,αmin=0.For convenience, the reason to determineXmin 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 particleih,(23)gx≡x2+y23−x2y3,denote A=x2+y2,B=x2y3, and substitute equation (23) into the shape function fx,α=A−α3−B=0, we obtain g=A3−A−α3, which means that the minimal g corresponds to the minimal α. A point xmin is in a series of points on the surface of particle jh which makes the value of g in equation (22) be the smallest gmin. In other words, when particle ih is separated from particle jh, particle ih is enlarged slowly, and the first common point xmin is on the equipotential surface of αmin, thereby meeting the condition gmin=A3−A−αmin3.Another way is to directly solve the equationfRiX−Xi,α=0 for the points on the surface of particle jh and finds xmin and αmin. However, computational efficiency and time have not been significantly improved in this way.The aforementioned process is used to deal with the example as shown in Figure7. For the case of separation, particle ih 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 ih is reduced to the surface of α=0.3674, and the contact point is located at X=2.5675,5.9296. The average computational time (approximately 19.2 ms) to detect the tangent point and the scaling factor is almost the same as that of superellipses.Figure 7
Contact detection of separated and overlapped heart-shaped particles. (a) Enlarged heart-shaped particles. (b) Reduced heart-shaped particles.
(a)(b)If particles move and rotate with the translational and angular velocities listed in Table2, then the collision process is shown in Figure 8. The collision time is t = 0.0088 s, and the collision point is at XT=3.1465,6.3711. It takes 0.173 s to calculate the collision process. The 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.Figure 8
Collision process of heart-shaped particles.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. The parameters of ellipses are listed in Table3. The contact detection is shown in Figure 9.Table 3
Parameters of ellipses.
Particlea/mmb/mmnmθ/radxc/mmyc/mmvx/mm·s−1vy/mm·s−1ω/rad·s−1i2.50.1222π/35−913020j1000.122π/471520−501k1000.122π/44−8———Figure 9
Contact detection of ellipses by using the grid method. (a) Separated ellipses. (b) Overlapped ellipses. (c) Enlarged ellipses. (d) Reduced ellipses.
(a)(b)(c)(d)Figures9(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).Figure 10
Collision process of two ellipses. The smaller ellipse is particlei and the bigger one is particle j. (a) Initial positions. (b) Final positions. (c) Process of collision. (d) Local zoom.
(a)(b)(c)(d)


## 4. Contact Detection Used in DEM

The grid method is used in the following superellipse system to further verify its validity. The contact forces and torques of simulation and parameters are listed in Tables4 and 5, more details can be found in previous studies [7, 35]. The particle system contains 8100 superellipses with shape function f1x≡x/1.010+y/1.010=1. The 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. The rotational and translational energies are zero. Then, the particles fall under gravity onto the bottom, which consists of particles with shape function f2x≡x/1.516+y/1.516=1. The periodic boundary condition is applied in the x direction. The initial and final states are shown in Figure 11.Table 4
Forces and torques between particles.
PropertiesEquationsNormal overlappingδnTangential displacement vectorδ∧t=δ⇀t/δ⇀tTangential maximum displacementδt,max=μs2−v/2−2vδnPoisson’s ratiovEffective Young’s modulusE∗=E/2−2v2Effective massmij=mimj/mi+mjEffective curvatureR∗=RiRj/Ri+RjNormal spring forcef⇀cn=4/3E∗R∗δn3/2n∧Normal viscous forcef⇀dn=−cn8mijE∗R∗δn1/2v⇀n,ijTangential spring forcef⇀ct=−μsf⇀cn1−1−δt/δt,max3/2δ∧tTangential viscous forcef⇀dt=−ct6μsmijf⇀cn1−δt/δt,max/δt,max1/2v⇀t,ijTorqueq⇀i=r⇀ic×f⇀ct+f⇀dtRelative sliding velocityv⇀ij=v⇀i−v⇀j+ω⇀i×r⇀ic−ω⇀j×r⇀jcNormal relative velocityv⇀n,ij=v⇀ij∙n∧n∧Tangential relative velocityv⇀t,ij=v⇀ij−v⇀n,ijRelative rolling velocityv⟶r=−ricnrjcn/rijnn∧×ω⟶i−ω⟶jRolling forcef⟶r=minμrfn,f⟶0rf∧0rf⟶0r=−krξ⟶r−γrv⟶rRolling torqueq⇀rolling=ricnrjcn/rijnn∧×f⟶rTable 5
Parameters of simulations.
ParametersValuesSI unitsLength-unit11 mmMass-unit11 mgTime-unit11 msParticle densityρ55 × 103 kg/m3Young’s modulusE1061 GPaFriction coefficientμs0.500.50Poisson’s ratiov0.250.25Normal viscosity coefficientcn0.40.4Tangential viscosity coefficientct0.40.4Rolling stiffnesskr106106 N/mRolling viscosityγr1031.0 kg/sRolling friction coefficientμr0.10.1Integration time-step1 × 10−31 × 10−6 sFigure 11
Random packing process of superellipses. (a) Initial state and (b) equilibrium state of the system.
(a)(b)The energy changes with time are shown in Figure12. 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. Thus, the final state is a stable equilibrium state.Figure 12
Energies change with time. (a) Translational energy, (b) rotational energy, (c) gravity potential energy, and (d) total energy.
(a)(b)(c)(d)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. Figure13 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. The stable arch is shown in Figure 6, and the results are consistent with that of a previous study [36].Figure 13
Snapshots of particle flow at different times. (a)t = 0 s. (b) t = 0.11 s. (c) t = 0.32 s.
(a)(b)(c)The contact detection proposed in this work can also be used to study the granular flow in a hopper similar to Reference [37]. The granular flows of the elliptical particle system shown in Figure 14 are simulated with the grid method under different outlet sizes. The relationships between the discharge rate and the outlet size are calculated and according to Beverloo’s law [38], as shown in Figure 15.Figure 14
Snapshot of granular flow in the hopper.Figure 15
Mass flow rate versus the outlet sizes for particles with aspect ratios about 3 : 1, the outlet sizes range from6dmaj to 16dmaj. Symbol: simulated data, line: fitted by using the Beverloo’s law.The two-dimensional hopper (Figure14) consists of a rectangular channel of height (H=800mm) and width (W=80mm). The outlet is situated at the center of the bottom and has a length (D), which is varied at will during the simulation process. The major principal axis and minor principal axis of the elliptical particles are dmaj=2mm and dmin=0.66mm. The 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. The 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.


## 5. Discussion and Conclusion

The 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. Therefore, 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, particlesi and j are in Figures 16(a) and 16(b) have the same shape function fx=x2/0.12+y2/1002=1. Their centroid coordinates are Xi=1.2,3.4 and Xj=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. This fails correct solution because the eccentricity is too large. But by using the grid method and enlarging particle j to satisfy fjx=x2/0.12+y2/1002=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).Figure 16
Contact detection between two ellipses with large eccentricity. (a) Initial size. (b) Local zoom of the initial size. (c) Contact after enlargement. (d) Local zoom of contact.
(a)(b)(c)(d)The grid method is not only suitable for 2D cases but also for 3D superellipsoids. In Figure17, two superellipsoids i and j have different shape functions of fix=x/46+y/46+z/412=1 and fjx=x/28+y/38+z/28=1, their centroid coordinates are Xi=−10,−15,−10 and Xj=−20,−25,−20, and the corresponding Euler angles are ψi,θi,φi=π/9,π/8,π/7 and ψj,θj,φj=π/6,π/5,π/4, respectively. The initial velocities and angular velocities are vix,viy,viz=−20,−10,−15, vjx,vjy,vjz=10,20,15, ωix,ωiy,ωiz=20,−30,2, and ωjx,ωjy,ωjz=−10,20,−2, respectively. By running the program in the MATLAB compiling environment on a PC with a CPU of Intel(R) Core (TM), I9-8950HK 2.9 GHz, and 32 GB of RAM, the NR method is employed to simulate the collision process of these two particles as shown in Figure 17.Figure 17
Collision process of two superellipsoids. (a)t = 0 s. (b) t = 0.02 ms. (c) t = 0.04 ms. (d) t = 0.06 ms. (e) t = 0.08 ms. (f) t = 0.10 ms. (g) t = 0.12 ms. (h) t = 0.1924 ms.
(a)(b)(c)(d)(e)(f)(g)(h)In the simulation, the time-, length-, and mass-unit values are ms, mm, and mg, respectively. If the time-step is set as1×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. The 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. The advantage of the method is that the inverse of the Jacobian matrix and the initial iteration point are unnecessary. The grid method will never diverge, and the problem wherein singular matrices appear in the NR method is avoided.



---

*Source: 2514023-2022-08-31.xml*