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: