Modelling the normal bouncing dynamics of spheres in a viscous fluid

. Bouncing motions of spheres in a viscous fluid are numerically investigated by an immersed boundary method to resolve the fluid flow around solids which is combined to a discrete element method for the particles motion and contact resolution. Two well-known configurations of bouncing are considered: the normal bouncing of a sphere on a wall in a viscous fluid and a normal particle-particle bouncing in a fluid. Previous experiments have shown the effective restitution coefficient to be a function of a single parameter, namely the Stokes number which compares the inertia of the solid particle with the fluid viscous dissipation. The present simulations show a good agreement with experimental observations for the whole range of investigated parameters. However, a new definition of the coefficient of restitution presented here shows a dependence on the Stokes number as in previous works but, in addition, on the fluid to particle density ratio. It allows to identify the viscous, inertial and dry regimes as found in experiments of immersed granular avalanches of Courrech du Pont et al. Phys. Rev. Lett. 90 , 044301 (2003), e.g. in a multi-particle configuration.


Introduction
Experiments on the bouncing motion of particles in a fluid have been carried out during the past decade.[1][2][3] and [4] studied the normal rebound in a fluid of a particle on a wall and of two particles respectively.In both configurations, the effective restitution coefficient ߝ = −ܸ ோ /ܸ ் , defined as the ratio between the relative velocity of rebound V R and the relative terminal velocity before contact V T , is lower than the restitution coefficient in the dry case, noted ߝ ௫ .The presence of the fluid squeezed between the objects during the bouncing motion tends to dissipate the initial kinetic energy of the particle system.The ratio ߝ/ߝ ௫ is shown to be a function of a single parameter, namely the Stokes number St which compares the inertia embedded in the bouncing motion, which is composed of the particle inertia and added mass of the fluid (see [5]), with the fluid viscous dissipation and reads ‫ݐܵ‬ = ൫ߩ + 0.5ߩ൯ܸ ் ‫)ߤ9(/ܦ‬ where ߩ and ߩ are the particle and fluid density, respectively, ‫ܦ‬ is the sphere diameter and ߤ is the dynamic viscosity of the fluid.Experiments of a dense avalanche in a fluid by [6] revealed that beyond the Stokes number, the particle to fluid density is also a pertinent dimensionless number to distinguish the observed granular regimes, namely viscous, inertial and dry.Interestingly, one may note that the dependence on the density ratio has not been explicitly reported in the bouncing configuration, for which the reported coefficient of restitution ߝ was shown to vary only with ‫.ݐܵ‬ In the present work, we simulate the bouncing motion of one and two spheres using a coupled fluidparticle method (see [7]).An immersed boundary method (IBM) resolves the fluid flow in the presence of moving solid grains while the discrete element method (DEM) models the contact and lubrication interactions between particles.Simulations allow to investigate fluid to particle density ratio in the range ࣩ (10) to ࣩ(10 ଷ ) difficult to explore in experiments due to the nature of the fluid and particulate elements.This article is structured as follows.First, the numerical technique is described.Then, simulations of normal bouncing of a sphere on a wall and of two spheres in a fluid are presented.We finally discuss the definition of the effective restitution coefficient of particles bouncing in a fluid.

Numerical methods
In this section, the coupled IBM/DEM is briefly presented.The reader may refer to [7] for more details and validation about the numerical approach.

Fluid flow calculation
The IBM (see [8]) solves the following modified Navier-Stokes equations (1)-(2) for a Newtonian fluid in the entire computational domain, e.g. in the fluid domain as well as inside the particles, using a Cartesian system of coordinates, and ܲ are the local fluid velocity and pressure, is the gravity vector and = ߙ( − ‫ݐ∆/)‬ is a body-force source term formulated such as it accounts for the presence of the solid particle in the fluid.and ߙ denote the local velocity of the solid object and the solid volume fraction, respectively, and ‫ݐ∆‬ is the fluid time step.Note that in the fluid domain, = and the regular Navier-Stokes equations are solved.The hydrodynamic force applied by the fluid upon the particle ‫‬ of volume ߴ is calculated as

Particle motion calculation
In the case of normal bouncing, the particles motion is calculated through the resolution of the Newton's equations (4) for the particle linear momentum which reads where ‫ݑ‬ and ݉ are the translational velocity and particle mass of the particle ‫,‬ respectively.The total force applied on the solid particle can be divided in four contributions, namely the particle weight; the hydrodynamic force , defined in (3); the normal contact force and a lubrication force which are defined in the following.Solid-solid interactions are modelled via a DEM method, namely the soft-sphere method [9].Let us consider the normal unit vector n which links the centre of mass of the two solids embedded in a contact interaction.An overlap ߜ between the solid objects during the contact is allowed and permits to compute the normal contact force = ‫ܨ‬ with a damped massspring model as follows For a contact interaction between two objects ݅ and ݆ of mass ݉ and ݉ , ݇ and ߛ are the normal stiffness and the damping coefficient, respectively, and are defined as a function of the dry restitution coefficient ߝ ୫ୟ୶ and the contact time ‫ݐ‬ by ߛ = −2݉ * ݈݊(ߝ ௫ )/ ‫ݐ‬ and ݇ = ݉ * ߨ ଶ ‫ݐ/‬ ଶ + ߛ ଶ /4݉ * where ݉ * = ݉ ݉ /(݉ + ݉ ) is the effective mass.If one considers a contact interaction of the particle ݅ with a wall, ݉ * is reduced to ݉ .
The flow structure may not be entirely captured by the IBM when two solid objects approach or separate from each other.As shown by several authors (see [10][11] for instance), this issue is overcome using a lubrication force = ‫ܨ‬ which is defined as where ܴ * = ܴ ܴ /(ܴ + ܴ ) is the effective radius of the particles ݅ and ݆ and ߟ is a so-called effective roughness length which accounts for particle surface asperities and can range in [10 ି ܴ * ; 10 ିଷ ܴ * ] [11].Note that ܴ * is reduced to ܴ , if one considers the lubrication interaction of the particle ݅ with a wall.A time step ‫ݐ∆‬ = ‫ݐ‬ /100 is set to resolve the time integration of (4).Images are inspired from [11] and have been mirrored.
In the following, t denotes the dimensionless time for which the characteristic time ‫)݃/ܦ(‬ ଵ/ଶ has been used.A typical time evolution of the flow structure near bouncing is depicted in Fig. 1 while the particle velocity versus time is plotted in Fig. 2. At ‫,91~ݐ‬ the particle reaches the near wall area at a constant velocity ܸ ் .Then at ‫ݐ‬ = 19.89, the particle starts to decelerate due to the fluid film drainage characterized by some vorticity appearing near the wall surface.The particle inertia is large enough for the particle-wall rebound to occur with the considered particle roughness.The bouncing occurs at ‫ݐ‬ = 20.02.Afterwards, the particle takes off from the wall but is decelerated by (i) its own weight, (ii) the drag induced by the fluid entrained in the particle wake and (iii) the lubrication force stemming from the fluid film recovering from the particle and the wall.On a short time scale just after the bouncing the deceleration is mostly controlled by (ii) and (iii).
Fig. 2. Time evolution of the sphere velocity (same case as in Fig. 1).The terminal velocity ܸ ் and the rebound velocities ܸ ோ and ܸ ௌ are represented.The case of a sphere labelled 1 approaching another one labelled 2 initially at rest in a fluid is now investigated with the same numerical tool.The grid resolution is ‫ݔ∆‬ = ‫.02/ܦ‬Here, gravity is set to zero.An example with ߩ /ߩ = 8, ‫ݐܵ‬ = 140 and ܴ݁ ≅ 175 is shown in Fig. 3 and Fig. 4 where the structure of the flow and the time evolution of the velocity of the two particles are presented near the bouncing event, respectively.At ‫ݐ‬ = 0.24, particle 2 is not influenced by particle 1. Vorticity is generated on particle 2's surface at ‫ݐ‬ = 0.71 as particle 1 is approaching at a distance which is less than a radius.The bouncing starts at ‫ݐ‬ = 1.18 after which part of the kinetic energy of particle 1 is transfered to particle 2. Afterwards, the separation of particle 2 from particle 1 entrains fluid inbetween the two particles.At later times, particles move with slowly decreasing velocities during the time of the simulation, due to viscous dissipation at the particle surface.

Discussion about the definition of the effective coefficient of restitution
During particle rebound, several characteristic velocities may be defined.First, as for the bouncing of one particle on a wall, we classically define the terminal velocity ܸ ் as the velocity when the particle is not influenced by the wall presence and the rebound velocity ܸ ோ as the particle velocity just after contact (see Fig. 2).A coefficient of normal restitution of bouncing can be defined as ߝ = −ܸ ோ /ܸ ் and only depends on the Stokes number [1,2].Additionally, one may remark a drastic change in the particle acceleration after contact which arises from the various forces acting on the particle during taking off.A velocity ܸ ௌ may then be defined as a second characteristic velocity of rebound (see Figs. 2).The corresponding effective restitution coefficient can then be defined as ݁ = −ܸ ௌ /ܸ ் .In the case of the bouncing motion of two spheres, [4] show that the normalised coefficient of restitution also depends on the Stokes number defined as ‫ݐܵ‬ = ߩ ‫ܸ(ܦ‬ ்ଵ − ܸ ்ଶ )/(9ߤ) with ܸ ் the velocity of the ݅-th particle before bouncing.Again, no influence of the density ratio was reported in this study.As in the case of the normal bouncing of a single particle on a wall, one can define in this case a restitution coefficient of bouncing of two particles as ߝ = −(ܸ ோଵ − ܸ ோଶ )/(ܸ ்ଵ − ܸ ்ଶ ) where ܸ ோ is the rebound velocity after contact (see Fig. 4).Again, we remark a change in the particle acceleration when particle 2 takes off after contact.At that time, we define the velocities ܸ ௌଵ and ܸ ௌଶ which stand for a second characteristic velocity of rebound of particle 1 and 2 (see Fig. 4).We here propose an alternative definition of the effective restitution coefficient of bouncing for two particles, namely ݁ = −(ܸ ௌଵ − ܸ ௌଶ )/(ܸ ்ଵ − ܸ ்ଶ ).Note that same notations for ߳ and ݁ are used between the different bouncing configurations for clarity.
As shown in Fig. 5a, the present IBM/DEM simulations are in good agreement with available experimental data in both configurations confirming the exclusive dependence of ߝ/ߝ ௫ with the Stokes number.Nevertheless, this exclusive Stokes dependency is not observed with the new effective coefficient of restitution ݁, proposed here, scaled by ߝ ௫ (see Fig. 5b).One may observe some plateaus for Stokes numbers higher than approximately 300 which depend on the particle to fluid density ratio ‫ݎ‬ ௗ = ߩ /ߩ.A new scaling of ݁ proposed here is ߝ = ߝ ௫ /(1 + ‫ݎ‬ ௗ / ‫ݎ‬ ௗ ) with ‫ݎ‬ ௗ = 13.This scaling makes the data collapse on a master curve for ݁/ ߝ as a function of the Stokes number (see Fig. 5c).With this new definition of ݁, the energy dissipation during bouncing is shown to be dependent not only on the dry restitution coefficient ߝ ௫ and the Stokes number ‫,ݐܵ‬ as is the case for ߝ, but also on the particle to fluid density ratio ‫ݎ‬ ௗ .Note that the quite simple influence of ߝ ௫ reported here is simply generalized from the results obtained from ߝ similar to the ones obtained in the experiments.However, this dependency would probably deserve a deeper investigation (both experimental and numerical) regarding the new definition of the effective restitution.].Experiments of [1] and [4] are represented by • and •, respectively.The solid line is the analytical model proposed in [11] for / which was straightforwardly used for / .This new definition of the coefficient of restitution and in particular its new scaling allow to make an analogy with the regimes observed in immersed granular avalanches [6].Firstly, when ‫ݐܵ‬ ≪ ‫ݐܵ‬ where ‫ݐܵ‬ is some critical Stokes number of the order of 10, there is no rebound and ݁ = 0 regardless of ߩ /ߩ.This may be interpreted as the viscous regime.Alternatively, when ߩ /ߩ ≫ ‫ݎ‬ ௗ , there is no influence of the density ratio on ݁ which depends only on the Stokes number.In this case, when ‫ݐܵ‬ ≫ ‫ݐܵ‬ we have ݁ = ߝ ௫ and the regime is the so-called dry or free-fall regime.Finally, when ߩ /ߩ~ ‫ݎ‬ ௗ and ‫ݐܵ‬ ≫ ‫ݐܵ‬ , the regime would presumably be inertial as in the immersed granular avalanches of Courrech du Pont et al. [6].The present estimation of the density ratio ‫ݎ‬ ௗ = 13 corresponding to the transition between the inertial and dry regimes is close to the experimental value of 16 predicted in [6].

Conclusion
The bouncing of a sphere on a wall or another sphere in a viscous fluid has been simulated with an IBM/DEM method.The method is shown to adequately reproduce experimental results for the whole range of investigated parameters.A novel effective restitution coefficient is proposed which accounts for both the solid contact and the whole hydrodynamics effects just before and after the bouncing.With this coefficient both ߩ / ߩ and the Stokes numbers are observed to play a role in the dynamics of bouncing in line with the viscous, inertial and dry regimes found in avalanche experiments of Courrech du Pont et al. [6].The analytical model proposed in [11] for ߝ/ߝ ௫ has been straightforwardly extended to account for this influence of the density ratio by defining a new effective scaling factor such as: ݁/ߝ with ߝ = ߝ ௫ /(1 + ‫ݎ‬ ௗ ‫ݎ/‬ ௗ ), ‫ݎ‬ ௗ = ߩ /ߩ and ‫ݎ‬ ௗ = 13.

Fig. 4 .
Fig. 4. Time evolution of the particle velocities (same case as in Fig. 3) during the bouncing of two spheres.The terminal velocity ܸ ் and rebound velocities ܸ ோ and ܸ ௌ of particle ݅ are represented.