Introduction

\label{introduction}
Particle laden flows are common phenomena appearing in natural and industrial systems: fog drops, blood, sediments, spray dryers, and fuel injectors are just a few of their extensive applications. While numerical study of these flows is challenging due to fluid-structure interaction of several particles with the carrier fluid, the interaction of particles with each other and with the enclosing surfaces adds to the complexity of the simulations. The current collision models mostly include definition of case dependent parameters that can affect the results. Moreover, the modeling becomes significantly more complicated when the particles have non-spherical shapes. This includes particles with arbitrary shapes such as sediments and dust storms or particles with non-spherical standard shapes such as biconcave shape of red blood cells, or the cylindrical shape of fibers. Thus, this study presents a Simplified Spring Collision Model (SSCM) for the contact treatment of spherical and non-spherical particles. It shows an inexpensive and robust model with parameters easily obtainable from the flow and particle motion.
The literatures studying various collision models span very diverse approaches [1-5]. The direct numerical simulation of collision of finite size particles as immersed objects mostly utilizes the discrete element method (DEM) [6] in conjunction with an immersed boundary method [7]. DEM considers the flow in an Eulerian framework and tracks the particles individually in the Lagrangian framework based on the Newton’s second law of motion. DEM soft sphere model (SSM) [8] is used for the cases with limited values of particle Stokes number in which the effect of flow on two nearby particles are not negligible i.e. when wet collisions occurs. SSM uses the differential format of particle dynamic motion and conservation laws update the particles’ velocities and positions in time [9, important parameters in the wet collision based on SSM.

1-1- Collision Forces

\label{collision-forces}
The simplest implementation of SSM obtains the collision forces by defining a repulsive potential for each particle. Glowinski et al [11] and also Wan & Turek [12] suggested that the collision force for spherical particles is proportional to the square of the amount of overlapping. The overlapping is a parameter representing how close the two particles are to each other. A similar approach uses a Lenard-Jones type potential to apply a repulsion force [13]. Even though the method is simple to implement, there exists several problems with it. Firstly, lack of physical reasoning for implementing the model can lead to unrealistic results [14]. While the method can prevent particles’ overlapping, it can enforce extra repulsion leading to unphysical results. Therefore, the accuracy of the results highly hinges on choosing case dependent parameters [15-17]. Secondly, the method ignores energy dissipation due to material damping [14]. Thirdly, definition of potential for non-spherical shape particles can be difficult and limited to classical shapes such as ellipsoids [18].
A common treatment in the particulate systems relies on modeling forces using some form of spring and potentially a type of damper. Several methods are suggested to find the spring system characteristics. One of the common methods defines the coefficients based on the particle’s solid structure characteristics. For the spherical particles Hertz theory defines the spring force as \(k_{\text{Hz}}\delta^{\frac{3}{2}}\) where \(\delta\) is the displacement calculated from the amount of overlapping [8, \(k_{\text{Hz}}\) is the Hertz spring constant obtained from the two particles material and geometrical properties. The Hertz theory as an elastic collision model agrees well with the experiments in dry collisions [8, but in a wet collision it neglects the viscous dissipation of energy. A dashpot can enforce the dissipation of energy during the collision to coincide the physical behavior. The damping coefficient of the dashpot is suggested to be proportional to some exponent of displacement [23, includes the solid object characteristics in the calculations, the implementation of the method to model wet collision can be inefficient even for spherical particles. This is because the time scales related to the actual solid contact mechanics is usually several orders of magnitude smaller than the fluid due to very high spring coefficient values [25].
Decrease of spring coefficient avoids problem of extra diverse time scales as it extends the collision time to a scale comparable to the flow scale. Consequently, while the momentum exchange is preserved, the problem changes from high impulse in a short time to low impulse in a long time. This situation is very helpful as it does not change the problem in the large time scale in addition to the fact that the problem becomes practically solvable. However, the question becomes how to set the spring and dashpot variables to reach the actual physical behavior.
Kempe and Fröhlich [14] proposed the adaptive collision model (ACM) by stretching the collision time to ten times of the flow solver’s time step size. The model inspired by a Hertz-like force definition for the spring, represents a second order ODE as the spring and dashpot system. Using each particle’s rebound velocity from the solid dry coefficient of restitution, the ODE coefficients were obtained from an iterative solver. When the gap between two particles is small the viscous forces become significant. Therefore, the lubrication model of Ladd [26] empowered the model to account for the non-negligible viscous dissipation. By proposing an explicit formulation to find the ODE coefficient instead of an iterative approach, Ray et al [4] simplified the implementation of the ACM. Both approaches ignore the hydrodynamic forces on the particles when the particle is very close to the wall. To avoid the unphysical behaviors during the enduring contacts, Biegert et al [27] limited the ignorance of hydrodynamic forces by including a critical Stokes number below which the hydrodynamic forces should not be ignored.
The idea of preservation of momentum exchange not only can help to increase the minimum time scale, it can also justify the implementation of a linear spring. When the velocity of the bounced particle is important and not the way this velocity is achieved, a linear ODE such as harmonic damped oscillator can also represent the collision. This concept inspired Izard et al. [28] and Costa et al. [29] and to use the idea of adjusting the spring and dashpot based on the contact duration. As suggested by Costa et al. [29] obtaining the contact during from the flow time scale softens the collision and saves the computational cost significantly. A modified lubrication force can enforce the dissipation [29]. Applying the dissipation through the lubrication forces as the way used in [14, 28, Firstly, the implementation of lubrication forces requires resolved grid at the places near the contact place which reduces the efficiency of the simulation. Secondly, as shown in [27] the dissipation can depend on the lubrication gap thickness. Biegert et al. [27] could decrease the dependency on the lubrication gap by using higher order sub-stepping for the particle motion. To avoid the complications related to implementation of a lubrication model, this paper offers a collision model with simplified definition of collision parameters to make the study of collision physics easier. With no implementation of lubrication forces, the SSCM model leads to having less expensive simulations and accurate results. This model overcomes the difficulties and shortcomings of the available collision models e.g. iterative schemes, necessity of using case sensitive parameters, and disparities between the time scales.

1-2- Collision of Non-Spherical Particles

\label{collision-of-non-spherical-particles}
Study of particle collision becomes much more challenging when the particles’ shapes are not spherical. Difficulties with studying these particles collision stems from the fact that finding the related collision parameters i.e. contact point and the collision forces is a tedious task. While for the spherical particles the contact point is located on the intersection of the particles’ interface and the center to center connecting line, for non-spherical objects the contact point location depends also on the the more complex geometry. Besides, for non-spherical objects the contact normal vector mostly does not pass through the gravity center; therefore, a moment is induced to the particle. Moreover, the collision forces depend on the amount of overlapping and the radii of curvature which are also more difficult to find when the particles do not have a spherical shape.
When the particles have a classical shape e.g. ellipsoids and polygons, the geometrical features help finding the contact/collision point. There are several ways to find the collision point of the non-spherical particles which mostly require an iterative approach to find the contact point [30-37]. Unfortunately, there is no common way to define the collision point for non-spherical particles. The collision point can be defined as the geometric center of the intersection volume [38]. Lin and Ng [30, line connecting the surface of the two particles as the collision point. The connection line can connect the deepest part of each particle into the other one’s similar point [39, common normals of the two particles are located [30]. The collision point of the ellipsoids also could be obtained from the closest distance between two particles as described by [32, approach sets up two balls inside each particle and tangent to the surface and updates the positions of the balls by sliding the balls to the intersection of the ellipsoids and the line connecting the ball centers. For the polygonal particles, Nezami et al. [33] used the common plane [34] between the particles to detect the collision point. Other contact point detection methods include polyhedrons [35, [37, [43] as well as “potential particles” [44, particles by defining a potential field which can morph from a sphere to cube. The resulted shape can have a shape of sphere, cube or a transition of sphere to cube. The method was extended to other classical shapes such as tetrahedral-shaped, prism and rounded tetrahedral-shapes.
For arbitrary particles, no special feature of the geometry can help to find the collision point. However, it is common to compose the particle from small geometric elements to form the particle. Thus, the original particle with arbitrary shape is a composite of the elements that are glued. The geometric elements can be spheres [46-53], cylinders [54] , ellipsoids [55] or a combination of several shapes [56]. The problem with combining different shapes is that the accuracy depends on the number of used elements. While fewer number of elements are easier and faster to form the final clumped geometry, the resulted surface will be rough and inaccurate and multiple collision points might be captured. To avoid this problem, recently use of “spherosimplices” for complex objects has gained popularity. Spherosimplices constructs the shapes by moving a cylinder or sphere over a frame named skeleton [57-61]. The method is shown to be faster than the composed particle method. However, the generation of the shape becomes complicated for the complex objects especially in 3D as it depends on how accurate the shape of skeleton is set up. Moreover, the sharp edges will be smoothed because of using a cylinder or sphere to generate a shape.
The other way to find the contact point of the arbitrary shape particles is to classify the simulation grid points into fluid and each individual particle. Wachs [62] implemented this concept by using a spatial algorithm to define the grid points in the context of fictitious domain method [63]. We believe classification of grid point is the most general method in definition of arbitrary shape particles as no analytical representation will restrict the definition of the particles’ shapes. Moreover, finding the contact point which is a challenging task for non-spherical particles becomes easy. The method presented by Wachs [62] suggests a range for the spring/dashpot stiffness based on the problem and an initial estimation is required for the for their values. Moreover, the direction of normal forces needs to be defined clearly.
To overcome the shortcoming of the current collision models, this paper proposes a soft collision model for the simulation of arbitrary shape particles by using a field variable to find the overlapping places. This helps us to avoid many difficulties related to finding the collision points, collision normals, local particle interface curvature and the amount of overlapping. Moreover, as a non-iterative scheme, the method saves time in finding the collision parameters. We also propose a novel spring method for the collision which facilitates implementation of the model significantly. A particle tagging system is also presented to conduct the model efficiently for high particle numbers. The method does not include any case dependent parameters and can handle any particle shape. The results show that the model is able to simulate particle-particle and particle-wall collision with proper accuracy.