2- Methods

\label{methods}
The remedy to implementation of very small time scales is engaged with one interesting feature of the wet collisions: In a wet collision the material properties do not act as the dominant parameter because the collision depends mostly on flow properties represented by the Stokes number [64, \(St=\frac{\rho_{p}}{9\rho_{f}}\text{Re}_{p}\), where \(\rho_{p}\),\(\ \rho_{f}\) and \(\text{Re}_{p}\) are particle and fluid densities and the particle’s relative Reynolds number. \(St\cong 10\) is the critical Stokes number below which no rebound happens [64, which is compatible with the experiments [67]. When \(St>500\), the effect of fluid becomes negligible and dry collision models become applicable [66]. Having the correlation between the restitution coefficient, \(e=\frac{u_{\text{out}}}{u_{\text{in}}}\), and the Stokes number established, the question now is how the spring-dashpot constants should be defined for a general shape particles so that the rebound velocity could be calculated correctly.
The method starts with introducing a field variable i.e. overlapping function with values that vary like a potential based on the distance from the particle interface. The definition of the overlapping function depends on how the particles are defined. For arbitrary shapes the value of the particle level-set, the signed distance from the interface, can define the overlapping function. The overlapping function itself does not directly contribute to the collision forces; however, it helps with evaluating the forces. A cut off value for the overlapping function could be used because the far field grid points for each particle will not contribute to the calculations. One can also use a \(\tanh\) profile of the normalized distance variable to lift the necessity of using a cut off:
\(\psi_{i}(\mathbf{x},t)=\frac{1}{2}\left[\tanh{\left(-\frac{d_{i}\left(\mathbf{x},t\right)}{\xi_{i}}\right)+1}\right]\) (1)
where \(d_{i}\left(\mathbf{x},t\right)\) is the distance from the interface of particle \(i\) and \(\xi_{i}\) is the thickness of the interface that can be adjusted for each problem. We will discuss this variable later in section (3-3). \(\psi_{i}\) in Equation (1), varies from zero outside the particle to one inside the particle. On the interface \(\psi_{i}=0.5\) and the transition between solid to fluid depends on the interface thickness.
For arbitrary shapes, the value of the particle level-set, i.e. the signed distance from the interface, can define the overlapping function. For classical particle shapes the overlapping function could be defined in an easier way by using the analytical definition of the particle. A super-quadric equation [37, can represent a new distance function for a broad range of shapes as
\(d_{i}\left(\mathbf{x},t\right)=\left(\left(\frac{x}{a}\right)^{\frac{2}{\varepsilon_{2}}}+\left(\frac{y}{b}\right)^{\frac{2}{\varepsilon_{2}}}\right)^{\frac{\varepsilon_{2}}{\varepsilon_{1}}}+\left(\frac{z}{c}\right)^{\frac{2}{\varepsilon_{1}}}-1\) (2)
\(\varepsilon_{1}\ \)and \(\varepsilon_{2}\) are parameters that control the roundedness or blockiness of the shape. For an ellipsoid, \(\varepsilon_{1}=\varepsilon_{2}=1\) but for a rectangular prism \(\varepsilon_{1}=\varepsilon_{2}=0.2\).
The total overlapping function is the linear summation of each particle’s overlapping function:
\(\psi\left(\mathbf{x},t\right)=\sum_{i=1}^{N_{p}}\psi_{i}\) (3)
Here \(N_{P}\) is the total number of particles. As long as the particles are far apart, their overlapping function fields does not overlap and \(\psi_{i}=\psi\) on the related grid point. When two particles are close, \(\psi>\psi_{i}\ \)and the amount of overlapping signifies how much the particles are close. At each node near the vicinity of the particle interface, the degree of overlapping is
\(\text{dψ}_{\text{ij}}\left(\mathbf{x},t\right)=\psi-\psi_{i}\) (4)
where the index ij denotes the overlapping between particle \(i\) and particle \(j\). Closer to the particle the value of \(\text{dψ}_{i}\)will be higher. For two nearby particles, there is a region where the value of \(\text{dψ}_{i}\) is more than a minimum cut off \(\varepsilon\). The volume of this region is
\(V_{c,ij}=\int{H(\text{dψ}_{\text{ij}}-\varepsilon)}\text{dV}\) (5)
Here \(H\) is the step function. As we will discuss in section 3-2, \(\varepsilon\) is a small cut off value for the overlapping function and is used because the far field grid points of each particle do not contribute to the collision. The coordinates of the center of collision then is measured as
\(\mathbf{X}_{cp,ij}=\frac{\int{\mathbf{x}H(\text{dψ}_{\text{ij}}-\varepsilon)}\text{dV}}{V_{c,ij}}\) (6)
The algorithm proceeds by calculating the normal vector field using the gradient of the particle overlapping function.
\(\mathbf{n}_{i}(\mathbf{x},t)=\frac{\nabla\psi_{i}}{\left|\nabla\psi_{i}\right|}\) (7)
Integration of the normal vector field on the overlapping region leads to finding the collision normal vector and the unit normal vector:
\(\mathbf{N}_{cp,ij}\left(\mathbf{X}_{cp,ij},t\right)=\int\mathbf{n}_{i}H(\text{dψ}_{\text{ij}}-\varepsilon)dV\) (8)
\({\hat{\mathbf{N}}}_{cp,ij}\left(\mathbf{X}_{cp,ij},t\right)=\frac{\mathbf{N}_{cp,ij}}{\left\|\mathbf{N}_{cp,ij}\right\|}\) (9)
The volume of overlapping region also indicates the amount of overlapping or the overlapping distance as
\(\mathbf{\delta}_{N,ij}=\sqrt[3]{V_{c,ij}}{\hat{\mathbf{N}}}_{cp,ij}\) (10)
\(\mathbf{\delta}_{N,ij}\) is\(\ \)the overlapping distance of the particle \(i\) in proximity with the particle \(j\).
The next step is definition of a particle field velocity to determine the local velocity at each grid point related to the paritcle. For each particle, using the particle’s velocity and angular velocity the particle field velocity \({\tilde{\mathbf{u}}}_{i}\ \)is defined as
\({\tilde{\mathbf{u}}}_{i}\left(\mathbf{x},t\right)=\left(\mathbf{V}_{i}+\mathbf{\Omega}_{i}\times\left[\mathbf{x}-\mathbf{x}_{i}\left(t\right)\right]\right)\psi_{i}\) (11)
Then the relative particle field velocity is obtained:
\({\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{x},t\right)={\tilde{\mathbf{u}}}_{i}-\sum_{j\neq i}^{N_{P}}{\tilde{\mathbf{u}}}_{j}\) (12)
This definition becomes computationally intensive when the number of particles increases because it requires \(N_{P}(N_{P}-1)\) operations. When the particle number density is low and each particle experiences only one collision at time, the concept of total particle field velocity can be used to reduce the number of operations:
\(\tilde{\mathbf{u}}\left(\mathbf{x},t\right)=\sum_{i=1}^{N_{p}}{\tilde{\mathbf{u}}}_{i}\) (13)
Therefore, the relative particle field velocity becomes:
\({\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{x},t\right)={\tilde{\mathbf{u}}-\tilde{\mathbf{u}}}_{i}\) (14)
This definition saves the computational time by reducing the number of operations to \({2N}_{P}\). Then the relative velocity between particle \(i\) and particle \(j\) at the collision point \({\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{X}_{cp,ij},t\right)\) is used to define the Stokes number:
\(\text{St}_{\text{ij}}=\frac{m_{\text{ij}}\ {\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{X}_{cp,ij},t\right)}{6\pi\mu r_{\text{ij}}^{2}}\) (15)
In equation (15), \(m_{\text{ij}}=\left(\frac{1}{m_{i}}+\frac{1}{m_{j}}\right)^{-1}\)is the reduced mass and \(\mu\) is the carrier fluid’s viscosity. \(r_{\text{ij}}\) is the reduced radius obtained from the amplitude of curvature. Considering \(\kappa_{i}=\nabla.\mathbf{n}_{i}\mathbf{(}\mathbf{X}_{cp,ij},t\mathbf{)}\) as the local curvature of particle \(i\), in 2D the reduced radius is \(r_{\text{ij}}=\left(\kappa_{i}+\kappa_{j}\right)^{-1}\). In 3D \(\kappa_{i1}+\kappa_{i2}=\nabla.\mathbf{n}_{i}\mathbf{(}\mathbf{X}_{cp,ij},t\mathbf{)}\) and thus (\(r_{\text{ij}}=2\left(\kappa_{i1}{+\kappa}_{i2}+\kappa_{j1}+\kappa_{j2}\right)^{-1}\)) is the reduced radius.
The Stokes number helps to define the restitution coefficient. At low Stokes numbers, hardly any rebound happens. On the other hand, at very high Stokes number values the collision becomes a dry collision with wet restitution coefficient equal to one. Based on [69], the following correlation can be used to calculate the restitution coefficient \(e_{i}\) for the entire range of the Stokes numbers:
\(e_{\text{ij}}=\left\{\par \begin{matrix}1-8.65\text{St}_{\text{ij}}^{-0.75}\text{\ \ \ \ \ \ }\text{for}\text{\ \ \ \ }\text{St}_{\text{ij}}>18\par \ \ \ \ \ \ \ \ 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{for}\text{\ \ \ \ \ }\text{St}_{\text{ij}}\leq 18\end{matrix}\right.\ \) (16)
Then, the normal coefficient of restitution \(e_{N,ij}\) and tangential coefficient of restitution \(e_{T,ij}\) become:
\(e_{N,ij}=\frac{\left\langle{\tilde{\mathbf{u}}}_{\text{ij}}\mathbf{(}\mathbf{X}_{cp,ij},t\mathbf{)}.\mathbf{N}_{cp,ij}\right\rangle}{\left\|{\tilde{\mathbf{u}}}_{\text{ij}}\mathbf{(}\mathbf{X}_{cp,ij},t\mathbf{)}\right\|}e_{\text{ij}}\) (17)
\(e_{T,ij}=\sqrt{e_{\text{ij}}^{2}-e_{N,ij}^{2}}\) (18)

2-1- Normal component of the collision forces

\label{normal-component-of-the-collision-forces}
The collision force is applied by assuming a mass and spring system. The spring stiffness should stop the normal component of the particle motion in a course equal to the distance that the collision is detected and the physical particle interface. This distance is function of minimum cut off \(\varepsilon\) in Equation (5) which depends on the value of the interface thickness in Equation (1).
Unlike other collision models, we do not consider any dashpot effect here. Therefore, all of the kinetic energy of the particle in the normal direction should be transformed to the spring potential energy within the collision course \(\sigma_{\text{ij}}\):
\({\frac{1}{2}k}_{N,ij}\sigma_{\text{ij}}^{2}=\frac{1}{2}m_{\text{ij}}\left({\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{X}_{cp,ij},t\right)\right)^{2}\) (19)
Consequently, the spring stiffness in the normal direction \(k_{N,i}\ \)is
\(k_{N,ij}=m_{\text{ij}}\left(\frac{{\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{X}_{cp,ij},t\right)}{\sigma_{\text{ij}}}\right)^{2}\) (20)
and the normal collision forces become
\(\mathbf{F}_{N,i}\left(\mathbf{X}_{cp,ij},t\right)\mathbf{=-}k_{N,i}\mathbf{\delta}_{N,ij}\) (21)
We will discuss the effect of \(\sigma_{\text{ij}}\) on the collision outcomes in the section 3-2.
To include the damping effect based on the restitution coefficient it is very common to use a dashpot as a damper. However, implementation of a dashpot is engaged with the challenge of the definition of the spring stiffness and the damping coefficient [14]. A new approach is suggested here to overcome the difficulties related to implementation of a dashpot model. In this approach to consider the damping effect after the particle reaches the final collision course, the spring’s stiffness is reduced based on the normal restitution coefficient:
\(k_{N,ij,b}=e_{N,ij}^{2}k_{N,ij}\) (22)
where \(k_{N,ij,b}\) is the spring stiffness during the bounce. The collision force as an external force is added to the forces applied to the particle in each time step until there is no overlapping between the two particles. It should be mentioned that during the collision, \(k_{N,ij}\) and \(k_{N,ij,b}\) remain unchanged while \(\mathbf{\delta}_{N,ij}\) changes at each time step.

2-1- Tangential component of the collision forces

\label{tangential-component-of-the-collision-forces}
The tangential component of the contact (friction) can happen in two modes: One mode is the rolling condition when no slip condition at the contact point in tangential direction is applied by assuming a zero relative velocity and using the static friction coefficient. The other mode, sliding motion, uses the kinetic friction coefficient to find tangential forces.
(A small Literature Survey for tangential component???)
Here we follow a scheme inspired by Biegert et al. [27] which chooses the friction coefficient based on the contact mode. At first, the vector in tangential direction of the collision forces is obtained from finding the tangential component of the relative velocity:
\(\mathbf{T}_{cp,i}\left(\mathbf{X}_{cp,ij},t\right)\mathbf{=}\tilde{\mathbf{u}}\left(\mathbf{X}_{cp,ij},t\right)\mathbf{-}\left\langle{\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{X}_{cp,ij},t\right)\mathbf{.}{\hat{\mathbf{N}}}_{cp,i}\right\rangle{\hat{\mathbf{N}}}_{cp,i}\) (23)
Normalizing the tangential vector finds the unit tangential vector as
\({\hat{\mathbf{T}}}_{cp,i}\left(\mathbf{X}_{cp,ij},t\right)=\frac{\mathbf{T}_{cp,i}\left(\mathbf{X}_{cp,ij},t\right)}{\left\|\mathbf{T}_{cp,i}\left(\mathbf{X}_{cp,ij},t\right)\right\|}\) (24)
If \(\text{\ \ }\left\|{\tilde{\mathbf{u}}}_{\text{ij}}^{k+1}\left(\mathbf{X}_{cp,ij},t\right)\right\|=0\ \), \({\hat{\mathbf{T}}}_{cp,i}\left(1\right)\mathbf{=-}{\hat{\mathbf{N}}}_{i}\left(2\right)\ \)and \(\ {\hat{\mathbf{T}}}_{cp,i}\left(2\right)\mathbf{=}{\hat{\mathbf{N}}}_{i}\left(1\right)\). Next we use the static friction coefficient \(\mu_{s}\) to find friction mode as
\(\frac{\left\langle{\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{X}_{cp,ij},t\right)\mathbf{.}{\hat{\mathbf{T}}}_{cp,i}\right\rangle}{\left\langle{\tilde{\mathbf{u}}}_{\text{ij}}\left(\mathbf{X}_{cp,ij},t\right)\mathbf{.}{\hat{\mathbf{N}}}_{cp,i}\right\rangle}\left\{\par \begin{matrix}<\mu_{s}:\ \ \ \ No\ slip\par >\mu_{s}:\ \ \ \ \ Slip\end{matrix}\right.\ \) (25)
Because this model implements the normal forces form the soft collision model, the tangential part of the collision also will be a soft collision model. For the no slip mode, the friction force equals the tangential component of the hydrodynamic forces. If the slip happens
\(\mathbf{F}_{T,i}\left(\mathbf{X}_{cp,ij},t\right)\mathbf{=-}\mu_{k}\ \mathbf{F}_{N,i}\left(\mathbf{X}_{cp,ij},t\right)\) (26)
where \(\mu_{k}\) is the kinetic friction coefficient.
Having the normal and tangential components of the collision forces, the moment induced by the collision is obtained from
\(\mathbf{T}_{c,i}\mathbf{=}\left(\mathbf{X}_{cg,i}-\mathbf{X}_{cp,ij}\right)\mathbf{\times}\left(\mathbf{F}_{T,i}\mathbf{+}\mathbf{F}_{T,i}\right)\) (27)
where \(\mathbf{X}_{cg,i}\) is the coordinate of the particle’s center of gravity.

2-3- Fluid-Structure Interaction

\label{fluid-structure-interaction}
To allow small amount of solid-solid overlapping of the interfaces, we use a diffuse interface method to apply the fluid-structure interaction. As described in [70-72], we see the Smoothed Profile Method (SPM) as a proper choice for the simulation of flows laden with particles because of higher stability and efficiency. The reader is referred to [13,
After SPM finds the forces and moments from the flow, the update of velocity and angular velocity will be based on
\(\frac{d\mathbf{V}_{i}}{\text{dt}}=\left(\mathbf{F}_{h_{i}}+\mathbf{F}_{\text{ext}_{i}}+\mathbf{F}_{N,i}+\mathbf{F}_{T,i}\right)/M_{i}\) (27)
\(\frac{\mathbf{\text{dΩ}}_{i}}{\text{dt}}=\frac{\left(\mathbf{T}_{h_{i}}+\mathbf{T}_{\text{ext}_{i}}+\mathbf{T}_{c,i}\right)}{I_{i}}\)
Here, \(M_{i}\ \)and \(I_{i}\) are mass and moment of inertia. Subscript \(h_{i}\) represents the force and moment from the fluid, and the subscript \(\text{ext}_{i}\) is for external force. Similarly, particle position and orientation are updated by integration of velocity and angular velocity: