3- Results and Discussion

\label{results-and-discussion}
In this section at first the collision model is validated by comparison with the experimental results for a particle wall collision. Then, the robustness of the method is studied for various cases of non-spherical particles including in collisions with friction and without friction. Finally, the model is tested for the simulation of particle-particle cases.

3.1. Validation: Particle-Wall Collision

\label{validation-particle-wall-collision}
In order to validate the method for the collision with a stationary wall, the experimental results of Gondret et al. [76] which were also implemented later in [27, of steel with \(D_{p}=3\ mm\), \(\rho_{p}=7800\ kg/m^{3}\), \(\nu=0.3\) and \(E=2.4\times 10^{11}\) is released for a free fall in silicon oil RV10 with \(\mu=0.01\ Pa.s\) and \(\rho_{f}=935\ kg/m^{3}\). Here, \(d_{p}\), \(\rho_{p}\), \(\nu\), \(E\) are respectively solid particle’s diameter, density, Poisson’s ratio, and Young modulus. \(\rho_{f}\) and \(\mu\) are the fluid’s density and viscosity. The simulation is performed in non-dimensional form by normalizing the distance with the particle diameter and the non-dimensional parameters \(Re=164\), \(Fr=11.8\) and \(\frac{{\gamma=\rho}_{p}}{\rho_{f}=8.342}\). The simulation is made in a \(10\times 10\times 10\ \)domain with different minimum grid sizes to study the effect of grid refinement. The cut off value for the overlapping function is set to 10% of the particle diameter. Consequently, when a particle approaches the wall, the course of deceleration is \(\sigma=0.1\).
The particle starts to fall and reaches the final steady velocity before hitting the wall. It starts bouncing and losing the kinetic energy and finally it stays motionless on the wall. The very low or zero spring stiffness values related to the final stage of the collision means that the particle may not have enough collision force to hold it over the wall after ultimate reduction of spring coefficient. This happens when the bouncing height is so low that the particle cannot escape the overlapping area. Therefore, the spring coefficient is not updated to higher \(k_{N,ij}\) values related to the next bounce as it still has the low bounce spring coefficient value \(k_{N,ij,b}\). To avoid penetration of the particle to the wall in these cases, an extra upward force equal to negative of weight is applied.
The results are presented in Figure (1) and they show the overall agreement of the method with the experiment. Moreover, the independency of the results from the grid also is approved. There is a small discrepancy among the current results and the experiment which is more noticeable in the velocity figure after the second collision. The maximum velocity just after the collision is smaller in the simulation than the experiment. This is because in the current soft collision model, the effect of viscous dissipation is applied through the restitution coefficient. In reality there is also viscous dissipation through the computations which becomes significant at low Stokes number values. For this reason, Kempe and Fröhlich [14] ignored the hydrodynamic drag on the particle during collision course. The discrepancy it is not apparent in the first collision because of the short duration of collision at high impact velocity which results in a harder collision. However, it grows in the subsequent collisions with diminution of the impact velocity and lower Stokes numbers. Moreover, the accumulation of errors by time also is the other reason for increase in the discrepancy.

3-2. Effect of Collision Course

\label{effect-of-collision-course}
In the presented collision method, the collision course is the only parameter that is not defined. We defined the collision course, \(\sigma\), as the distance between the point with \(\left.\ \text{dψ}_{i}\left(\mathbf{x},t\right)\right|_{\text{cut\ off}}=\varepsilon\) and the particle interface. This section studies the dependency of the results to \(\sigma\). A problem like the one in the section 3.1 is considered but with different collision courses (Figure 2). The results are related to minimum grid size of \(\Delta x=0.075\) and \(CFL=0.25\).
For \(\sigma=0.3\), the error related to the simulation is not significant as it is mostly limited to the instants of impact. After the collision .i.e. when particle leaves the overlapping region, it reaches the actual path and the error is partially compensated. This occurs in spite of the fact the collision now has become too soft and the sudden changes in velocity are smoothed. Moreover, after the second collision, the particle does not escape the very wide overlapping region and the error increases.
Variation of height with time for the first collision of the particle is represented in Figure (2c) for the particle using different collision courses. When \(\sigma=0.3\) the particle does not reach the wall completely leading to emerge of the small error. For \(\sigma=0.05\), the difference between the simulation and experiment becomes noteworthy by showing a higher restitution coefficient for the collision. In this case, the collision duration is so short that the simulation becomes sensitive to the small errors. As it could be observed in Figure (2c), even though the particle path is very close to the experiment, after the collision it bounces with higher velocity.

3-3- Collision of an Elliptical Particle with a Frictionless Wall

\label{collision-of-an-elliptical-particle-with-a-frictionless-wall}
When a particle is not spherical, the rotation becomes important even if the walls are frictionless. This is because in general the collision normal vectors do not point to the center of gravity. Here, we consider the collision of a free-falling ellipse on a frictionless wall. The ellipse with \(D_{p}=(2.0,0.5)\) and \(\frac{\rho_{p}}{\rho_{f}=10}\) is located at \((5.0,7.0)\) in a \(10\times 10\) channel with four no-slip walls. The particle is released and settled under the gravity. At Froude number \(Fr=0.51\), we study two Reynolds numbers of \(Re=20\) and \(Re=2000\). The Reynolds and Froude numbers are measured with respect to the major diameter of the particle just before the collision happens. The initial angle of the particle is \(\theta=\frac{\pi}{4}\) with respect to the horizontal wall. Due to high asymmetry of the problem, the particle should start to rotate after the start of motion. However, to avoid more complexities, the rotational motion of the particle is restricted until the collision occurs. During the falling time, the particle is forced to fall with a unit downward vertical velocity. The minimum grid size is \(x=0.0375\) and \(CFL=0.2\).
Figure (3) shows the sequence of the particle position in the channel. Due to the tilted placement of the particle at the time of impact, the particle reaches the wall on its sharp edge. Then it bounces slightly and rotates towards the final position of the particle with minimum gravitational potential on the flat side of the particle.
Figure (4) represents the time history of the particle kinematical variables. It is clear from Figures (3) and (4) that at \(Re=20\) after the impact of the particle with wall and a small vertical bounce, the particle settles down to the final vertical position very fast. This is because of the high dissipation at low Re. The other hand, the particle continues its motion in the horizontal direction by sliding over the frictionless wall after gaining a horizontal momentum. This momentum damps due to dissipation in the fluid. At \(Re=2000\), the lower dissipation of energy leads a number of subsequent collisions after the first collision as well as a horizontal motion in the other direction.

3-4- Tangential Collision Model

\label{tangential-collision-model}
To validate the tangential collision model, we study a case similar to Biegert et al. [27]. A sphere with \(D_{p}=0.125\) and \(\frac{\rho_{p}}{\rho_{f}=2.5}\) is exposed to a linear shear flow in a \(1\times 1\times 1\) domain. Two cases in different fluids with \(\nu=0.1\) and \(\nu=0.02\) are considered. The top wall is driven with \(u_{w}=1\) leading to \(\text{Re}_{w}=10\) and \(\text{Re}_{w}=50\). The static and dynamic friction coefficient as \(\mu_{s}=0.8\) and \(\mu_{k}=0.15\) related to silicate materials are considered. The linear velocity profile in the domain varies with height and is applied to all boundaries except the outlet which satisfies the condition of zero gradient for the variables. The particle is placed on the bottom wall and is subject to gravity of \(g=9.81\). A minimum grid size of \(0.05\) is considered with \(CFL=0.5\).
The particle is initially at stationary and with the start of the simulation the collision is detected. It is held for \(T=\frac{tu_{w}}{H}=0.01\ \) for the flow to develop slightly and avoid the spikes in the forces at the start of the simulation. Then the particle is released leading to acceleration and motion to reach the steady condition. We compared our results with Biegert et al. [27] in Figure (5). While the results are very similar, there is a small discrepancy which is probably because of using a different diffuse interface solver. Figure (5) also shows that at similar to Biegert et al. [27], at \(\text{Re}_{w}=10\) the translational velocity (\(U_{p}\)) is more than the rolling velocity (\(R_{p}\omega_{p}\)) which means that the particle is sliding on the wall. However, at \(\text{Re}_{w}=50\), the particle does not slide and a rolling condition happens because \(U_{p}=R_{p}\omega_{p}\). Here, we did not distinct between \(R_{p}\omega_{p}\) and \(U_{p}\) as they are very similar.
Figure (6) represents the streamlines and the vortex contours on the symmetric vertical plane parallel to the flow. It shows that the contours are very similar for the two Reynolds number values but at \(\text{Re}_{w}=50\) there is a small asymmetry due to the wake effect. The streamlines also show that the particle is in rolling condition at \(\text{Re}_{w}=50\) while it slips significantly at \(\text{Re}_{w}=10\).

3-5- Collision of Rectangular Particle with Wall

\label{collision-of-rectangular-particle-with-wall}
To study the effect of tangential part of the collision in this section we consider a rectangular block of copper with \(2\times 0.5\) dimensions standing vertically on a wall. We consider the bottom wall with zero velocity and in two situations: one case without friction and the other one with static and dynamic friction coefficients similar to copper on steel friction i.e. \(\mu_{s}=0.53\) and \(\mu_{k}=0.36\). A uniform unit velocity is applied for the input and the top wall in the \(x\) direction. The particle is located vertically at \((3,1)\) in a \(15\times 6\) channel. With \(\frac{\rho_{p}}{\rho_{f}=8.96}\), \(Re=\frac{U_{\infty}L}{\nu}=100\) and \(Fr=\frac{U_{\infty}^{2}}{\text{gL}}=0.102\) related to \(\nu=0.01\), \(g=9.81\) and \(L\) as the average dimension of the rectangle. The rectangle is fixed at first until \(T=0.1\) to have the flow developed over the particle and avoid a sudden spike in force at the start of the simulation.
The sequence of particle motion is shown in Figure (7) for both friction and frictionless collisions. Figure (8-9) display the quantitative variation in particle’s motion for the two cases. With friction, the particle shows a limited motion which leads to a final stop at about \(T\sim 2.2\). Without friction, it continues to move to reach the end of the channel with a steady velocity of \(u_{\text{final}}\sim 0.56\). Even though the contact is frictionless, the particle does not reach the freestream velocity because of the boundary layer on the wall. Moreover, the particle is losing momentum by time as it has less flow blockage in comparison with the start of the motion.
The differences between the two cases is not limited to the translational motions as the two particles have different rotational motions. With friction, the particle pivots around the front collision point and rotates until it contacts the wall on the longer side and becomes stable. On the other hand, without friction, the particle is free to rotate around any of the bottom edges. Therefore, when the applied moment on the particle is positive due to more blockage of the flow on the bottom, the particle rotates on the other direction leading a positive angular velocity. However, for both cases the rotation is restricted because after the particle reaches the long edge, the stabilizing effect of the weight does not allow to particle to rotate more.
The effect of grid refinement on the variation of angle is studied in Figure (10) and it shows that with refinement of the grid, the results follow the same overall trend. However, there is a small discrepancy for the small fluctuations on the angle which belongs to the time that the particle hits the wall. This time is the time that the collision forces become dominant and the results are more sensitive to the initial conditions and also the coarse grid.
Figures (11-12) display the contour plots of vorticity for the frictionless and frictional collision at different time instants. The evolution of vorticity from the tip of the rectangle is noticeable which in the case of the frictionless collision is attached to the tip and moves with the flow and the particle. For the frictional collision, it detaches the tip due to fast downward motion while another vorticity is generated on the other side of the particle.

3-6- Collision of Two Particles

\label{collision-of-two-particles}
Now, we study a drafting, kissing, and tumbling case similar to Glowinski et al. [77] who used a repulsion force to apply the normal collision for the disk shape particles. Two similar particles start to fall under the gravity of \(g=981\) in a vertical channel with size \(2\times 6\) containing a fluid with \(\nu=0.01\). The particles with \(d_{p}=0.25\) and density ratio of \(\frac{\rho_{p}}{\rho_{f}=1.5}\ \)are initially located at \((1,4.5)\) and \((1,5)\). Having the same horizontal position, motion of the leading particle generates a wake which results in a drag reduction on the trailing particle. Therefore, the trailing particle will have greater velocity and reach the leading particle. After the collision, the instability in the system leads to deviation of the particles from the centerline of the channel.
Here, a minimum grid size of \(x=0.0125\) with \(CFL=0.5\) is considered. Figures (13-14) shows the results for a collision course equal to 10% of the particle diameter. It is clear from the results that having a diffuse interface and a soft sphere collision model with a thick collision course results in an error in velocity during the collision period. This error will be compensated for the vertical velocity after the collision period. For the horizontal velocity, the small errors in collision accumulates and leads to greater deviation because of high sensitivity to the initial conditions. In overall, the error in the position of the particles is not significant because the change in the direction of velocity components results in compensation of error in the position as the integral of the velocity.