...

Simeonov2012.pdf

by user

on
Category: Documents
1

views

Report

Comments

Transcript

Simeonov2012.pdf
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 flow
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 flow 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 fluid–particle simulations on Cartesian grids. The unresolved
part of the lubrication pressure/shear force in the subgrid gap was predicted using theoretical Stokes flow
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 first objective is to include a soft-sphere model (Walton and
Braun, 1986) for mechanical contact in Direct Numerical Simulations (DNS) of finite 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
modified 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 coefficient 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.ijmultiphaseflow.2012.05.008
restitution (Hu et al., 2001; Ardekani and Rangel, 2008); however,
a rotational coefficient 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 flow
which has a tangential velocity discontinuity at the boundary;
the initially infinite boundary layer gradients correspond to a singular viscous force which cannot be properly resolved with a finite
numerical discretization. Second, the infinite particle acceleration
during such instantaneous collisions corresponds to an infinite
added mass drag (of vanishing duration) and the corresponding
fluid impulse is indeterminable. Collisional added mass effects
apparently are not accounted for in hard-sphere models but may
be important. The significance 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 fluid–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 efficiency 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 flow in the gap could be resolved explicitly by using
very fine 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-fitted 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 fluid–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 findings 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 fluid flow and pressure generated by
the motion of spherical particles of radius, a, over a regular
Cartesian grid. The fluid 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 field across the particle boundaries. Noting that there is
complete freedom in the manner of specifying the fictitious interior pressure, we make the simplest choice where the interior pressure field 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 specified 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 finite difference Laplacian
(so that second order accuracy is preserved across discontinuities)
for stencils intersecting particle boundaries and the resulting finite-difference Poisson equation is solved with fast algorithms on
the Cartesian grid. This completes the solution for the pressure
field on the Cartesian grid and allows computing the pressure forcing terms in the fluid (2a) and the particle momentum equations
(below). The fluid momentum Eq. (2a) is reduced to an ODE system
using second order finite 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 efficiency 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 finite-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 justification
for the model based on laboratory tests and finite-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 coefficient 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Þ
Defining 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).
Defining 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 coefficients are
constant and that the unloading force reaches zero at finite 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 finally note
that for constant k1 and
pffiffiffiffiffiffiffiffiffiffiffiffi
k2, the coefficient 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
coefficient, kt, is constant and independent of the loading in the
Ft ¼ lf jFn j
Ft
;
jFt j
ð6eÞ
where lf is the friction coefficient. 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 sufficiently 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 flow 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 confined 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 defined 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 definitions 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 influence 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 first 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 flow
in the gap of two receding spheres in which the viscous shear of
the radial flow 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
rffiffiffiffiffiffiffiffiffiffi
3Rec
m=a
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ;
e ¼ ecav ¼
4
ðpouter pv Þ=q
ð15Þ
jXk
jXk
e
:
Xj þ
Xj þ
ejk ejk ln
4
4
eres
ð12bÞ
We first 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 flows (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
pffiffiffiffiffiffiffiffiffiffiffi
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 influence 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 verified 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 coefficients 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 influenced 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 fluid and the particles are initially at rest;
the subsequent flow 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 coefficient of a sphere.
Table 1
The main set of material parameters: normal stiffness k1, dry coefficient of restitution
edry, tangential stiffness kt, friction coefficient lf, specific 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 fluid 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 specific 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 fixed 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 influence
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 (defined 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 define 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 justifiable 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 significant 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 finally note that the resolved pressure force
(Fig. 2c, thick) significantly 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 artificial 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 fidelity 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 fields 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
flow 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 flow
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 fields 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 flow was
surrounded by thin boundary layers flowing 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 flanks. It is interesting to note that the thin boundary layer flows 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 significantly larger than
the precollision velocity of the projectile Uy = 18 mm s1. Fig. 6c
shows that the pressure field 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 coefficient 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 difficult 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 sufficient 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. Coefficient 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 coefficient 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 define the coefficient 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 specifically, 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 coefficient 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 coefficient 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 verified a simplified (hard-sphere) version of Walton’s model where
oblique collisions are fully determined by one friction coefficient
and two coefficients 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 coefficient
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 insignificant 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 coefficient, 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 coefficient 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 define
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 first
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 coefficient 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 defined
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 coefficient 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 coefficient, 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
coefficient, lf = 0.25. The dashed line is a linear fit (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 fits (dashed lines) to the laboratory data (FO94, Fig. 5). The linear fit for W1 < 1 corresponds to
the rolling solution W2 = bW1 in Walton’s hard sphere model
where b is the coefficient of tangential restitution, and the fit 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 significant y-velocity, 2 mm s1, due to strong normal lubrication
forces. The tangential lubrication forces are much weaker and do
not generate significant 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 field 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 field 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 coefficient, lf = 0.25. For clarity, the velocity field is plotted on every third/sixth grid node in the x/y-direction.
pressure field 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 coefficient,
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 coefficients 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 flow 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 flow relies on the
specification of the coefficient of restitution. Consequently, our discussion is focused on the hydrodynamic effects that account for the
differences between wet and dry coefficients of restitution, the effect of specified particle stiffness in our simulations, and the specification of the duration of collision for purposes of calculating the
effective coefficient 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 coefficient 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 coefficient 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 finally 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 coefficient 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 coefficient of restitution during contact. Being
proportional to the particle acceleration, the resolved pressure
force during contact only modifies the effective particle mass and
the loading/unloading particle interaction is controlled solely by
the stiffness coefficients 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 flow would affect the velocity of all other bubbles. In this
respect, it is interesting to investigate the influence 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 significant 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 coefficient 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
coefficient of normal restitution, e = edry + IN + IL. Also given are the resolved hydrodynamic contribution, IP, to the momentum change and the actual simulation
coefficient 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 coefficient of restitution e = 0.546 was again not significantly 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 coefficient 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
pffiffiffiffiffiffiffiffiffiffiffiffi 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 first becomes non-zero, thus the accuracy of contact detection and the
integration of mechanical forces depends critically on the time step
being sufficiently 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 flow
applications. For example, the normal stiffness used in Section 3
is significantly 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 flows 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 coefficient 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 coefficient 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 coefficient of normal restitution based on the velocities at
t = ±0.01 s is generally smaller than the corresponding dry coefficient due to significant 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 coefficient of tangential restitution,
W2/W1, only arises when hydrodynamic effects are incorporated
by choosing ‘‘impact’’ and ‘‘recoil’’ velocities at some finite 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 coefficients 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 verified in resolution tests. The
model ability to simulate the observed dependence of the coefficient 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
coefficient 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
verified in a comparison with experimental data for dry collisions
with friction coefficient, lf = 0.25 (FO94), and wet collisions with
friction coefficient, lf = 0.1 (YH06). Our oblique collision simulations showed that in wet collisions the coefficient 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 coefficient 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 coefficients 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 coefficient of restitution during mechanical contact. While this opens the possibility
to use hard-sphere collision models in DNS of particle-laden flow
(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 Office 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 fluid. J. Fluid Mech. 596, 437–466.
Barnocky, G., Davis, R.H., 1988. Elastohydrodynamic collision and rebound of
spheres: experimental verification. 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 fluid
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 flow. 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
flows. J. Appl. Mech. 62, 131–135.
Drake, T.G., Calantoni, J., 2001. Discrete particle model for sheet flow 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/fictitious domain method for particulate flows. 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 flows 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 fluid-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 flow, 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-filled 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 fluid 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 flows 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
fluid-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 flow with collisions. Int. J. Multiphase Flow 29, 495–509.
Tsuji, T., Yabumoto, K., Tanaka, T., 2008. Spontaneous structures in threedimensional bubbling gas-fluidized 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.
Fly UP