Collision dynamics of granular particles with adhesion

We investigate the collision of adhesive viscoelastic spheres in quasistatic approximation where the adhesive interaction is described by the Johnson, Kendall, and Roberts JKR theory. The collision dynamics, based on the dynamic contact force, describes both restitutive collisions quantified by the coefficient of restitution as well as aggregative collisions, characterized by the critical aggregative impact velocity gcr. Both quantities and gcr depend sensitively on the impact velocity and particle size. Our results agree well with laboratory experiments.


I. INTRODUCTION
Granular systems are abundant in nature and examples range from sands and powders on Earth ͓1,2͔ to more dilute systems, termed as granular gases ͓3-5͔, in space; planetary rings and other astrophysical objects may be mentioned in this respect ͑see, e.g., ͓3,6͔͒.Granular matter exhibits highly unusual properties: In a dense state it behaves similar to a liquid or solid, while in a rarefied state it resembles in some respect common molecular gases.This rich behavior stems not only from the density variation and related lack of scale separation, but to a large extent from specific particle interactions ͓2͔.The interaction force is a combination of elastic rebound, dissipation due to viscous deformations, and adhesion caused by the molecular van der Waals forces.The interplay between these three basic contributions leads to a rich collision behavior of adhesive, dissipative particles, which in turn determines the unusual macroscopic properties of granular systems.
In the present study, we address the collision dynamics of dissipative particles with adhesive interactions.Particles in real granular systems may have a complicated nonspherical shape and differ in mass, size, and composition.Here we consider a simplified case where all grains are smooth spheres of the same material.Moreover, we consider exclusively mechanical interactions thereby disregarding longrange forces such as electrostatic forces or gravity.
We apply the viscoelastic collision model ͓7,8͔ and incorporate adhesive interactions based on the Johnson, Kendall, and Roberts theory ͑JKR͒ ͓9͔ which was shown to be appropriate for systems of practical interest ͓10,11͔.To justify the assumption of quasistatic deformations for the viscoelastic model we assume that the impact speed is small compared to the speed of sound.In contrast to numerous studies of viscoelasticity with adhesion for elastomeric materials ͑e.g., ͓12-14͔͒ where the memory of the viscous material was described by relaxation ͑or creep͒ functions, we neglect memory effects.That is, we assume that the characteristic times of the collision process is much larger than the dissipative relaxation time of the particles' material.This assumption allows one to express the dissipative part of the stress through the elastic stress and to develop an analytical theory.We derive an expression for the total force acting between colliding spheres and demonstrate that it may not be written as a simple superposition of elastic, adhesive, and viscous forces; instead, it includes a viscoadhesive cross term.We study numerically and analytically the collision dynamics and compute the coefficient of restitution as a function of impact rate and particle size.Finally, we compute the maximal impact speed at which aggregation may occur, separating the domain of restitutive collisions where particles rebound, from the domain of aggregative collisions where particles constitute a joint aggregate.

II. PARTICLE INTERACTION FORCES A. Elastic force
To introduce the notation first let us consider the contact of elastic spheres of radii R 1 and R 2 that are made of the same material.Their material properties are described by a single parameter, D ϵ͑3/2͒͑1− 2 ͒ / Y, where Y and are Young modulus and Poisson ratio, respectively.The compression ϵ R 1 + R 2 − ͉r ជ 1 − r ជ 2 ͉ characterizes the deformation of and a the radius of the contact area between the particles.
The normal pressure distribution at the contact plane is a function of the compression ͓15͔ where x 2 + y 2 ഛ a 2 and x = y = z = 0 denotes the center of the circular contact area.Stress ij and strain u ij are linearly related ͓16͔ where Here, u ជ͑r ជ͒ is the displacement field and the indices i , j , l denote Cartesian coordinates ͑Einstein's summation rule ap-plies͒.The elastic coefficients E 1 and E 2 read The repulsive force F H is obtained by integrating the stress distribution over the contact region and reads where The compression depends linearly on the contact area, a 2 = R eff .

B. Elastic adhesive force
At given compression the contact area of adhesive spheres is larger than the contact area of elastic spheres as predicted by Hertz's theory.Among several theories describing adhesive contact ͑e.g., ͓9-11,13,17-28͔͒ the JKR theory was shown to be accurate enough for a wide range of applications of practical interest ͑see, e.g., ͓10͔͒ and is, moreover, suited for an analytical analysis due to its simplicity.The applicability of the JKR theory is characterized by the value of the Tabor parameter,

͑7͒
where z 0 is the characteristic atomic scale ͓21͔.As a rule of thumb, the JKR theory is reliable for T տ 5 ͓10,11͔, otherwise the DMT theory ͓17͔ is preferable.Throughout this article we refer to applications where T Ͼ 5.
The JKR theory combines two basic solutions for the pressure distribution of an elastic contact, the Hertz solution ͓15͔ for a compressed sphere, Eq. ͑1͒, and the Boussinesq solution for the uniform displacement of a circular area in a plane, oriented normally to the surface ͑e.g., ͓18,29͔͒ where the Boussinesq force F B depends on the corresponding displacement B and the contact radius a ͓9͔ In the JKR theory an effective Hertzian force is introduced which yields the observed contact of radius a and would lead to the compression H in a lack of adhesion.The actual compression is, however, smaller than in the purely elastic case, Ͻ H .The difference H − is attributed to the Boussinesq force F B and the corresponding displacement B , which is related to the contact radius as where the adhesive coefficient ␥ is twice the surface free energy per unit area of a solid in vacuum.F B acts against the Hertzian force and, thus, reduces the compression.Therefore, the force between adhesive spheres is also reduced so that the total force and compression read where Eqs.͑5͒, ͑9͒, and ͑10͒ were used.Correspondingly, the total ͑static͒ stress ij st ͓which still obeys the general form of Eq. ͑3͔͒, is now a sum of two components, due to the Hertzian and Boussinesq's displacements in the material, u ជ H and u ជ B .From Eqs. ͑11͒ and ͑12͒ we find the contact radius and compression a eq for unforced adhesive spheres ͑F =0͒, If the particles in contact are pulled apart from one another the radius of the contact area a decreases.As a becomes smaller than the equilibrium value a eq the contact force becomes negative and a decreases further until it reaches its minimum at ͓9͔

͑15͒
where ͑16͒ A yet stronger negative force causes contact failure.The absolute value of F sep can hence be interpreted as the strength of an adhesive bond between the spheres.

Relation between static and dissipative stress tensors
The viscous deformation of particles leads to a dissipative contribution to the stress tensor.For a small deformation rate u ˙ij ͑r ជ , t͒ it reads

͑17͒
where 1 and 2 are the viscous constants, analogous to the elastic constants E 1 and E 2 in Eq. ͑3͒.If the impact velocity g is significantly smaller than the speed of sound of the particles' material and if the characteristic relaxation time of the dissipative processes in the bulk of the material is much shorter than the duration of the impact, the approximation of quasistatic deformation ͓7͔ is eligible.It assumes that the system passes through a sequence of equilibrium states, that is, the slowly varying displacement field coincides at each time instant with the corresponding static field, The superposition of the displacements u ជ H = u ជ H ͑r ជ , H ͒ and u ជ B ͑r ជ , B , a͒, as described in the previous section, leads to a total displacement u ជ st = u ជ H + u ជ B = u ជ͑r ជ , a͒, and in the quasistatic approximation to the displacement rate Then, the dissipative stress reads Comparing the expression for the static stress originating from elastic and adhesive interaction and the previous expression for the dissipative stress, we find that in quasistatic approximation the dissipative stress may be obtained from the corresponding static stress by using viscous constants in place of the elastic ones ͓7͔ and applying the operator a ‫ץ˙‬ / ‫ץ‬a,

͑21͒
Note that although the JKR theory was derived as a static one, it may be applied to slow collisions in terms of the quasistatic approximation addressed here as denoted in Eq. ͑21͒ ͓30͔.

Dissipative force
As it follows from Eq. ͑21͒, the dissipation part of the stress tensor ij dis is expressed through the nondissipative part ij st ; hence, we first need to quantify the latter.In particular, we need the normal component of the static stress tensor at z = 0, which defines the interparticle force.
The dissipative particle interaction is characterized by the normal component of the dissipative stress tensor and thus by the static stress tensor at z = 0, see Eq. ͑21͒.Applying Eqs.͑2͒ and ͑3͒ to both nondissipative contributions of pressure we obtain zz st ͑x,y,z = 0͒ = P H ͑x,y͒ − P B ͑x,y͒ where P H ͑x , y͒ and P B ͑x , y͒ are given by Eqs.͑1͒ and ͑8͒.
The dissipative force acting between two spheres may be computed by integrating the dissipative stress over the contact area at z = 0 using Eqs.͑21͒ and ͑22͒.Instead of this direct computation we apply the method of Refs.͓7,31,32͔, where we transform the coordinate axes as Accordingly, the rescaled contact radius reads a = ␣aЈ, ͑25͒ and yields zz st with substituted coefficients.Thus, using the new coordinates we may write

͑26͒
The expressions in the square brackets on the right-hand sides of Eqs.͑22͒ and ͑26͒ are very similar: While the displacement field u ជ͑r ជ͒ is the same, the coordinates r ជ and r ជ Ј are related by the transformation given in Eq. ͑23͒.Both expressions have the structure of the normal stress at the plane of contact.This similarity may be used to compute the dissipative force.Consider two elastic semispaces ͑z Ͼ 0͒ of the same material.Let the coordinates of the first system r ជ 1 be related to the coordinates of the second system r ជ 2 by an invertible transformation, r ជ 2 = r ជ 2 ͑r ជ 1 ͒.Assume the pressure P 1 ͑x 1 , y 1 ͒ is applied on a domain at z = 0 of the first system and a corresponding pressure P 2 ͑x 2 , y 2 ͒ is applied on a respective domain of the second system.Assume further the displacement field of the first system u ជ 1 ͑r ជ 1 ͒ in the point r ជ 1 coincides with that of the second system at the corresponding point r ជ 2 ͑r ជ 1 ͒, that is, u ជ 1 ͑r ជ 1 ͒ = u ជ 2 (r ជ 2 ͑r ជ 1 ͒).Then the pressures in two systems are related by and the total forces, exerted at z = 0 on the two systems, are equal ͑see Appendix A for the proof͒.Taking into account Eqs.͑19͒, ͑22͒, and ͑25͒ and applying Eq. ͑27͒ to the righthand side of Eq. ͑26͒, separately to each part of the elastic stress tensor with adhesion, we obtain

͑28͒
Then, according to Eq. ͑21͒ the dissipative stress reads By integrating the stress over the contact area we obtain the dissipative force acting between two colliding spheres: where Using Eq. ͑12͒ for F͑a͒ we finally arrive at

͑32͒
For ␥ = 0 this result reduces to the known viscoelastic model ͓7͔.Note the appearance of a cross term in Eq. ͑32͒ which depends on both dissipative and adhesive constants ϳAa ˙ͱ␥a.An earlier model neglected this term thereby overestimating dissipation ͓23,24͔.

A. Equation of motion
The head-on collision of spheres with impact speed g imp is described by the equation of motion where m eff = m 1 m 2 / ͑m 1 + m 2 ͒ is the particles' reduced mass.
Since the total contact force F = F H − F B + F dis is given by Eqs.͑12͒ and ͑32͒ as a function of the contact radius a, it is convenient to write Eq. ͑33͒ in terms of a͑t͒, Here, the prime denotes the derivative with respect to a.For the time evolution of a collision and the according initial a init and final conditions a final we refer to the discussion in the next section.

B. Initial and final collision conditions
The approaching and departing stage of a collision between adhesive particles differ in contact formation and failure ͑for experimental results see, e.g., ͓9,33,34͔͒.While approaching particles collide at = 0 they, if at all, separate at Ͻ 0. This hysteresis is well described in terms of the JKR theory and reflects in the size of the initial a init and final contact radii a final and its relation to the compression ͓see Eq. ͑11͔͒.
While particles approach each other ͑ Ͻ 0͒, the interaction force is zero.Upon contact at = 0, they start touching each other and a strong adhesive force causes a rapid deformation of their surfaces that leads to the formation of an initial contact area of radius a init ͑e.g., ͓34,35͔͒.Thus, the formation of the initial contact area with radius a init occurs very quickly and this strongly nonequilibrium process is accompanied by dissipation of mechanical energy.The energy gained by the reduction of free surface during this process is transformed into the following: ͑a͒ elastic energy of the compressed material, ͑b͒ kinetic energy of the particles which are, thus, accelerated toward each other, and ͑c͒ into heat by viscous material deformation.Once contact is established, adhesive forces continue to accelerate the particles while the Hertzian repulsion and viscous damping decelerate this motion.This marks the beginning of the actual collision with relative impact speed g imp , where the interaction force is described by Eqs.͑12͒ and ͑32͒.
For sufficiently large g imp the particles move toward each other until reaching a maximum compression max .This maximum deformation corresponds to a contact area whose radius is significantly larger than a eq , where the sum of the repulsive Hertzian force and the attracting adhesive force vanishes.While unloading elastic and ͑negative͒ adhesive energy is restored into kinetic energy of relative motion until reaching the critical contact area of radius a final = a sep where the particles separate, see Eq. ͑15͒.The separation corresponds to a negative compression, sep Ͻ 0 ͓30,36͔.This scenario is labeled restitution where the coefficient of restitution falls between 0 Ͻഛ1 ͑see below͒.A general collision process is sketched in Figs.1͑a͒-1͑g͒.
For smaller impact rates a different scenario is observed.While the initial loading stage is similar to the one described above, the unloading stage differs from that of restitution.If the kinetic energy ceases sufficiently during unloading before the particles are completely separated, the compression oscillates until the remaining kinetic energy is dissipated by viscous deformations.Eventually, the oscillation ceases at eq and contact radius a eq where the elastic F H and adhesive forces F B balance.In this case, the particles stick together.Thus, this scenario is labeled aggregation where = 0 and is a consequence of the asymmetry between initial and final contact ͑a final Ͻ a init ͒ in combination with viscous dissipation.
While the Hertzian force, the viscous force, and the adhesive force, which compose the interaction force, are described above, the analysis of the far-from-equilibrium process during the formation of contact is difficult.The rate of material deformation is large here and certainly not well described within the quasistatic approximation.Therefore, for the analysis in this paper we assume that the duration of the initial rapid process of contact formation is negligible compared to the duration of contact, that is, we assume that the initial contact formed instantly.We refer to this process as jumping.
The process of jumping is intimately related with the initial conditions of the collision process.Generally two initial contact radii are possible and are sketched in Fig. 1͑b͒ as types I and II.
For contact formation as in type I, jumping itself does not affect the compression but may be understood as a deformation of the particles' surfaces while the distance of the centers of mass of the particles after jumping is given by the where the initial contact area a init = a 0 corresponds to a vanishing initial compression init = 0 as obtained from Eq. ͑11͒ that reads The initial contact area a 0 is thus smaller than during equilibrium as indicated by experiments ͓33,34,36,37͔.Physically, such an initial condition is justified by the fact that involved time and length scales in jumping are in most cases negligible compared to the macroscopic collision scales.The formation of a relatively small contact area ͑a 0 Ӷ R͒ is caused by strong surface forces.This process, however, occurs very fast ͑about fs͒ and happens in an avalanchelike scenario ͓35͔.Also, if the kinetic energy of the particles at impact is much larger than the surface energy gained in jumping, then this energy does not contribute significantly to enhance the relative particle velocity.Hence neither g imp , nor , can noticeably alter during this very short time interval.For contact formation as in type II, we allow for the extreme case of jumping into an equilibrium contact where adhesive and elastic forces balance at a eq and correspondingly eq as indicated on a microscopic scale ͓35͔.The mathematical formulation of the type II initial condition reads where a eq and eq are given by Eqs.͑13͒ and ͑14͒.Physically, also here we assume an instantaneous process where no additional acceleration is gained although the particle centers shifted.Remarkably, the energy gained from adhesion while 0 ഛ ഛ eq is dissipated by viscous deformation.This follows from the fact that at a = a eq the adhesive force and the elastic Hertz force cancel while the relative particle velocity g imp does not change in this process.
While the first realization ͑type I͒ is supported by experiments ͓33,34,36,37͔, the latter ͑type II͒ is discussed for completeness.A qualitative analysis of the process of contact formation with respect to its duration, particles' acceleration, and effective initial contact radius can be found in Appendix B.
Before closing the section we wish to discuss the range of validity of the quasistatic approximation in the collision dy- 1. Sketch of a collision of adhesive particles.The left column represents the type-I and the right column the type-II initial condition: ͑a͒ The collision starts prior to jumping at time t =−0 when the distance R 1 + R 2 corresponds to the compression =0 at relative velocity g imp .͑b͒ Right after jumping at time t = + 0, a finite contact area of radius a init is formed ͑type I: a init = a 0 ↔ = 0, left; type II: a init = a eq ↔ eq Ͼ 0, right͒.In both cases the impact velocity g imp is preserved, see Eqs. ͑35͒ and ͑37͒.͑c͒ Subsequently, both particles further deform each other until the maximum compression is reached.͑d͒, ͑e͒ Then, they depart and reach = 0 with a nonzero contact radius a. ͑f͒ The adhesive force leads to the formation of a neck and although Ͻ 0 the contact is not broken.͑g͒ Finally, at sep the particles lose contact and the bond is broken.The energy stored in the neck is dissipated by subsequent viscous oscillations of the particles.
namics.As it has been already mentioned, this approximation is applicable if the impact velocity is much smaller than the speed of sound in the particle's material and the collision duration is much shorter than the viscous relaxation time, which may be estimated from the material constants, vis −1 ϳ Y / ͓7͔ ͑here ϳ 1 ϳ 2 ͒.Since Eq. ͑31͒ implies A ϳ / ͑Y 2 ͒, we obtain an estimate for the viscous relaxation time, vis = 2 A. ͑38͒ This is to be compared with the impact time c .While it is difficult to derive c for viscoelastic and adhesive collisions, the impact time of elastic particles may be used as a lower bound ͑e.g., ͓5͔͒, where is the material density of the particles.Moreover, c,el also provides a relevant time scale ͑compression and decompression duration͒ even for sticking collisions, which formally have a divergent collision time.Hence, in what follows we assume the condition vis Ӷ c,el ͑40͒ to be fulfilled, which guarantees the validity of the quasistatic approach.
In the next sections we solve the equation of motion ͑34͒ numerically for given initial, Eqs.͑35͒ and ͑37͒, and final conditions ͑15͒ and study main quantities characterizing the collision, such as the coefficient of restitution and the threshold velocity g cr below which colliding particles aggregate.

IV. COEFFICIENT OF NORMAL RESTITUTION
The coefficient of ͑normal͒ restitution defined by characterizes the loss of kinetic energy of the translational motion due to a collision.In the case of head-on collisions where no energy is transferred from translational into rotational motion, this coefficient characterizes the collision completely, e.g., ͓5,7͔.It may be obtained from the solution of Newton's equation of motion ͑34͒ via where, as before, c denotes the duration of the collision.The coefficient of restitution ranges in the interval ͓0,1͔, where = 0 stands for aggregation and = 1 for elastic restitution.
An interesting application of our theory could be the evolution and dynamics of Saturn's rings.They mainly consist of icy particles whose collision dynamics profoundly affects ring characteristics such as, e.g., their vertical extent ͓6͔.Therefore, for the numerical analysis of Eq. ͑34͒ we chose material parameters of ice at low temperature, Y =7 ϫ 10 9 Pa, = 0.25, =10 3 kg m −3 , and ␥ = 0.74 Nm −1 ͓38͔.The dissipative constant A containing the viscous parameters 1 and 2 is not directly available.Its value A =10 −4 s is the result of a best fit of a viscoelastic noncohesive collision model to experimental data ͓7,39͔.Using these parameters and Eq.͑38͒ we obtain for the viscous relaxation time of ice particles vis Ϸ 6 ϫ 10 −6 s and analyze the collisions which satisfy the condition ͑40͒.
Figure 2 shows the coefficient of restitution as a function of the impact speed g imp for R eff = 1 cm and the above given material parameters for different model realizations.A purely viscoelastic collision, i.e., ␥ = 0, naturally coincides with the underlying viscoelastic model ͓7͔ ͑black dashed line in Fig. 2͒.The behavior of the coefficient of restitution for adhesive particles addressed in this article is qualitatively different from the behavior of nonadhesive viscoelastic particles: While for ␥ = 0 restitution occurs at all impact speeds ͑in particular → 1 for g imp → 0͒, for ␥ 0 particles agglomerate at sufficiently low impact rates ͑ → 0 for g imp → 0͒.In general, the coefficient of restitution increases with decreasing impact rate until reaching a threshold rate g cr where the coefficient of restitution drops to zero.Due to the extra dissipation caused by adhesive forces, adhesive, viscoelastic particle contacts are more dissipative than purely viscoelastic collisions.Thus, even in the absence of observable sticking ͑g imp ഛ g cr practically infeasible͒, the restitution coefficient is overall reduced ͑cf.the dashed black and the solid and dashed-dotted lines in Fig. 2͒.For both types of initial conditions ͑type I, black solid and type II, dashed-dotted lines͒ similar results are obtained.Our results are in good qualitative agreement with experimental studies of icy spheres of various sizes and surface frosts ͓39-43͔ where critical impact rates of maximal 4 mm/ s were reported ͓40-42͔.For exemplary comparison, we plot results of ͓39͔ ͑gray dotted line in Fig. 2͒.Note that here the systematic error in effective mass present in almost all aforementioned experiments due to a utilized disk pendulum as well as a slightly different geometry of an impact into a flat wall has not been accounted for.Further, ͓43͔ concluded that data in ͓39͔ is valid for frost-covered ice spheres of 2.75 cm.Therefore, in order to apply the model presented in this article, careful fits to the existing data are essential.
Figure 3 shows the coefficient of restitution as a function of the effective particle radius for fixed impact rate g imp .Similar to the dependence on the impact rate there are domains of restitution and of aggregation.The critical impact rate g cr below which particles aggregate is larger for smaller particles.Thus, small particles are more likely to stick together while larger ones overall experience more elastic collisions.Moreover, for a given R eff the most elastic collisions occur between same-sized particles ͑black solid or black dotted curve = R 2 / R 1 =1͒, whereas size asymmetry significantly reduces the restitution coefficient ͑black dash-dotted curve, =10͒ shifting g cr by half an order of magnitude.This behavior may be easily understood if one writes the effective mass as The function f͑͒ is maximal for =1 ͑equally sized spheres͒ and decays monotonically with .The larger the inertia, the more elastic the collision and the smaller the relative importance of dissipation and adhesion in collisions.
Hence, for a fixed effective radius R eff the inertia effects, quantified by m eff , would be maximal for = 1.This corresponds to the maximal coefficient of restitution.With increasing size asymmetry the collisions become less restitutive and thus more aggregative.This strongly suggests that the size ratio of the colliding particles, besides a generally low impact rate, is an important parameter favoring a possible aggregation.This particular illustration bears an interesting insight into the collision dynamics in general.In the case of viscoelastic particles where ␥ = 0 the coefficient of restitution increases with the effective radius to the purely elastic limit of =1 ͑black dashed line͒.Contrary, whenever adhesion is involved collisions are always dissipative due to the hysteresis related to adhesive interactions.Collisions of types I and II initial conditions are plotted as black and gray lines, respectively, and clearly illustrate this effect.The larger the difference between a init and a final the more energy is lost during contact.Note that the influence of either type I or II initial condition is rather minor.An adhesive collision is always less elastic than a corresponding viscoelastic one in agreement with experimental studies.
Overall, the differences between viscoelastic ͑␥ =0͒ and adhesive ͑␥ 0͒ viscoelastic collisions are most pronounced for small particles and rather slow impacts as shown in Figs. 2 and 3.

V. STICKING VELOCITY
The agglomeration of slowly colliding particles is a salient property of the viscoadhesive interaction model.The threshold impact velocity g cr distinguishes restitutive ͑g imp Ͼ g cr ͒ from aggregating collisions ͑g imp ഛ g cr ͒ where a joint agglomerate is formed from the constituents.The outcome of a collision sensitively depends on material parameters, impact rate, and size of the colliding particles.From the numerical solutions, as shown in Figs. 2 and 3, we are able to infer the critical velocity for any collision setup ͑see Fig. 4͒.
An analytical estimate of the threshold velocity can be obtained by considering the work done in order to form the neck prior to separation.Since a init a final for ␥ 0, a significant amount of energy will be invested into this hysteresis as discussed above ͑cf.Fig. 3͒.Thus, neglecting energy losses due to the material's viscoelasticity, we estimate the critical impact velocity as ͓44͔ where the adhesive work W ad reads Using Eqs.͑12͒ and ͑11͒ the adhesive work is finally where q = 1.57for a init = a eq and q = 0.09 for a init = a 0 .
Figure 4 shows the critical velocity g cr as a function of the effective particle size R eff .Black solid and dashed-dotted 3. The coefficient of restitution as a function of the effective particle radius R eff for icy particles colliding at a fixed impact rate g imp .Restitution becomes dominant with increasing effective particle radius.Slower impacts and size ratios different from one favor aggregation.The limit of purely elastic collision for large particles is only reached in case of nonadhesive collision.Thus, an adhesive collision is always dissipative.lines denote the results obtained by numerically solving the equations of motion ͑34͒, while gray solid and dashed-dotted lines refer to the corresponding analytical equations ͑44͒ and ͑46͒.Additionally, the qualitative analytical result for the critical velocity obtained by a different approach in a previous study ͓37͔ is plotted for comparison.For large particles the same size dependence is obtained.
Generally, the smaller the particles the higher the critical velocity for restitution to occur.Moreover, g cr is larger for particles of different size ͑ 1͒ than for equally-sized particles with the same R eff .This is easy to understand: Indeed, the relative importance of surface forces that cause aggregation is larger for smaller particles, whereas inertia forces, counteracting aggregation, are smaller as seen before.Aggregation is most likely in these regimes.Numerical results for either type-I or type-II initial conditions are still fairly similar and even identical for small particles.For large R eff , the analytical estimate according to Eqs. ͑44͒ and ͑46͒ reproduces the size dependence of g cr remarkably well.Differences between numerics and analytics are most profound for small particles and are due to the neglect of viscous dissipation in the analytical estimate.Clearly, the analytics underestimates energy losses during collision and thus of the critical velocity g cr .Note that the initial conditions of type II agree much better with the corresponding numerics of g cr in particular in the large particle limit.This is due to the fact that a part of the loading stage between a 0 and a eq , present in type I, is skipped for type-II initial conditions.Skipping this part is equivalent to an assumption that the energy gain, due to the decreasing free surface for the interval a 0 Ͻ a Ͻ a eq , is completely compensated by viscous losses, that is, to an indirect account for dissipation.In the case of small particles, however, the analytics of both type-I and type-II initial conditions differs significantly from the numerics indicating the importance of viscous dissipation.

VI. CONCLUSION
We present an analytical model for the collision dynamics of adhesive viscoelastic spheres which accounts for aggregation and restitution.It is based on the viscoelastic model ͓7͔ where surface interactions based on the JKR theory ͓9͔ have been incorporated.The applicability of this model is based on the quasistatic assumption which implies a relative impact rate much smaller than the material's speed of sound.This allows us to treat the collision process as a sequence of equilibrium states and further justifies the use of the otherwise static JKR and Hertz's theory.We wish to stress, however, that due to the assumed low-velocity collisions no internal shock waves are induced and thus neither plastic deformation nor fragmentation of any kind are covered by this model.
We derive the total, dynamic contact force acting between two colliding spheres.Interestingly, it cannot be expressed as a simple superposition of elastic, adhesive, and dissipative force, since a viscoadhesive cross term arises.Formulating the equation of motion for normal, head-on collisions with appropriate initial conditions we numerically solve the collision dynamics.In particular, we perform simulations of collisions between centimeter-sized icy particles at low temperatures as motivated by an application to the Saturnian main rings.Characteristic parameters, such as the coefficient of restitution and critical impact speed, distinguish between aggregative and restitutive collisions, have been obtained quantitatively, and are in good agreement with experiments.
We observe that aggregation is likely for small and slow particles, whereas restitution occurs otherwise.Further, it is very likely for small particles to stick to larger ones.Adhesive collisions are generally less elastic than corresponding purely viscoelastic impacts but significantly differ only for small R eff , differently sized collision partners 1, or very slow impact rates g imp .In these cases, the relative kinetic energy is completely dissipated leading to an aggregation of particles.This aggregation is possible due to the hysteresis of adhesion and thus the asymmetry of loading and unloading stages.Viscous dissipation further promotes the formation of aggregates in collisions.An analytical estimate for the critical velocity agrees well for large effective particle radii but noticeably deviates for smaller ones.This clearly indicates that viscous dissipation plays an important role in aggregative and restitutive collisions of small viscoelastic particles with adhesion.
Aggregation ͑as well as its counterpart-fragmentation͒ is one of the main driving mechanisms of evolution of multiparticle systems, such as, e.g., planetary rings ͑see a recently developed kinetic model ͓23͔͒.Therefore, an accurate description of particle aggregation is very important.The collision model of adhesive viscoelastic particles ͓23,24͔, adopted in Ref. ͓23͔, yields overall qualitatively similar results ͓23,24͔.It does not incorporate the exact ͑a͒ relation between the contact area and the compression as given in Eq. ͑11͒, but merely traces the collision along the compression and not the contact area.Moreover, it overestimates viscous dissipation by neglecting the viscoadhesive cross term in Eq. ͑32͒.It yields so far higher critical velocities and thus plays a role of an upper bound for the results of the more complete model presented here.FIG. 4. Critical velocity g cr as a function of the effective particle radius for the collision of same-sized icy particles.The analytical estimate given by Eqs.͑44͒ and ͑46͒ is compared to the exact numerical results given by Eq. ͑34͒.The critical impact rate distinguishes restitutive from aggregative collisions, where the domain of restitution is above the respective and the domain of aggregation below each respective curve.The results from Refs.͓23,37͔ show a qualitatively similar particle size dependence of g cr .
In general, the model presented here is capable of qualitatively and quantitatively reproducing the contact dynamics of adhesive viscoelastic particles.For a particular application, this model can easily be fitted to existing experimental results, thus providing an essential tool and addition for the analytical as well as numerical study of multiparticle systems.

APPENDIX A: RELATION BETWEEN PRESSURES AND FORCES IN TWO SYSTEMS
In this section we relate pressure and total force of two systems, whose coordinates and displacement field are related by a linear transformation.Let the pressure P ជ 1 ͑x , y͒ act on the surface of an elastic semispace, z Ͼ 0, of the first system and give rise to a displacement field in the bulk ͓16͔,

͑A1͒
where G ik ͑x , y , z͒ is the corresponding Green's function.
Similarly, one can write for the second system,
Hence the total forces in both systems are equal.

Derivation of the rate equation
Here we present a qualitative analysis of the rapid processes of contact formation.Since both theories, the JKR theory and our theory of the dissipative force, are valid in the quasistatic regime only, we have to apply a more general approach.Therefore we use the energy balance equation supplemented by general expressions for the elastic energy of a deformed sphere and dissipation function of the viscoelastic medium.We start from the first law, where dQ is the heat produced due to the viscous dissipation, dW is the mechanical work, equal to the sum of increments of mechanical potential energy ͑elastic energy in our case͒ and kinetic energy, dW = dU elas + dE kin , and dE is the increment of the internal energy of the particles' material.Hence we recast the above equation into a form The internal energy is related to the free energy F as E = F − T‫ץ‬F / ‫ץ‬T, where T is temperature.For the case of interest only the part of the free energy, related to the varying contact area F =−a 2 ␥, is relevant; this yields that is, we ignore the variation of temperature of the particles.The elastic energy of a sphere, which deformation is obtained by a superposition of Hertzian and Boussinesq's deformations ͑see Sec.II A and II B͒, and characterized by a compression and a contact radius a, reads ͓18͔ where, as previously, H = a 2 / R eff ; note that Eq. ͑B4͒ refers to the general case, where and a are not related by Eq. ͑11͒ as in the JKR theory ͑see, for the derivation, Ref. ͓18͔͒.Strictly speaking, the above expression corresponds to a static case, when the external parameters such as and a completely determine stress and/or strain distribution in a bulk.Nevertheless, it is still a good approximation if the displacement rate is significantly smaller than the speed of sound in the material-this enables the stress and/or strain to relax to their static values.It is also worth it to note that since the static JKR theory is essentially based on two relations Eqs.͑B3͒ and ͑B4͒, for E͑a͒ and U elas ͓27͔, its validity extends to dynamic applications, as long as the above condition is met.We assume that this is the case and prove the self-consistency of the assumption.
The kinetic energy may be written as where u ជ ˙͑r ជ͒ is the displacement rate at a point r ជ, and is the material density, unaltered by small deformations; the integration is performed over the particles' volumes.Finally, the dissipation rate Q ˙may generally be written as where u ˙ij is the strain rate and R is the dissipation function of the viscoelastic material ͓16͔; the integration is performed over the same domain.The negative sign in front of Q ˙follows from the fact that the heat is produced by the system, and not transmitted to it from outside.We start the analysis for the initial conditions of the first kind, Eq. ͑35͒, and consider for simplicity a collision of a sphere with a plane, that is, R eff = R.In this case = 0 and Eq.͑B4͒ yields U elas = ͑3/20͒a 5 / DR 2 .The total configurational energy in this case U tot = U elas + E = ͑3/20͒a 5 /DR 2 − a 2 ␥ ͑B7͒ is always negative for small a and has a ͑negative͒ minimum, dU tot / da =0, at a = a 0 .Once brought into a contact, the system "slides" down the energy gradient, enlarging the contact area, and tends to reach the potential minimum at a = a 0 .To obtain the rate of this process one needs to determine E kin and dQ / dt and use the kinetic equation ͑B2͒.However, it is not easy to find these quantities, hence we apply the qualitative estimates where ϳ 1 ϳ 2 is the characteristic viscous constant, u ānd u ¯ij are respectively characteristic displacement and strain, and ⌬V is the characteristic volume of the domain where the deformation mainly occurs.It follows from a simple geometric analysis ͑see Fig. 5͒ that For the characteristic displacement u ¯we take this quantity for the point on the sphere surface on the boundary of the contact zone, where u = a 2 / R ͑see Fig. 5͒.To estimate the strain u ¯ij ϳ ⌬u / ⌬x we take into account that the displacement of the material on the sphere surface varies from zero, u =0, at the center of the contact zone to u = a 2 / R on the contact boundary of radius a, that is, ⌬u = a 2 / R and ⌬x = a.The volume ⌬V is estimated as a volume difference between the truncated cone with the two bases of 2a and a and height h =2a 2 / R and a spherical cap of the same height and basis radius 2a, see Fig. 5.
Using Eq. ͑B10͒ we recast Eqs.͑B8͒ and ͑B9͒ into the form where c is the characteristic speed of sound in the material.
In new variables the energy balance reads
FIG. 5. Illustrates the geometry of the contact formation for the initial conditions ͑35͒.For simplicity the collision of a sphere with a plane is considered.The radius of the sphere is R and the current radius of the contact area is a.The displacement of the sphere's surface is zero at the center of the contact zone, u͑0͒ = 0, while at the boundary of the contact zone, that is, on the distance a from the center u͑a͒ = a 2 / R.This yields the estimate for the characteristic strain, u ¯ij ϳ͑a 2 / R͒ / a ϳ a / R. The shadowed region corresponds to the domain of the most intensive deformation.Its volume is equal to the volume difference of the truncated cone with basis radii 2a and a and height h = a sin ␤, and the spherical cap of the same height h and basis radius 2a.For a / R Ӷ 1, from simple geometry, follows tan͑␤ /2͒Ӎ␤ /2=a / R or tan ␤ Ӎ sin ␤ Ӎ 2a / R. The velocity of a material point on the sphere's surface u ˙z͑a͒ is related to the rate of the contact enlargement a ˙͑which is the speed of a geometric point͒ as u ˙z͑a͒ = a ˙tan ␤ Ӎ a ˙͑2a / R͒.
We consider separately two limiting cases, v Ӷ 0 ͑negligible dissipation͒ and the opposite case, v ӷ 0 .

Contact formation for negligible dissipation
In the first case Eq. ͑B14͒ may be reduced to where we take into account that initially a ˆ=a ˆ= 0. We assume that the contact is practically formed when a = ␣a 0 ͑or a ˆ= ␣͒ where ␣ ഛ 1; here we choose ␣ = 0.99.Solving Eq. ͑B15͒ and returning to the dimensional time we obtain the duration of the contact formation, ˆ3 0 .͑B16͒ The integral may be easily performed with the result ͑5/3͒ ϫ͑1− ͱ 1−2␣ 3 /5͒, which yields, using the above value of ␣ = 0.99, Eq. ͑36͒ for a 0 , and the definition of D ϳ 1/Y, For other values of ␣, different numerical coefficients are obtained.Naturally, we skip these for our qualitative estimate to 1 .Note that according to Eq. ͑B17͒ the formation of a contact of radius a 0 occurs ͱ R / a 0 times faster than ͑a 0 / c͒-the time that the sound travels the distance a 0 .That is, the boundary of the contact zone effectively moves with the speed ͱ R / a 0 times larger than the speed of sound.This, however, is not a speed of a material particle, but rather of a geometric point.Indeed, as it is illustrated in Fig. 5, the velocity at which the surface of the sphere moves toward the plane, u ˙z͑a͒, taken at a contact radius a is related to the speed of the contact boundary a ˙as u ˙z͑a͒ = a ˙tan ␤, where tan ␤ =2a / R, see Fig. 5. Since, as noticed previously, a ˙= c ͱ R / a, we obtain u ˙z͑a͒ = ͑2a / R͒c ͱ R / a =2c ͱ a / R or u ˙z͑a͒ Ӷ c.That is, in spite of the supersonic spread of the contact boundary, the material itself moves with a speed significantly smaller than the speed of sound.This proves the self-consistency of our approach.

Contact formation with dissipation
Next we consider the case of large viscosity, v ӷ 0 , which seems to be of a wider application ͑this condition is fulfilled for the ice particles addressed here͒.In this case the first term on the left-hand side of Eq. ͑B14͒ is ͑ 0 / v ͒ smaller than the term on the right-hand-side, and hence may be neglected.Without this term this equation may be recasted into the form a ˆ= 3 4 ͑1 − a ˆ3͒a ˆ−3 ͩ 0 R v a 0 ͪ.

͑B18͒
Solving this equation in the same way as Eq.͑B15͒, with the same condition for the contact formation, a ˆ= ␣ = 0.99, and transforming back to the dimensional time, we obtain the duration of this process, . ͑B19͒ As expected, contact formation takes longer for more viscous material.Nevertheless, it still occurs fast enough to assume that during this process the impact velocity persists and compression keeps constant, =0.

Analysis of validity of initial conditions
To prove the validity of the initial conditions we consider the attractive force which pulls the particles toward each other at the beginning of the contact at = 0.This may be done differentiating the total configuration energy U tot = U elas + E with respect to the "coordinate" , where U elas is generally given by Eq. ͑B4͒ and we take into account that E does not depend on .The force increases with increasing contact area, so that the increment of impulse m⌬g imp due to this force, acting during the time interval 2 , reads a ˆ͑a ˆ͒ da ˆ,

͑B21͒
where we transform to the dimensionless variables and change the integration variable from t ˆto a ˆ.Substituting Eq. ͑B18͒ for a ˆ͑a ˆ͒ into the last equation and performing the integration we obtain for ␣ = 0.99, . ͑B22͒ The validity of the initial conditions requires the prerequisite ⌬g imp Ӷ g imp .For the centimeter-sized particles we find ⌬g imp ϳ 10 −5 m / s, where the relation ϳ AY, which follows from Eq. ͑31͒, has been used.Next, we estimate the increment of the compression ⌬ for the time of the contact formation 2 .Using m ¨= F with the force F Ӎ F͑ =0,a͒ as it is given by Eq. ͑B20͒, we write Performing calculations similar to those leading to Eq. ͑B22͒, we finally arrive at ⌬ = 3.02 It is worth it to compare ⌬ with the value of eq , which characterizes the starting compression for the initial conditions of the second type.Using Eqs.͑14͒ and ͑B24͒ we find where in the last equations we again exploit ϳ AY.For the icy particles addressed here this ratio is estimated as 10 −5 , which justifies the initial conditions ͑35͒ used in our study.
sum of their radii.A mathematical formulation of the type I initial t→+0 a͑t͒ = a init = a 0 Ͻ a eq , ͑35͒ FIG.3.The coefficient of restitution as a function of the effective particle radius R eff for icy particles colliding at a fixed impact rate g imp .Restitution becomes dominant with increasing effective particle radius.Slower impacts and size ratios different from one favor aggregation.The limit of purely elastic collision for large particles is only reached in case of nonadhesive collision.Thus, an adhesive collision is always dissipative.

Ε Bridges et al. 39 Brilliantov et al. 7 Albers and Spahn 24 Eq. 34,42 : a init a eq Eq. 34,42 : a init a 0
FIG.2.Coefficient of restitution as a function of the impact rate g imp for same-sized icy particles of radius R =2 cm ͑R eff =1 cm͒ for different model realizations.As expected and in general agreement with experiments ͓40,41͔, adhesive collisions are more dissipative than purely viscoelastic ones ͑cf. the dashed black and the solid and dashed-dotted lines͒.Moreover, below a certain impact speed, here g cr Ϸ 2 mm/ s, particles stick together instead of rebouncing in the case of restitution for g imp Ͼ g cr .Results of ͓39͔ are shown as an exemplary comparison.