International Journal of Multiphase Flow 46 (2012) 38–53 Contents lists available at SciVerse ScienceDirect International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w Modeling mechanical contact and lubrication in Direct Numerical Simulations of colliding particles Julian A. Simeonov ⇑, Joseph Calantoni Marine Geosciences Division, Naval Research Laboratory, Code 7434, Bldg. 1005, Stennis Space Center, MS 39529, United States a r t i c l e i n f o Article history: Received 2 February 2012 Received in revised form 4 May 2012 Accepted 9 May 2012 Available online 21 June 2012 Keywords: Particle-laden ﬂow Direct Numerical Simulations Lubrication hydrodynamics Contact mechanics a b s t r a c t We developed a model for inexpensive Direct Numerical Simulations of particle-laden ﬂow by fully resolving the hydrodynamics at all times except when the gap between colliding particles becomes comparable to the grid step. The resolved hydrodynamics were obtained with a previously developed pressure boundary integral method for direct ﬂuid–particle simulations on Cartesian grids. The unresolved part of the lubrication pressure/shear force in the subgrid gap was predicted using theoretical Stokes ﬂow models. Singular lubrication forces were avoided by postulating that contact begins when the gap distance between the particles becomes equal to the size of micro-asperities. The mechanical contact between particles is dynamically resolved and particle interactions are both inelastic and frictional. The proposed numerical model was validated through resolution tests and comparison with experimental data for immersed binary collisions. Published by Elsevier Ltd. 1. Introduction The properties of dense particle suspensions can be strongly affected by inelastic and frictional losses during particle collisions. Our ﬁrst objective is to include a soft-sphere model (Walton and Braun, 1986) for mechanical contact in Direct Numerical Simulations (DNS) of ﬁnite size colliding particles. Frictional contact effects are essential for adequate modeling of oblique collisions (Foerster et al., 1994) but do not appear to have been included in previous DNS models. In earlier DNS models (Glowinski et al., 1999; Hu et al., 2001) the mechanical interaction was completely prevented by adding an ad-hoc action-at-a-distance repulsive force that keeps the minimum particle distance greater than a grid step. The reason for this restriction was to avoid the problem of choosing which of the two boundary conditions (of two colliding particles) to enforce at overlapping grid points. More recently, a modiﬁed collision method was proposed that allows particle overlap and replaces the action-at-a-distance force with an elastic force that activates only when the overlap is non-zero (Singh et al., 2003). The latter approach is similar to the soft-sphere model for normal contact but the values of the stiffness parameter were motivated by the numerical objective of keeping the overlap smaller than a grid step. A frictionless hard-sphere DNS model was introduced where upon contact the normal velocity component changes instantaneously to its rebound value based on the dry coefﬁcient of ⇑ Corresponding author. Tel.: +1 228 688 4206; fax: +1 228 688 4853. E-mail address: [email protected] (J.A. Simeonov). 0301-9322/$ - see front matter Published by Elsevier Ltd. http://dx.doi.org/10.1016/j.ijmultiphaseﬂow.2012.05.008 restitution (Hu et al., 2001; Ardekani and Rangel, 2008); however, a rotational coefﬁcient of restitution would be needed to realistically model oblique collisions. While the hard-sphere approach is very attractive because of the low computational cost associated with not having to resolve the extremely short duration of real collisions, it presents numerical issues related to singular hydrodynamic effects. First, after collisions the hard-sphere approach gives rise to an (initially irrotational) impulsively started ﬂow which has a tangential velocity discontinuity at the boundary; the initially inﬁnite boundary layer gradients correspond to a singular viscous force which cannot be properly resolved with a ﬁnite numerical discretization. Second, the inﬁnite particle acceleration during such instantaneous collisions corresponds to an inﬁnite added mass drag (of vanishing duration) and the corresponding ﬂuid impulse is indeterminable. Collisional added mass effects apparently are not accounted for in hard-sphere models but may be important. The signiﬁcance of various hydrodynamic effects (e.g. added mass and wake interaction) for dissipating the particle momentum during collisions will be examined in detail with the present model. Our second objective is to develop an inexpensive Cartesian grid model for direct ﬂuid–particle simulations where an ‘‘outer’’ solution is computed by fully resolving the hydrodynamics on the particle scale, except in the gap between colliding particles when the latter becomes smaller than a few grid steps. The numerically computed outer solution is determined with our pressure boundary integral method (Simeonov and Calantoni, 2011), which exploits the computational efﬁciency of fast Poisson solvers on uniform grids. In a different approach based on non-uniform grids, the J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 lubrication ﬂow in the gap could be resolved explicitly by using very ﬁne meshes in the gap, but except for Johnson and Tezduyar (1996) such simulations have generally been limited to 2D (e.g. Hu, 1996; Ardekani and Rangel, 2008) due to high computational costs. The body-ﬁtted mesh approach will also be computationally prohibitive for 3D simulations with a large number of particles. Proper account for the micro-scale lubrication effect is essential to correctly predict the particle viscous momentum losses during collisions characterized with a low collisional Stokes number, St ¼ m ðU 2 U 1 Þ ; 6pla2 ð1Þ where m⁄ and a⁄ are the reduced particle mass and radius, respectively, (U2 U1) is the relative particle velocity, and l is the dynamic viscosity. It is well known (e.g., Davis et al., 1986; Barnocky and Davis, 1988; Yang and Hunt, 2006) that a critical Stokes number exists (St 10) below which there is no rebound after collision in viscous liquids. In the present model, the lubrication forces/torques due to the unresolved ‘‘inner’’ solution will be approximated with analytical models (Cooley and O’Neill, 1969; O’Neill and Majumdar, 1970). We have not attempted to include elastohydrodynamic effects (Davis et al., 1986) but our phenomenological lubrication model includes the effect of cavitation during rebound. We should also mention previous work (Maury, 1997; Dance and Maxey, 2003) where lubrication forces were added as a correction to numerical models of particle suspensions; however, these models did not allow mechanical contact and are limited to low Reynolds number applications. The next section begins with a short introduction of our numerical method for modeling ﬂuid–particle hydrodynamics on Cartesian grids. The force models for mechanical contact and lubrication are also presented in Section 2. In Section 3, we evaluate the robustness of our collision model with resolution tests and comparisons with the normal collision experiments of Yang and Hunt (2006), hereafter YH06. DNS of oblique collisions also were tested against the corresponding experiments of Foerster et al. (1994), hereafter FO94, and YH06. The results are discussed in Section 4 and our ﬁndings are summarized in the Conclusions section. 2. Methods 2.1. Hydrodynamics solver In this section we provide a short account of the pressure boundary integral method (Simeonov and Calantoni, 2011) for computing the particle scale ﬂuid ﬂow and pressure generated by the motion of spherical particles of radius, a, over a regular Cartesian grid. The ﬂuid velocity, u, and pressure, p, are determined by solving the incompressible Navier-Stokes equation in the pressure Poisson equation formulation, @u rp þ u ru ¼ þ mDu; @t q ð2aÞ Dp ¼ qr ðu ruÞ; ð2bÞ 39 where Uj and Xj are the particle translational and rotational velocities. In (3a) ir (sinhcos/, sinhsin/, cosh) is the local unit radial vector and (h, /) are the local spherical coordinates on the particle surface. The boundary condition for (2b) is a Neumann boundary condition given by the normal component of the momentum equation evaluated on the boundary. On the particle surface exterior, the Neumann boundary condition is 1 @p @ 2 ur dUj ¼ m 2 ir þ aðX2j ðXj ir Þ2 Þ; q @r @r dt ð3bÞ where ur u ir is the radial velocity component. Using fast Poisson solvers for the pressure implies solving over the whole Cartesian grid and therefore requires extending the pressure ﬁeld across the particle boundaries. Noting that there is complete freedom in the manner of specifying the ﬁctitious interior pressure, we make the simplest choice where the interior pressure ﬁeld is harmonic (i.e. Dp = 0). With this choice, the Laplacian of pressure is discontinuous across the boundaries of particles, except when the particle rotational velocity vanishes. Furthermore, either the pressure or its gradient will be discontinuous because only the function or its derivative can be speciﬁed as a boundary condition for the Poisson equation. In our implementation we assume that the pressure is continuous across particle boundaries and therefore the pressure gradient is discontinuous. The discontinuity may be represented by a single layer potential, where the unknown single layer potential density (in this case the pressure gradient jump) is related to the Neumann boundary condition (3b) via boundary integral equations (Simeonov and Calantoni, 2011). We solve the system of boundary integral equations by expanding the single layer density distribution on the particle surface into surface spherical harmonics. The expansion reduces the system of boundary integral equations to a diagonally dominant linear algebra problem, which is solved with Gauss–Seidel iterations. Having determined the pressure gradient jump and knowing the jump in Dp by construction, we can evaluate the discontinuities in the pressure and its Cartesian derivatives up to second order derivatives. Following Mayo (1984) the computed discontinuities are then used to correct the second-order ﬁnite difference Laplacian (so that second order accuracy is preserved across discontinuities) for stencils intersecting particle boundaries and the resulting ﬁnite-difference Poisson equation is solved with fast algorithms on the Cartesian grid. This completes the solution for the pressure ﬁeld on the Cartesian grid and allows computing the pressure forcing terms in the ﬂuid (2a) and the particle momentum equations (below). The ﬂuid momentum Eq. (2a) is reduced to an ODE system using second order ﬁnite difference approximations of the spatial derivatives; central differences are used for Cartesian nodes far from the particle surface and one-sided differences are used for nodes next to the boundaries. The ODE’s are then integrated in time using a 4th order Runge-Kutta scheme. The second-order spatial accuracy of this method was established previously (Simeonov and Calantoni, 2011). 2.2. Particle motion where m is the viscosity and q is the density. The above formulation was motivated by the desire to exploit the numerical efﬁciency of fast Poisson solvers on uniform Cartesian grids with a grid step h. The boundary condition for the momentum Eq. (2a) is the no-slip condition on solid boundaries or periodic on periodic boundaries (if any). On the jth particle surface the no-slip boundary conditions implies, u ¼ Uj þ aXj ir ; ð3aÞ The translational and rotational equations of motion for a spherical particle with uniform density, qp, mass, m, and moment of inertia, I = 2ma2/5, subject to hydrodynamic forces, a body force, g, inter-particle mechanical contact, FM, and lubrication forces, FL, are given by m dUj q ¼ mg 1 qp dt þ FM þ FL ; ! þ I @uu @uh pir l ih l iu dS @r @r ð4aÞ 40 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 I @uu dXj @uh 8p3 lXj I ¼ al iu ih dS þ TM þ TL ; dt @r @r 3 ð4bÞ dqj ¼ Uj ; dt ð4cÞ where qj is the position vector of the jth particle center. A detailed description of the mechanical contact force, FM, the lubrication force, FL, and the corresponding torques TM and TL is given in the following sections. In order to evaluate the hydrodynamic force/torque integrals we need numerical estimates for the surface distribution of pressure, the radial derivatives of uh and u/. In the present implementation variables on the particle surface are evaluated on the isolatitude HEALPx spherical grids (Gorski et al., 2005). The derivatives on the particle surface are evaluated with third-order one-sided ﬁnite-differences utilizing three additional nodes on the radial axis passing through a given HEALPx node. The velocity at the auxiliary nodes on the radial axis is interpolated from 64 nearby Cartesian nodes using cubic Lagrange interpolation. The particle equations of motion (4) are reduced to an ODE system by evaluating the integrals with quadratures. The mechanical contact model used here is based on the work of Walton and Braun (1986); there exists extensive justiﬁcation for the model based on laboratory tests and ﬁnite-element simulations of colliding spheres and disks. Similar variations of the mechanical contact model used here have been previously used in discrete particle/element simulations of sediment transport processes (e.g. Drake and Calantoni, 2001; Calantoni and Thaxton, 2008). The mechanical contact model is a variant of the so-called ‘‘soft-sphere’’ models where the mechanical forces are functions of the overlap distance at the contact point. The normal force model is motivated by the experimentally observed hysteresis of the stiffness coefﬁcient where the stiffness under loading of the sample, k1, is smaller than the stiffness during unloading, k2. The overlap distance at the contact point between jth particle with radius, aj, and kth particle with radius, ak, is d ¼ maxf0; aj þ ak jqk qj jg: ð5aÞ Deﬁning the unit vector from the center of the jth to the kth particle as qk qj ; jqk qj j ð5bÞ the normal force on the jth particle is 8 > < k1 d Fn ¼ ejk k2 d ðk2 k1 Þdmax > : 0 ds ¼ Vjk ðVjk ejk Þejk ; dt ð6aÞ where Vjk ¼ Ujk 0:5ðaj þ ak dÞXjk ejk ; ð6bÞ is the total relative velocity at the contact point. The elastic tangential force on the jth particle, Ft, is obtained by time integration of dFt ds ¼ kt ; dt dt ð6cÞ assuming that Ft = 0 at the beginning of contact. Because the tangential plane may also vary with time, after each time step the result of (6c) must be projected back on the current tangential plane as Ft ðtÞ ¼ Ft ðFt ejk Þejk : ð6dÞ If the elastic force computed by ((6a)–(6d)) exceeds the friction limit then the tangential force is simply given by 2.3. Mechanical contact model ejk ¼ normal direction. The elastic tangential force (sticking) is motivated by experimentally observed reversal of the tangential velocity for certain impact conditions (see Section 3.3). Deﬁning the relative translational, Ujk = Uk Uj, and rotational velocities, Xjk = Xk + Xj, the rate of change of the tangential slip, s, at the contact point is ; dd=dt > 0 ; dd=dt < 0; d > dmax ð1 k1 =k2 Þ ; ; dd=dt < 0; d 6 dmax ð1 k1 =k2 Þ ð5cÞ where dmax is the maximum overlap at the end of the loading cycle. The above formulation assumes that the stiffness coefﬁcients are constant and that the unloading force reaches zero at ﬁnite overlap and remains zero upon further unloading. The formulation also accounts for plastic losses as the work during loading exceeds the work during unloading. We ﬁnally note that for constant k1 and pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ k2, the coefﬁcient of restitution, edry ¼ k1 =k2 , is also constant; this assumption is acceptable as long as the range of impact velocities is less than about 1 ms1 (see Walton, 1993). For the tangential component, we assume that the particles are either in sticking (with elastic tangential forces) or frictional-sliding contact (see Drake and Walton, 1995). The tangential stiffness coefﬁcient, kt, is constant and independent of the loading in the Ft ¼ lf jFn j Ft ; jFt j ð6eÞ where lf is the friction coefﬁcient. Finally the total mechanical contact force is FM ¼ Ft þ Fn ; ð7aÞ and the corresponding torque is TM ¼ aj FM ejk ; ð7bÞ 2.4. Lubrication and cavitation The hydrodynamic solver introduced in Section 2.1 can be used for fully resolved simulations of colliding particles by making the grid step of the same order as the particle surface roughness (a fraction of a micron). Unfortunately, such resolution is prohibitive for simulations with a large number of particles. Our goal is coarse grid simulations where the grid step is much larger that the surface roughness but is sufﬁciently small to resolve hydrodynamic effects on the particle scale (i.e. the outer solution). Therefore, for the purpose of such coarse-grid simulations it is necessary to introduce a model for the lubrication force when the interstitial gap becomes smaller than the grid step. The analytical models for the lubrication force considered here are based on the theoretical works of Cooley and O’Neill (1969) and O’Neill and Majumdar (1970). Cooley and O’Neill (1969) used the matched asymptotic expansion method to solve the Stokes ﬂow problem of two normally colliding spheres in the limit of vanishing clearance distance, ae, between the spheres, where e is a small dimensionless parameter. Their solution consists of an outer solution independent of e and a singular inner solution, which is conﬁned to a very small region along the gap having a lateral extent of about ae1/2. It was also found that the force (along the line of centers) on a sphere approaching an identical stationary sphere with velocity, U, is given by fL ¼ FL 1 9 lnðeÞ þ Oð1Þ; ¼ 6plaU 4e 40 ð8Þ where the outer solution contributed only to the O(1) term. Using the same approach, O’Neill and Majumdar (1970) determined the singular (/ ln(e)) tangential forces and torques for the case where the particle motion consists of relative translation along or rotation J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 about an axis that is perpendicular to the line of centers. Because of the linearity of the Stokes equations all these analytical models can be superimposed to predict the singular lubrication forces/torques for an arbitrary relative translation and rotation of two particles. Note, that we are excluding the case of relative rotation about the line of centers (considered by Jeffrey and Onishi (1984)), which results in a non-singular torque proportional to eln(e). We now present the general expressions for the forces and torques on the jth particle with translational and rotational velocity, Uj and Xj, due to a nearby kth particle having translational and rotational velocity, Uk and Xk, and radius, ajj, where j = ak/aj. The normal lubrication force on the jth particle due to relative translation along the line of centers (Cooley and O’Neill, 1969) is Fjn 6plaj j2 1 1 jð1 þ 7j þ j2 Þ e ln 2 eres 5ð1 þ jÞ3 ð1 þ jÞ e eres ¼ kðeÞ ðUjk ejk Þejk : ! ð9Þ Note that we have added terms that gradually increase the analytical force from zero at e = eres 3h/aj when the solution is fully resolved to approximately one as the clearance distance aje approaches the surface roughness length scale r ajeres. The factor kðeÞ is a step function weight, kðeÞ ¼ 1; eres > e > r=aj 0; otherwise introduced to keep the analytical force zero at large e and during contact. The above formulation includes the assumption that the particles make contact through micro-asperities of size, r, much smaller than the grid step. We also note that the numerically computed hydrodynamic force via the integral in (4a) represents the outer solution when e 6 eres (i.e. the O(1) term in (8)) and no weight factor is used for this term. The surface roughness parameter in (10) essentially controls the maximum value of the lubrication force (9) at contact; its effect on the collisional dynamics will be investigated in Section 3.2. Following O’Neill and Majumdar (1970), the tangential force and torque on the jth particle due to relative translation that is perpendicular to the line of centers are deﬁned as Fjtt 6plaj ¼ kðeÞ 4jð2 þ j þ 2j2 Þ 15ð1 þ jÞ 3 ½Ujk ðUjk ejk Þejk ln e ; eres ð11aÞ 8pla2j ¼ kðeÞ jð4 þ jÞ e ; ðejk Ujk Þ ln eres 10ð1 þ jÞ2 ð11bÞ where the weight function kðeÞ is used again to ensure that the analytical force does not contribute when the solution is fully resolved. Keeping in mind the deﬁnitions of Ujk and ejk from the previous section, (11a) and (11b) imply that Tktt ¼ Tjtt and Fktt ¼ Fjtt . For the case of relative rotation about an axis perpendicular to the line of centers (O’Neill and Majumdar, 1970), the tangential force and torque on the jth particle are given by Fjtr 6pl a2j ¼ kðeÞ 2j2 15ð1 þ jÞ 2 FjL ¼ Fjn þ Fjtt þ Fjtr ; ðXjk þ 4j1 Xj þ 4jXk Þ ejk ln e ; ð12aÞ eres and TjL ¼ Tjtt þ Tjtr : 8pl a3j ¼ kðeÞ 2j 5ð1 þ jÞ ð13bÞ The expressions for a particle interacting with a wall can be obtained from ((9), (11) and (12)) by taking the limit j ? 1. Here and throughout the remainder of the manuscript we only consider the case of identical sized particles, j = 1 (i.e. aj = ak = a). In the presence of nuclei such as microscopic bubbles of undissolved gas or gas pockets trapped on the particle surface, cavitation effects may inﬂuence the minimum pressure generated in the gap by the receding particle surfaces. Generally, the lubrication pressure will be minimum when the gap clearance is the smallest and will increase with increasing clearance. However, if cavitation is initiated, the minimum pressure in the gap p0 cannot drop below the vapor pressure, ð14Þ Essentially, cavitation may limit the magnitude of the negative dynamic pressure and the corresponding attractive lubrication force on rebound. To evaluate this effect, we assume that the cavitation effect is formally equivalent to a ‘‘roughness’’ distance, ecav, below which the attractive lubrication force is held constant at its maximum value. As a ﬁrst step in modeling the effect of cavitation, we propose to obtain an expression relating p0 and ecav using the well-known axisymmetric lubrication model for the radial ﬂow in the gap of two receding spheres in which the viscous shear of the radial ﬂow is maintained by an inward radial pressure gradient. The usual integration of this lubrication theory equation, double integration between the particle boundaries followed by an integration in the radial direction, yields the following result for the pressure, p0, in the center of the gap, p0 pouter ¼ 3lUjk ejk 1 ; 2a e2 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3Rec m=a pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ; e ¼ ecav ¼ 4 ðpouter pv Þ=q ð15Þ jXk jXk e : Xj þ Xj þ ejk ejk ln 4 4 eres ð12bÞ We ﬁrst note that the expression in the square brackets of (12b) does not have a component along the line of centers direction. As a consistency check, we also note that the total torque Ttr + Ttt ð16Þ where Rec = 2a(U2 U1)/m is the collision Reynolds number, computed at every time step for each potential collision pair. If the minimum pressure separation distance, aecav, is less than the surface roughness, r, then mechanical contact still controls the minimum pressure and the lubrication force model is unaffected by the cavitation effect. However, for aecav > r, cavitation may be accounted for by holding the attractive lubrication force constant at its maximum value Fjn 6pla and Tjtr ð13aÞ where pouter is the pressure at radial distances much greater than ae1/2. Eqs. (14) and (15) suggest that the cavitation pressure limit is reached when and Tjtt vanishes for a rigid-pair rotation where Xk = Xj = X, Uj = 0, and Ujk = Uk = aj(1 + j)X ejk. In our numerical implementation we have omitted the non-singular term eln(e) of (12a) which is negligible as e becomes very small. Combining (9), (11), and (12), the total lubrication force and torque on the jth particle are pv ¼ p0 : ð10Þ ; 41 ¼ 1 1 1 9 ecav ðUjk ejk Þejk ; ln 4 ecav eres 40 eres ð17Þ achieved at e = ecav; note, there would still be some weak dependence on e through changes in the relative velocity. We now compare aecav to realistic surface roughness scales r for low Rec collisions and natural water ﬂows (velocity less than 1 m/s) where the pressure outside particles will differ little from the ambient static pressure. For illustrative purposes we also assume shallow depths so that the hydrostatic pressure is much smaller than 42 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 1 atm and therefore pouter will not differ much from the astmospheric pressure pa. The water vapor pressure, pv, is a strong function of temperature that increases from 1 kPa to 100 kPa as the temperature increases from 0 to 100 °C. Thus, the largest vapor pressure is still two orders of magnitude smaller than the atmospheric pressure pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ and can safely be neglected in (16). Noting that pa =q ¼ 10 m s1 , 4 (16) suggests a cavitation scale, aecav ðRe1=2 mm or c Þ 10 aecav 1 lm for Rec = 100. Since realistic particle roughness is comparable to 1 lm cavitation effects will have little inﬂuence in the low Stokes number, St < 100, simulations considered here. Cavitation will also be negligible at larger St where hydrodynamic interactions do not seem to affect particle collisions. According to the present discussion, cavitation effects may only be important in the collisions of very smooth particles at low and intermediate St. 3. Results We demonstrated and veriﬁed our new hybrid model for collisions with available laboratory data. The total hydrodynamic force in the proposed hybrid collision model consists of a numerically determined contribution in the resolved ‘‘outer’’ region and a theoretical approximation for the force contribution due to the unresolved ‘‘inner’’ region (in the gap). The hybrid force model was tested by varying the resolution in a set of baseline simulations of normal collision showing that the dependence of the resolved hydrodynamic force on the interstitial distance matches the analytical dependence of e1. Simulations utilizing the optimal grid step h = a/16 were performed to validate the model prediction of experimentally determined coefﬁcients of restitution in normal collisions (YH06), and dry (FO94) and wet (YH06) oblique collisions. To facilitate comparison with the pendulum experiments of YH06, our numerical setup consists of a stationary target particle, not inﬂuenced by body forces, and a projectile that is accelerated by a body force. For simplicity, we assume that the body force acting on the projectile (representing the horizontal component of the tension on the pendulum string) is spatially uniform and acts solely in the y-direction in the simulation coordinate system. In all simulations, both the ﬂuid and the particles are initially at rest; the subsequent ﬂow is then purely the result of the projectile particle motion. The magnitude of the body force, g, for given ‘‘terminal’’ Reynolds number is determined from the balance of forces, mg ¼ qpm2 Re2 8 C D ðReÞ; ð18Þ where C D ðReÞ ¼ 24 ð1 þ 0:15Re0:687 Þ; Re ð19Þ is a standard empirical formula for the drag coefﬁcient of a sphere. Table 1 The main set of material parameters: normal stiffness k1, dry coefﬁcient of restitution edry, tangential stiffness kt, friction coefﬁcient lf, speciﬁc density qp/q, radius a, surface roughness r, vapor pressure pv and the corresponding cavitation length scale ecav. Note that friction and tangential stiffness affect only oblique collisions (Section 3.3). Also given are some references discussing the material parameters in relation to particle collisions. k1 (N m1) edry kt/k1 lf r (mm) a (mm) qp/q pv (kPa) ecav (mm) 104 0.97 0.7 0.25 0.0001 1.0 6.0 1.0 0.0001 Tsuji et al., 2008; Schafer et al., 1996 YH06; Drake and Walton, 1995 Drake and Walton, 1995 Drake and Walton, 1995 YH06 In all simulations, the ﬂuid had properties of water with density, q = 1 g cm3, and viscosity, m = 1 mm2 s1. The particles were spheres with radius, a = 1 mm, and have a large speciﬁc density, qp/q = 6, (e.g. steel) to minimize momentum losses during free particle motion. The material properties governing particle–particle interactions are given in Table 1. The physical dimensions of the computational domains were ﬁxed at an x–y–z size of 8 32 8 mm for all simulations. Periodic boundary conditions were used throughout all simulations to eliminate the inﬂuence of walls in domains of limited size. 3.1. Resolution tests We begin with a detailed analysis and comparison of the hybrid model for collision dynamics at three different grid resolutions, h = 1/16, 1/32, and 1/64, expressed as fractions of the particle radius, a. One of the goals of the resolution tests is to show that for h = 1/16 and Re < 100 the numerical model provides an acceptable ‘‘outer’’ solution even when the gap region is unresolved. The ability to fully resolve the ‘‘outer’’ dynamics is especially important during rebound when boundary layers may develop as the particle velocity changes direction over a very short period of time. The second goal is to demonstrate that our hybrid model gives a consistent prediction for the total hydrodynamic force as the gap hydrodynamics become unresolved. Both goals were achieved by comparing the partially resolved solution of the h = 1/16 case to the fully resolved solution of the h = 1/64 case (e.g. 4/64 < e < 4/16). For each different grid resolution a single baseline simulation was performed with the body force acceleration for the projectile having Cartesian components g = (0, 51.4, 0) mm s2, corresponding to a terminal Reynolds number, Re = 45. The actual impact Reynolds number (deﬁned below) will be somewhat smaller because the target is mobile and the collision takes place prior to reaching terminal velocity conditions. The initial Cartesian coordinates of the projectile and the target are (0, 6, 0) mm and (0, +6, 0) mm, respectively. In the low resolution case (h = 1/16) the computational grid has 128 512 128 nodes and the nominal time step was Dt = 0.0003 s; a smaller time step Dt = 5 ls was used to resolve the stiff mechanical response during contact. The collision kinematics of the particles when they are less than one diameter apart is shown in Fig. 1. During the initial acceleration stage not shown in Fig. 1, the projectile and the target have accelerated to speed U1 = 18 mm s1 and U2 = 1 mm s1, respectively; the pre-collision target motion is purely due to hydrodynamic interaction with the projectile. Here we deﬁne the impact Reynolds number, Rei, as the collision Reynolds number, Rec, computed at time, t = 0.01 s (time, t = 0 s, corresponds to maximum particle overlap during contact). Then, the impact Reynolds number is Rei = 34 and the corresponding impact Stokes number is Sti = (2qpRei)/9q = 45.3. Fig. 1 shows that for the relatively small normal stiffness, k1 = 104 N m1, the impact duration was short (100 ls) and resulted in extremely large accelerations. Even prior to contact the particle acceleration was several orders of magnitude larger than the imposed body force acceleration magnitude, g = 51.4 mm s1. The effect of these accelerations and the corresponding added mass forces on the particle dynamics is omitted in instantaneous collision models (hard-sphere); the question if this is justiﬁable will be investigated below. The total force on the projectile particle and its hydrodynamic and mechanical components are presented in Fig. 2. The robustness of the numerical method is demonstrated in Fig. 2a which shows little difference in the history of the total force in the three different resolution runs; note that, since the acceleration is proportional to the total force the particle kinematics was essentially the same in all three cases. We further note that the mechanical force (Fig. 2b thin) exceeded the total hydrodynamic force 43 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 (a) force [10−6 N] y [mm] 8 6 4 20 15 5 acceleration [m/s2] 0 (c) 1000 100 10 1 0 −1 −10 −100 (b) 10 force [10−6 N] (b) Uy [mm/s] 2 100 10 1 .1 .01 0 −.01 −.1 −1 −10 −100 (c) −0.1 −0.01 −0.001 0 0.001 0.01 0.1 force [10−6 N] (a) 1000 100 10 1 0 −1 −10 −100 1000 100 10 1 0 −1 −10 −100 time [s] −0.01 −0.001 0 0.001 0.01 0.1 time [s] Fig. 2. (a) Total force on the projectile for three different resolutions, h = 1/16 (thick black), h = 1/32 (thin black), and h = 1/64 (dashes), corresponding to the baseline simulation conditions, Re = 45. (b) Mechanical (thin) and hydrodynamic (thick) components of the total projectile force in the h = 1/16 simulation. (c) Components of the total hydrodynamic force in (b): analytical lubrication (thin), resolved viscous stress (dashes), and resolved pressure force (thick). 2500 h=1/16 h=1/16, gap unresolved h=1/32 h=1/32, gap unresolved h=1/64 h=1/64, gap unresolved 2000 x 10 5 2 1500 1.5 total force (Fig. 2b thick) by at least an order of magnitude. The total hydrodynamic force reversed its sign on contact and reached its largest value during contact. Interestingly, the components of the total hydrodynamic force in Fig. 2c showed that the hydrodynamic force maximum at contact was due to the resolved pressure. In the 4 ms interval prior to contact, the total hydrodynamic force was dominated by the singular nature of the analytical lubrication force (Fig. 2c, thin), but there was also a signiﬁcant contribution from the resolved pressure (Fig. 2c, thick), which opposes the lubrication force. A similar opposing pressure force also occurred during rebound. As will be discussed below, this large resolved pressure force was due to added mass effects associated with the large particle accelerations. The large accelerations result from mechanical forces during contact and lubrication forces immediately before and after contact. We ﬁnally note that the resolved pressure force (Fig. 2c, thick) signiﬁcantly exceeds any viscous drag contributions (Fig. 2c, dashes) in the 4 ms prior to contact and the 100 ls immediately after rebound. The nature of the large viscous forces that arose on rebound will be discussed in detail below. To examine the ability of our hybrid numerical–analytical model to make consistent predictions of the total hydrodynamic force even after the numerical solution becomes unresolved in the gap between approaching particles, we compared the asymptotic behavior of the total hydrodynamic force during approach as a function of inverse gap distance e1 at three different resolutions (Fig. 3). For e1 < 5 the numerical solution was fully resolved and the predicted force was purely numerical. In all three simulations there was a clear transition in the asymptotic behavior of the force such that the force became proportional to e1 for e < 1 (Fig. 3), a dependence expected from the analytical model (8) for the lubrication force. More importantly, Fig. 3 shows that for e1 > 12 when the h = 1/16 and h = 1/32 simulations were no longer fully resolved, −0.1 total hydrodynamic force [10−9 N] Fig. 1. Projectile (thin) and target (thick) position (a), velocity (b), and acceleration (c) in a low resolution simulation (h = 1/16) for terminal projectile Re = 45. To accommodate the extremely short collision time scale we use a special non-uniform time axes whose origin corresponds to maximum contact overlap and whose scaling is linear for 0.001 s < t < 0.001 s and logarithmic outside this interval. A similar scaling is used for the ordinate in (c) where the acceleration varies over 5 orders of magnitude. The transition from linear to logarithmic vertical scaling in (c) produces an artiﬁcial slope discontinuity (also Fig. 2). 1 1000 0.5 0 500 0 2 4 6 0 8 10 −1 ε 3000 1000 2000 ε−1 [mm−1] 12 14 16 18 20 −1 [mm ] Fig. 3. Total hydrodynamic force as a function of e1 during approach for three different grid resolutions, h = 1/16 (thick black), h = 1/32 (thin black), and h = 1/64 (thick gray), corresponding to the baseline simulation conditions, Re = 45. Solid lines denote fully resolved numerical solution; dashed lines denote unresolved gap dynamics. the difference between the predicted hybrid numerical–analytical force (thick and thin dashes) and the fully resolved h = 1/64 one (gray solid) does not exceed 5%. The inset of Fig. 3 shows that all three simulations remain very close and maintain the initial linear 44 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 (a) u(x,y,0) at ε=0.166 8 (b) h=1/16 h=1/64 7 7 6 6 5 5 y [mm] y [mm] h=1/16 v(x,y,0) at ε=0.166 8 4 3 15 contour interval=1 mm/s −2 −1 5 0 1 velocity [mm/s] h=1/16 h=1/64 10 0 −3 4 3 2 1 h=1/64 2 0 −3 0 3 2 contour interval=1 mm/s 1 −2 −1 0 1 2 3 x [mm] x [mm] (c) p(x,y,0) at ε=0.166 8 h=1/16 h=1/64 7 6 y [mm] 5 4 h=1/16 h=1/64 2 1 1500 1000 contour interval=50 mPa 500 0 −3 −2 −1 0 1 2 3 pressure [mPa] 3 0 x [mm] Fig. 4. Contoured is the velocity in the x-direction, u(x, y, 0) (a), the velocity in the y-direction, v(x, y, 0) (b) and pressure (c) at t = 0.01 s and e = 0.166 prior to maximum contact loading at low resolution, h = 1/16 (left panels) and high resolution, h = 1/64 (right panels). Positive and negative contour values are black and gray, respectively. The insets in (a) and (c) show the u velocity component and the pressure as a function of x along the middle of the gap. dependence down to very small gap clearances e = 3 104 comparable to the size of surface micro asperities. We now consider the ﬁdelity of the modeled outer solution in the interval, 0.01 s < t < 0.01 s, when the gap is under-resolved in the low resolution simulation. The hydrodynamics of approaching particles is illustrated by the velocity and pressure ﬁelds in Figs. 4 and 5 where the gap clearance is e = 0.166 and e = 0.063, respectively. As expected, the main feature consists of a viscous ﬂow driven out of the gap and the associated pressure maximum in the gap. Clearly, the low resolution case (Figs. 4c and 5c, left panels) underpredicts the maximum gap pressure compared to the fully resolved case (right panels). Nevertheless, the low and high resolution simulations differ very little in the predicted ﬂow and the pressure at distances jxj P a outside the gap. The ability to correctly model the outer solution when the solution in the gap is not resolved is further illustrated by the insets in Figs. 4 and 5 showing virtually no difference in the lateral velocity and pressure at distances greater than 1 mm from the center of the gap. The velocity and pressure ﬁelds during maximum particle overlap are shown in Fig. 6; the low and high resolution cases gave 45 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 (a) h=1/16 v(x,y,0) at ε=0.06 8 h=1/16 h=1/64 7 7 6 6 5 5 y [mm] y [mm] (b) u(x,y,0) at ε=0.06 8 4 h=1/64 4 3 3 1 20 15 contour interval=1 mm/s 10 5 0 −3 velocity [mm/s] 25 h=1/16 h=1/64 2 2 1 0 −2 −1 0 1 2 contour interval=1 mm/s 0 −3 3 −2 −1 0 1 2 3 x [mm] x [mm] (c) p(x,y,0) at ε=0.06 8 h=1/16 h=1/64 7 6 y [mm] 5 4 h=1/16 h=1/64 2 1 6000 4000 contour interval=50 mPa 2000 0 −3 pressure [mPa] 3 0 −2 −1 0 1 2 3 x [mm] Fig. 5. Contoured is the velocity in the x-direction, u(x, y, 0) (a), the velocity in the y-direction, v(x, y, 0) (b) and pressure (c) at e = 0.063 prior to maximum overlap at low resolution, h = 1/16 (left panels) and high resolution, h = 1/64 (right panels). Positive and negative contour values are black and gray, respectively. The insets in (a) and (c) show the u velocity component and the pressure as a function of x along the middle of the gap. essentially the same ‘‘outer’’ solution. Along the symmetry plane through the center of the gap the velocity in the x-direction, u(x, y, 0), was still outwards (Fig. 6a) but this outward ﬂow was surrounded by thin boundary layers ﬂowing in the opposite direction. The velocity in the y-direction, v(x, y, 0), at contact (Fig. 6b) was also similar to that at t = 0.01 s, except for thin boundary jets on the particle ﬂanks. It is interesting to note that the thin boundary layer ﬂows moved opposite to their attached particle and therefore do not seem to be generated by viscosity. However, the large shears associated with these thin jets must be responsible for the strong viscous drag force (Fig. 2c, dashes) that is observed in the interval 0 < t < 0.01 s after rebound. We should also point out that the maximum velocity of the jets around the projectile reached v(x, y, 0) = 24 mm s1, which was signiﬁcantly larger than the precollision velocity of the projectile Uy = 18 mm s1. Fig. 6c shows that the pressure ﬁeld was no longer dominated by the singular viscous dynamics in the gap but has a simpler dipolar structure and magnitude, which was 3 orders of magnitude larger than that at t = 0.01 s (Fig. 4c). This strong dipolar structure was consistent with the dipole source term in (3b) and the extremely large 46 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 (a) u(x,y,0) at ε=−0.0003 8 h=1/16 h=1/64 7 6 6 5 5 4 3 6 4 contour interval=1 mm/s −2 −1 0 1 2 3 velocity [mm/s] h=1/16 h=1/64 2 0 −3 h=1/64 4 3 2 1 v(x,y,0) at ε=−0.0003 8 7 y [mm] y [mm] h=1/16 (b) 2 contour interval=1 mm/s 1 0 0 −3 −2 −1 x [mm] 0 1 2 3 x [mm] (c) p(x,y,0) at ε=−0.0003 8 h=1/16 h=1/64 7 6 y [mm] 5 4 2 1 −100 −200 contour interval=50 mPa −300 h=1/16 h=1/64 0 −3 −2 −1 0 1 2 −400 pressure [Pa] 3 3 x [mm] Fig. 6. Contoured is the velocity in the x-direction, u(x, y, 0) (a), the velocity in the y-direction, v(x, y, 0) (b) and pressure (c) at maximum overlap, t = 0, at low resolution, h = 1/ 16 (left panels) and high resolution, h = 1/64 (right panels). Positive and negative contour values are black and gray, respectively. The insets in (a) and (c) show the u velocity component and the pressure as a function of x along the middle of the gap. particle accelerations (up to 100 ms2, see Fig. 1c) during contact. It is this inviscidly generated pressure and the associated pressure gradients that drive the thin boundary jets on the particle sides during contact (Fig. 6a and b). To further illustrate the added mass effect, we plotted the resolved pressure force on the projectile as a function of its acceleration (Fig. 7). When the particles are well separated (Fig. 7, close to the origin), the resolved pressure force increases with dUy/dt. However, when the particle acceleration becomes large (at small separations), the resolved pressure force becomes proportional to dUy/dt, which is consistent with an added mass force; the coefﬁcient of proportionality is approximately the same during approach (circles) and rebound (pluses). We also note that the dependence of the pressure force on the acceleration is not exactly linear because the added mass effect also depends on the gap distance. However, comparing our result to analytical predictions for the added mass forces between two accelerating particles is difﬁcult because the inviscid solution breaks down at very small gap distances. J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 Fig. 8 shows that the low-resolution baseline simulation was in good agreement with the high resolution case at time t = 0.01 s after the collision. Flow was being drawn inside the opening gap (Fig. 8a) while there was a pronounced velocity maximum behind the projectile (Fig. 8b) corresponding to the fact that the wake continued to move in the original direction. Although the projectile was nearly stationary the moving wake generated a pressure maximum behind the projectile that was as strong as that in front of the moving target. Given the above results we conclude that a grid step, h = 1/16, was sufﬁcient to adequately capture the hydrodynamic interaction during particle collisions at impact Stokes numbers, St < 100. The remainder of the results presented below will use the low resolution grid step, h = 1/16. 3.2. Coefﬁcient of normal restitution Here, we present the main validation of our hybrid model for the total hydrodynamic force by comparing the simulation results with laboratory experiments of normal collisions in water-glycerol solutions (YH06). The analytical lubrication correction to the hydrodynamic force is necessary because low-resolution simulations (even at 128 points per particle) do not capture the absence of collision rebound at low Stokes numbers. The loss of particle momentum in the normal direction, n e12, during collisions of two particles was measured by the normal coefﬁcient of restitution, e¼ ðU02 U01 Þ n ; ðU2 U1 Þ n ð20Þ where U1 and U2 are the incident velocities of the projectile and the 0 0 target, respectively, and U 1 and U 2 are the corresponding rebound velocities. For the purpose of comparison with the laboratory experiments of YH06, here we deﬁne the coefﬁcient of normal restitution based on the velocities at 0.01 s before and after contact. The particle material properties and the initial conditions were identical to the baseline simulations described above in Section 3.1, except that the magnitude of the body force, g, was varied (Table 2) to achieve different impact Stokes numbers, St. To also investigate how the dependence of e on St changes with the particle surface roughness, we performed two sets of simulations corresponding to r = 0.1 lm and r = 1 lm; the smaller surface roughness value is similar to the characteristic roughness of the 1 x 10 4 0.8 pressure force [10−9 N] 0.6 0.001 0.4 steel and glass particles used in YH06. Fig. 9 shows that for r = 0.1 lm (squares) our simulations predicted e as a function of St in good agreement with the laboratory experiments. More speciﬁcally, our model correctly predicted the observed critical Stokes number for rebound (St 3) without explicitly resolving the expensive lubrication hydrodynamics. The simulations with the larger particle roughness, r = 1 lm, (triangles) generally predicted a larger coefﬁcient of restitution when compared with the r = 0.1 lm simulations; such behavior was expected since the size of r limits the maximum lubrication force (9), (10). The effect of the surface roughness on the coefﬁcient of restitution, e, decreases with increasing Stokes number, St, which agrees with previous DNS of immersed particle collisions (Ardekani and Rangel, 2008). In the laboratory experiments (YH06), the effect of surface roughness generally appears as scatter in the measured value of e due to the random nature of asperity contacts in real particle collisions. 3.3. Oblique collisions In this section we compare our collision model against experimental data for tangential restitution in dry (FO94) and wet (YH06) oblique collisions. The collision experiments of FO94 veriﬁed a simpliﬁed (hard-sphere) version of Walton’s model where oblique collisions are fully determined by one friction coefﬁcient and two coefﬁcients for the restitution of the normal and tangential components of the particle relative velocity. It will be shown below that dry collision simulations utilizing Walton’s soft-sphere model (Section 2.3) also agreed with the experimental data of FO94. The effect of lubrication hydrodynamics on the coefﬁcient of tangential restitution will be inferred by comparing the dry collision results with simulations of collisions in water. In Walton’s hard-sphere model, the tangential contact changes from elastic (rolling) for collisions at low angles of incidence to frictional (sliding) at higher angles of incidence. At incident angles corresponding to the transition between these two regimes, the conversion of linear to angular momentum during dry oblique collisions is maximized and initially non-rotating particles rebound with the largest possible rotational rates. However, experiments of oblique collisions immersed in water-glycerol solutions (YH06) found very little tangential restitution and insigniﬁcant rebound rotation at the transitional angles of incidence. Motivated by the seeming discrepancy between the wet and dry collision experiments, we considered two sets of DNS, one with friction coefﬁcient, lf = 0.25, corresponding to the acetate particle experiments of FO94, and another with lf = 0.1 corresponding to the unpolished steel particle experiments of YH06. Following FO94 we examine the different oblique collision regimes (rolling and sliding) by plotting the recoil angle as a function of the angle of incidence. The angles of incidence and recoil are measured by their respective tangents, 0.002 0.2 0.02 .05 0.005 .1 .1 0.002 .05 0 W1 ¼ W2 ¼ −0.4 0.001 −0.6 −0.8 0 1000 2000 V12 t ; V12 n ð21Þ V012 t ; V12 n ð22Þ and −0.2 −1 −4000 −3000 −2000 −1000 47 3000 4000 dUy/dt [mm/s2] Fig. 7. The pressure force on the projectile as a function of the acceleration during approach (circles) and recoil (pluses) for the baseline simulation at low resolution, h = 1/16. The numerical labels correspond to the gap clearance e. where t is a unit vector in the direction of the tangential slip (6a), and V12 and V012 are the total incident and the total recoil relative velocity, respectively, at the contact point (6b). We also note that the ratio W2/W1 is the coefﬁcient of tangential restitution, which is generally a function of the angle of incidence (see below). For comparison with the dry collision experiments of FO94, we deﬁne the incident and the recoil relative velocity at t = ±0.001 s relative to the time of maximum overlap, while for wet collisions we follow YH06 and use the relative velocities at t = ±0.01 s. 48 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 (a) u(x,y,0) at ε=0.111 8 (b) h=1/64 h=1/16 7 7 6 6 5 5 y [mm] y [mm] h=1/16 4 4 contour interval=1 mm/s 2 velocity [mm/s] 6 h=1/16 h=1/64 2 0 −3 h=1/64 4 3 3 1 v(x,y,0) at ε=0.111 8 2 1 0 −3 0 −2 −1 0 1 2 contour interval=1 mm/s 3 −2 −1 0 1 2 3 x [mm] x [mm] (c) p(x,y,0) at ε=0.111 8 h=1/16 h=1/64 7 6 y [mm] 5 4 2 1 −200 −400 contour interval=20 mPa −600 h=1/16 h=1/64 0 −3 −2 −1 0 1 2 −800 pressure [mPa] 3 3 x [mm] Fig. 8. Contoured is the velocity in the x-direction, u(x, y, 0) (a), the velocity in the y-direction, v(x, y, 0) (b) and pressure (c) at t = 0.01 s after maximum contact loading, at low resolution, h = 1/16 (left panels) and high resolution, h = 1/64 (right panels). Positive and negative contour values are black and gray, respectively. The insets in (a) and (c) show the u velocity component and the pressure as a function of x along the middle of the gap. Table 2 Body force and the corresponding collision Stokes number for the simulations of normal collisions in water (Section 3.2). Terminal Re g (mm s2) Impact St 6 3.4 2.805 10 6.49 6.841 15 11.05 11.388 30 28.7 23.584 45 51.4 44.986 100 170.6 98.655 To isolate the mechanical contact model (Section 2.3) we ﬁrst considered oblique collisions of initially non-rotating particles in vacuum with lf = 0.25 by turning off the hydrodynamic model and the horizontal body force. The projectile was assigned an initial velocity, U1 = (0, 20, 0) mm s1, and the x-offsets, xoff, given in Table 3, while the target particle was motionless. The particle material parameters were again those in Table 1. We note that the dry coefﬁcient of restitution, edry = 0.87, for the cellulose acetate particles of FO94 is somewhat smaller than that in Table 1. Also, in our mechanical contact model there is no a priori deﬁned tangential restitution as the contact may switch from rolling to 49 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 1 (a) V [mm/s] steel (YH06) glass (YH06) 0.8 DNS σ=1μ DNS σ=0.1μ 20 15 10 5 0.6 0 e (b) U [mm/s] 0.4 0.2 1 2 10 10 (c) St Fig. 9. Shown is the coefﬁcient of normal restitution as a function of the impact Stokes number for simulated (triangles and squares) and measured (+ and ) particle collisions. Ω [rad/s] 0 0 10 6 4 2 0 −2 −4 −6 6 projectile target 4 2 0 −0.1 1.5 1 0 −0.5 0.5 1 −0.001 0 0.001 0.01 0.1 Fig. 11. Plotted are the translational (a and b) and rotational velocity (c) of the particles in the DNS of an oblique collision with initial offset, xoff = 1.3 mm, and friction coefﬁcient, lf = 0.25. Ψ2 0.5 0 −0.01 time [s] Foerster et al. (1994) DNS vacuum at ± 1ms DNS water at ± 1ms DNS water at ± 10ms 1.5 2 2.5 3 3.5 Ψ1 Fig. 10. Plotted is the tangent of the recoil angle, W2, as a function of the tangent of the incident angle, W1, from our DNS of oblique collisions of particles with friction coefﬁcient, lf = 0.25. The dashed line is a linear ﬁt (see text) to the experimental data of FO94. sliding during a single collision. Fig. 10 (triangles) shows that the predicted soft-sphere model dependence of W2 on W1 for collisions in vacuum follows closely the linear ﬁts (dashed lines) to the laboratory data (FO94, Fig. 5). The linear ﬁt for W1 < 1 corresponds to the rolling solution W2 = bW1 in Walton’s hard sphere model where b is the coefﬁcient of tangential restitution, and the ﬁt at W1 > 1 is the corresponding sliding solution, W2 = W1 7l(1 + e)/ 2 (see FO94). We next investigate the effect of hydrodynamic interactions in a series of DNS of particle collisions in water. The particle material properties were again those in Table 1 with the particle roughness, r = 0.1 lm. For the DNS we use the same setup as in the Section 3.1 (h = 1/16), consisting of a free initially motionless target and a projectile accelerated by a body force with Cartesian components g = (0, 51.4, 0) mm s1. However, here the initial projectile particle position was varied through a series of offsets in the x-direction, xoff (Table 3). The kinematics and hydrodynamics of the collision with xoff = 1.3 mm is illustrated in Figs. 11 and 12. Based on the relative velocity calculated at t = 0.01 s, the impact Reynolds number was Rei = 34 (Fig. 11a). The short collision duration of about 100 ls (Fig. 11) gives rise to very large accelerations (not shown) similar to those in Fig. 1. Even before contact, the target acquires a signiﬁcant y-velocity, 2 mm s1, due to strong normal lubrication forces. The tangential lubrication forces are much weaker and do not generate signiﬁcant rotation prior to mechanical contact (Fig. 11c). However, after the mechanical interaction the projectile and the target acquire relatively strong angular velocities of about 6 rad s1. The rotational motion also tends to die off at a faster rate (50% at t = 0.1 s) compared to particle motion in the transversal direction (Fig. 11b). The velocity and pressure ﬁeld sequence in Fig. 12 illustrate the overall hydrodynamic interaction of the particles before, during and after contact. Initially, when the particle separation is 1.5 mm, there is minimal interaction between the particles (Fig. 12a). As the gap separation decreases to 0.5 mm a pressure maximum is generated between the particles, which corresponds to the numerically resolved lubrication effect. During contact (Fig. 12c) the pressure ﬁeld around each particle is approximately dipolar and results from the large particle accelerations. After the collision at a separation distance of 0.5 mm (Fig. 12d), there is again little interaction between the particles and the Table 3 Lateral offsets for the oblique collision simulations. Projectile offset, xoff (mm) DNS vacuum, lf = 0.25 DNS water, lf = 0.25 DNS water, lf = 0.1 0.3 0.7 0.3 0.5 0.8 0.5 0.7 0.9 0.7 0.8 1.0 0.9 0.9 1.1 1.1 1.0 1.3 1.3 1.2 1.4 1.5 1.4 1.5 1.7 1.5 1.7 1.6 1.8 1.9 50 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 (a) (b) Pressure [Pa] −0.1 0 Pressure [Pa] 0.1 −0.1 8 7 7 6 6 y [mm] y [mm] 8 5 5 4 4 3 3 2 2 −1 0 1 2 1 −2 3 −1 Pressure [Pa] (c) −50 0 1 2 3 2 3 Pressure [Pa] (d) −0.1 50 0 0.1 9 9 19 mm/ s 17 mm/ s 8 8 7 7 6 6 y [mm] y [mm] 0 x [mm] x [mm] 5 5 4 4 3 3 2 2 1 −2 0.1 18 mm/ s 18 mm/ s 1 −2 0 9 9 −1 0 1 2 3 x [mm] 1 −2 −1 0 1 x [mm] Fig. 12. Velocity (arrows) and pressure (color) during approach (a and b), contact (c) and recoil (d) in the DNS of an oblique collision with initial offset, xoff = 1.3 mm, and friction coefﬁcient, lf = 0.25. For clarity, the velocity ﬁeld is plotted on every third/sixth grid node in the x/y-direction. pressure ﬁeld around the projectile is several times stronger since the particle continues to be accelerated due to the body force. The considered wet collision run for xoff = 1.3 mm resulted in the maximum rebound rotation rate and the minimum recoil tangent, W2 = 0.383 (evaluated at t = 0.001 s), among all wet collisions with lf = 0.25 (Fig. 10, squares). The effect of hydrodynamic interactions in oblique collisions seems to be twofold. First, the incident tangent, W1 = 1.14, for the above case was larger than the incident tangent, W1 = 0.97, for the dry collision with xoff = 1.4 mm. We attributed the increase to the repulsive hydrodynamic forces during approach that will tend to increase the collision obliqueness relative to dry collisions with the same initial 51 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 numerically resolved (pressure and viscous) drag to the integral terms above are denoted by IL and IN, respectively, and are shown in Table 4 for two of the normal collision simulations described in Section 3.2 corresponding to low and intermediate Stokes numbers, St = 6.8 and St = 45, respectively. We also estimated the contribution of the inviscid pressure drag, IP, to the total momentum change during mechanical contact as 2 1.5 Ψ2 1 R t2 IP ¼ 0.5 0 −0.5 YH06, unpolished steel DNS water at ± 1ms DNS water at ± 10ms 0 0.5 1 1.5 2 2.5 Ψ1 Fig. 13. Plotted is the tangent of the recoil angle, W2, as a function of the tangent of the incident angle, W1, from our DNS (squares, circles) and the corresponding laboratory data of YH06 for oblique collisions of particles with friction coefﬁcient, lf = 0.1. particle offset. Second, it appears that the recoil tangent W2 at t = 0.001 s in wet collisions (Fig. 10 squares) was consistently larger than that in dry collisions. To compare with the wet collision experiments of YH06 we consider another series of DNS with mechanical friction, lf = 0.1 (Fig. 13). The DNS model predictions for W1 and W2 estimated at t = ±0.01 s (circles) fell within (but on the lower side) of the experimental data (pluses). 4. Discussion The coefﬁcients of normal and tangential restitution have been subject to numerous laboratory studies as they are the key parameters in models for particle collision. An essential element to accurately predicting ﬂow kinematics (e.g. bedload sediment transport) is proper characterization of the energy dissipation due to grain–grain interactions. Much of the physics of particle interactions in DEM simulations for granular ﬂow relies on the speciﬁcation of the coefﬁcient of restitution. Consequently, our discussion is focused on the hydrodynamic effects that account for the differences between wet and dry coefﬁcients of restitution, the effect of speciﬁed particle stiffness in our simulations, and the speciﬁcation of the duration of collision for purposes of calculating the effective coefﬁcient of restitution in both laboratory measurements and our simulations. 4.1. Hydrodynamic effects It is instructive to analyze the contribution of various hydrodynamic effects to the momentum losses that account for the difference between the dry and the wet coefﬁcient of normal restitution. Integrating the difference of the momentum equations for particles 1 and 2 in the three intervals (t = 0.01 s, t1), (t1, t2) and (t2, t = 0.01 s), where t1 and t2 are the incident and recoil contact times, it can be shown that the dry and the wet coefﬁcient of restitution are related by ewet ¼ edry þ edry R t1 t¼0:01 s ðF 2 F 1 Þdt R t¼0:01 t2 mU 21 ðt ¼ 0:01 sÞ s ðF 2 F 1 Þdt ; ð23Þ where F is the total hydrodynamic force in the y-direction; note that the contributions from the constant body force, g, essentially cancel out in the two integrals because edry 1 and t1, t2 0.01 s. The individual contributions of the analytical lubrication force and the t1 ðF 2 F 1 Þdt mðU 21 ðt 2 Þ U 21 ðt1 ÞÞ ð24Þ : It is important to note that during mechanical contact the total hydrodynamic force is dominated by the pressure term. As expected, the main contribution to the momentum loss at low Stokes number, St = 6.8, is due to the contribution from the analytical lubrication force, IL = 0.966. The small positive IN (Table 4) results mainly from added mass forces that oppose the large decelerations produced by the lubrication force in the interval (t = 0.01 s, t1). However, at intermediate Stokes number, St = 45, the contribution from the resolved hydrodynamics, IN = 0.156, was comparable to that of the analytical correction, IL = 0.235. The main contribution to IN (0.130) took place in the recoil interval (t2, t = 0.01 s) and resulted mainly from viscous drag (see Fig. 2c). Although the shear in the gap was comparable to that on the particle sides (Fig. 8a and b) only the latter would give rise to viscous drag losses in the y-direction. We have previously mentioned that large shears were generated on the particle sides by the large added mass pressure during contact. At later times, when the particle accelerations were much smaller, the lateral shear was maintained by inertial effects such as wake interaction (cf. Ardekani and Rangel, 2008). We ﬁnally note that at low and intermediate St the added mass effect, IP, contributed to 13% of the momentum change during mechanical contact. Therefore, the added mass force was not negligible when compared to mechanical forces (see also Fig. 2). On the other hand, the coefﬁcient of restitution estimate, econ, based on the contact velocities at t = t1, t2 differs little from the dry value, edry = 0.97 (Table 4). There is a simple explanation why the added mass force does not affect the coefﬁcient of restitution during contact. Being proportional to the particle acceleration, the resolved pressure force during contact only modiﬁes the effective particle mass and the loading/unloading particle interaction is controlled solely by the stiffness coefﬁcients of the prescribed mechanical force. It was suggested previously (Sangani and Didwania, 1993) that added mass effects associated with the collision of two bubbles in a bubbly ﬂow would affect the velocity of all other bubbles. In this respect, it is interesting to investigate the inﬂuence of a nearby third particle on the collision dynamics of the two particles in the baseline St = 45 simulation (Fig. 1). For this purpose, the baseline simulation was repeated by including a third particle one diameter away from the target, at y3 = 6 mm and x3 = ±4 mm; note that the third particle lies on the periodic boundary at x = ±4 mm. The effect of the third particle was negligible until t = 0.01 s prior to contact when the velocity of the projectile and the target differed only in the fourth signiﬁcant digit from the corresponding velocities in the original simulation. The presence of the third particle also had a minor effect on the collision as the effective coefﬁcient of restitution e = 0.550 was only 5% smaller that the one in the two-particle Table 4 The contribution of resolved, IN, and unresolved, IL, hydrodynamic losses to the overall coefﬁcient of normal restitution, e = edry + IN + IL. Also given are the resolved hydrodynamic contribution, IP, to the momentum change and the actual simulation coefﬁcient of restitution, econ, during contact. Impact St e IL IN IP econ 44.986 6.841 0.579 0.081 0.235 0.966 0.156 +0.077 0.131 0.133 0.967 0.969 52 J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 simulation. Like the target particle, the third particle begins to move in the y > 0 direction prior to collision, but at t = 0.01 s its velocity U3 = 0.08 mm s1 was an order of magnitude smaller compared to that of the target. During mechanical contact of the other two particles, the third particle acquired relatively large acceleration due to the large added mass pressure, which reduced the velocity to U3 = 0.05 mm s1. In a similar three-particle simulation where x3 = 2.1 mm and y3 = 5.0 mm, the third particle acquired a larger velocity U3 = 1 mm s1 at t = 0.01 s that was subsequently reduced by 25% due to the added mass efects; however, the effective coefﬁcient of restitution e = 0.546 was again not signiﬁcantly different from that in the baseline two-particle simulation. In the above three-particle simulations, the third particle and the projectile remain well separated (e > 0.2) and do not interact through strong lubrication forces. A larger reduction of the effective coefﬁcient of restitution might be expected in certain collisions where the initial position of the third particle results in a much stronger interaction with the projectile before the latter collides with the target particle. 4.2. Particle stiffness In the mechanical contact model used here the contact time is expected pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ to scale as the period of a simple harmonic oscillator, m=k1 . For realistic stiffness values the contact interval is extremely small and represents an important limitation for large-scale simulations. Contact is detected after the fact, i.e. when d ﬁrst becomes non-zero, thus the accuracy of contact detection and the integration of mechanical forces depends critically on the time step being sufﬁciently small. The value of the normal stiffness, k1 = 104 N m1, used for the collisions in Fig. 9 is rather low compared to the stiffness of actual materials in particle-laden ﬂow applications. For example, the normal stiffness used in Section 3 is signiﬁcantly lower than measured values for glass spheres (Drake and Walton, 1995) or estimates, Ea, based on the Young’s modulus, E, (typically 10–100 GPa). Previous work has shown that the dynamics of granular ﬂows is rather insensitive to the value of normal stiffness across a range of two orders of magnitude from about 103–105 N m1 (e.g. Jiang and Haff, 1993). However, both the type of contact law (e.g. Hertzian versus Hookean) and the value of the stiffness are critical for problems such as acoustic wave propagation in granular media where the wave speed has been shown to be proportional to the stiffness of interparticle contacts (e.g. Sadd et al., 2000). To evaluate the effect of the stiffness parameter on particle collisions, here we repeat the St = 45 normal collision simulation of Section 3.2 with several different values of k1 (Table 5). The results showed that there was less than 4% difference in the predicted coefﬁcient of restitution for stiffness values in the range 103 < k1 < 106 N m1. Therefore using a lower stiffness, k1 = 104 N m1, to reduce computing time should predict rebound velocities consistent with simulations with realistic values of k1. 4.3. Duration of collision Due to technical limitations, such as high-speed camera frame rates, it is generally not possible to measure the particle velocities immediately before and after contact. Typically, the measured incident and rebound velocities are those at about t = ±0.01 s before and after contact (YH06). For collisions in air (dry collisions) and typical laboratory size particles the hydrodynamic losses are Table 5 The coefﬁcient of normal restitution for impact St = 45, edry = 0.97 and different particle stiffness. k1 (N m1) 103 104 105 106 e 0.567 0.579 0.561 0.570 relatively small and the velocities at t = ±0.01 s are not much different from those just before and after contact. However, in wet collisions, a coefﬁcient of normal restitution based on the velocities at t = ±0.01 s is generally smaller than the corresponding dry coefﬁcient due to signiﬁcant hydrodynamic losses (Fig. 9). Likewise, the estimates of W1 and W2 will be insensitive to the particular choice of incident and rebound velocity for dry collisions, but will depend strongly on the choice of incident and recoil velocity in wet collisions. The latter is illustrated in Figs. 10 and 13 where the recoil tangent estimates based on impact and rebound velocities at t = ±0.01 s (circles) generally lie above the corresponding estimates based on the velocities at t = ±0.001 s (squares). On the other hand, when the impact and rebound velocities are those at t = ±0.0001 s the corresponding W1 and W2 (not shown) are essentially the same as in the dry simulations. For example, using the t = ±0.0001 s velocities, our simulations predict minimum W2 = 0.185 (not shown) that is comparable to the lf = 0.1 collision experiments (FO94, Fig. 4) where a minimum W2 = 0.2 is located at W1 = 0.4 0.5. The above results clearly illustrate that a difference between the dry and wet coefﬁcient of tangential restitution, W2/W1, only arises when hydrodynamic effects are incorporated by choosing ‘‘impact’’ and ‘‘recoil’’ velocities at some ﬁnite time interval before and after mechanical contact. Given the ambiguity in the choice of this time interval it seems more appropriate to treat hydrodynamic particle interactions separately rather than including them in the coefﬁcients of normal and tangential restitution. 5. Conclusions We developed a model for immersed particle collisions by coupling a numerical model for DNS of particle hydrodynamics with analytical models for mechanical contact and unresolved microscale lubrication forces. We also made a simple theoretical evaluation of cavitation effects on the lubrication forces during recoil and concluded that cavitation may affect only the collisions of extremely smooth particles. The consistency of the coupled numerical–analytical force model was veriﬁed in resolution tests. The model ability to simulate the observed dependence of the coefﬁcient of normal restitution on the Stokes number in wet collisions (YH06) was established in Section 3.2 and good agreement was found between the predicted and the observed critical Stokes number for collision without rebound. The increase of the predicted coefﬁcient of normal restitution with the particle roughness was consistent with previous simulations (Ardekani and Rangel, 2008). Our model predictions for tangential restitution, W2/W1, as a function of the angle of incidence in oblique collisions were veriﬁed in a comparison with experimental data for dry collisions with friction coefﬁcient, lf = 0.25 (FO94), and wet collisions with friction coefﬁcient, lf = 0.1 (YH06). Our oblique collision simulations showed that in wet collisions the coefﬁcient of tangential restitution depends strongly on the choice of incident and rebound velocity. It was also shown that in the limit when the incident and rebound velocity are those immediately before and after contact, the tangential coefﬁcient of restitution in wet collisions asymptotes to the corresponding dry values. Considering the ambiguity in choosing impact and recoil velocities other than those immediately before and after contact, we suggest that models of immersed particle collision should not include the hydrodynamic effects in the coefﬁcients of normal and tangential restitution. The latter should be used to describe only the mechanical interaction and the hydrodynamic interaction should be treated with force models such as those in Section 2.4. In relation to normal collisions, we analyzed the contribution of different hydrodynamic effects to momentum dissipation and J.A. Simeonov, J. Calantoni / International Journal of Multiphase Flow 46 (2012) 38–53 found that at low Stokes numbers the lubrication effect was responsible for most of the dissipation. However, at intermediate Stokes numbers viscous effects associated with particle–particle and particle–wake interaction were also found to be important. We also found that the strong added mass forces associated with large collisional accelerations did not affect the coefﬁcient of restitution during mechanical contact. While this opens the possibility to use hard-sphere collision models in DNS of particle-laden ﬂow (cf. Ardekani and Rangel, 2008), such models will not be useful in applications that involve enduring contact (e.g. bedload sediment transport). Acknowledgements J.A.S. was supported through the ASEE-NRL postdoctoral fellowship program and the Jerome and Isabella Karle Distinguished Scholar Fellowship Program at the Naval Research Laboratory. J.C. was supported under base funding to the Naval Research Laboratory from the Ofﬁce of Naval Research (PE#61153 N). This work was supported in part by a grant of computer time from the DoD High Performance Computing Modernization Program at the NAVY DSRC and the ERDC DSRC. We thank the anonymous reviewer for helpful comments. References Ardekani, A.M., Rangel, R.H., 2008. Numerical investigation of particle-particle and particle-wall collisions in a viscous ﬂuid. J. Fluid Mech. 596, 437–466. Barnocky, G., Davis, R.H., 1988. Elastohydrodynamic collision and rebound of spheres: experimental veriﬁcation. Phys. Fluids 31 (6), 1324–1329. Calantoni, J., Thaxton, C.S., 2008. Simple power law for transport ratio with bimodal distributions of coarse sediments under waves. J. Geophys. Res. 113, C03003. Cooley, M.D.A., O’Neill, M.E., 1969. On the slow motion generated in a viscous ﬂuid by the approach of a sphere to a plane wall or stationary sphere. Mathematika 16, 37–49. Dance, S.L., Maxey, M.R., 2003. Incorporation of lubrication effects into the forcecoupling method for particulate two-phase ﬂow. J. Comput. Phys. 189, 212–238. Davis, R.H., Serayssol, J.-M., Hinch, E.J., 1986. The elastohydrodynamic collision of two spheres. J. Fluid Mech. 163, 479–497. Drake, T.G., Walton, O.R., 1995. Comparison of experimental and simulated grain ﬂows. J. Appl. Mech. 62, 131–135. Drake, T.G., Calantoni, J., 2001. Discrete particle model for sheet ﬂow sediment transport in the nearshore. J. Geophys. Res. 106, 19859–19868. 53 Foerster, S.F., Louge, M.Y., Chang, H., Allia, K., 1994. Measurements of the collision properties of small spheres. Phys. Fluids 6 (3), 1108–1115. Glowinski, R., Pan, T.-W., Hesla, T.I., Joseph, D.D., 1999. A distributed Lagrange multiplier/ﬁctitious domain method for particulate ﬂows. Int. J. Multiphase Flow 25, 755–794. Gorski, K.M., Hivon, E., Banday, A.J., Wandelt, B.D., Hansen, F.K., Reinecke, M., Bartelmann, M., 2005. HEALPix: a framework for high-resolution discretization and fast analysis of data distributed on the sphere. Astrophys. J. 622, 759–771. Hu, H.H., 1996. Direct simulations of ﬂows of solid-liquid mixtures. Int. J. Multiphase Flow 22 (1996), 335–352. Hu, H.H., Patankar, N.A., Zhu, M.Y., 2001. Direct Numerical Simulations of ﬂuid-solid systems using the Arbitrary Lagrangian-Eulerian technique. J. Comput. Phys. 169, 427–462. Jeffrey, D.J., Onishi, Y., 1984. The forces and couples acting on two nearly touching spheres in low-Reynolds-number ﬂow, J. Appl. Math. Phys. (ZAMP) 35, 634– 641. Jiang, Z., Haff, P.K., 1993. Multiparticle simulation methods applied to the micromechanics of bed-load transport. Water Resourc. Res. 29 (2), 399–412. Johnson, A.A., Tezduyar, T.E., 1996. Simulation of multiple spheres falling in a liquid-ﬁlled tube. Comp. Meth. Appl. Mech. Eng. 134, 351–373. Maury, B., 1997. A many-body lubrication model,. C.R. Acad. Sci. Paris 325 (Ser. I), 1053–1058. Mayo, A., 1984. The fast solution of Poisson’s and the biharmonic equations on irregular regions. SIAM J. Numer. Anal. 21, 285–299. O’Neill, M.E., Majumdar, S.R., 1970. Asymmetrical slow viscous ﬂuid motions caused by the translation or rotation of two spheres. Part II: asymptotic forms of the solutions when the minimum clearance between the spheres approaches zero. Z. Angew. Math Phys. 21, 180–187. Sadd, M.H., Adhikari, G., Cardoso, F., 2000. DEM simulation of wave propagation in granular materials. Powder Technol. 109 (1–3), 222–233. Sangani, A.S., Didwania, A.K., 1993. Dynamic simulations of ﬂows of bubby liquids at large Reynolds-numbers. J. Fluid. Mech. 250, 307–337. Schafer, J., Dippel, S., Wolf, D.E., 1996. Force schemes in simulations of granular materials. J. Phys. I 6 (1), 5–20. Simeonov, J.A., Calantoni, J., 2011. A pressure boundary integral method for direct ﬂuid-particle simulations on Cartesian grids. J. Comp. Phys. 230 (5), 1749–1765. Singh, P., Hesla, T.I., Joseph, D.D., 2003. Distributed Lagrange multiplier method for particulate ﬂow with collisions. Int. J. Multiphase Flow 29, 495–509. Tsuji, T., Yabumoto, K., Tanaka, T., 2008. Spontaneous structures in threedimensional bubbling gas-ﬂuidized bed by parallel DEM–CFD coupling simulation. Powder Technol. 184 (2), 132–140. Walton, O.R., Braun, R.L., 1986. Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J. Rheol. 30 (5), 949–980. Walton, O.R., 1993. Numerical simulation of inelastic, frictional particle-particle interactions. In: Roco, M.C. (Ed.), Particulate Two-Phase Flow. ButterworthHeineman, Boston, pp. 84–911. Yang, F.-L., Hunt, M.L., 2006. Dynamics of particle–particle collisions in a viscous liquid. Phys. Fluids 18, 121506.