# Granular gas of viscoelastic particles in a homogeneous cooling state.

###### Abstract

Kinetic properties of a granular gas of viscoelastic particles in a homogeneous cooling state are studied analytically and numerically. We employ the most recent expression for the velocity-dependent restitution coefficient for colliding viscoelastic particles, which allows to describe systems with large inelasticity. In contrast to previous studies, the third coefficient of the Sonine polynomials expansion of the velocity distribution function is taken into account. We observe a complicated evolution of this coefficient. Moreover, we find that is always of the same order of magnitude as the leading second Sonine coefficient ; this contradicts the existing hypothesis that the subsequent Sonine coefficients , , are of an ascending order of a small parameter, characterizing particles inelasticity. We analyze evolution of the high-energy tail of the velocity distribution function. In particular, we study the time dependence of the tail amplitude and of the threshold velocity, which demarcates the main part of the velocity distribution and the high-energy part. We also study evolution of the self-diffusion coefficient and explore the impact of the third Sonine coefficient on the self-diffusion. Our analytical predictions for the third Sonine coefficient, threshold velocity and the self-diffusion coefficient are in a good agreement with the numerical finding.

###### keywords:

Granular gas; Viscoelastic particles; Velocity dependent coefficient of restitution; Velocity distribution function; Self-diffusion^{†}

^{†}journal: Physica A

## 1 Introduction

An ensemble of macroscopic particles, which move ballistically between dissipative collisions, is usually termed as a granular gas, e.g. book ; Goldhirsh_AnRev2003 ; hinrichsen ; book1 ; book2 and the loss of energy at impacts is quantified by a restitution coefficient ,

(1) |

Here and are relative velocities of two particles before and after an impact and is the unit vector, connecting their centers at the collision instant. In what follows we assume that the particles are smooth spheres of diameter , that is, we do not consider the tangential motion of particles. Then the restitution coefficient yields the after-collision velocities of particles in terms of the pre-collision ones,

(2) |

The simplest model for the restitution coefficient, , facilitates significantly the theoretical analysis and often leads to qualitatively valid results, e.g. book ; Goldhirsh_AnRev2003 . This assumption, however, contradicts the experimental observations, e.g. goldsmit ; bridges ; kuwabara and basic mechanical laws, e.g. tanaka ; ramirez . A first-principle analysis of a dissipative collision may be performed, leading to a conclusion that must depend on the relative velocity of the colliding particles kuwabara ; bril96 ; Morgado97 . If the impact velocity is not very high (to avoid plastic deformation of particles material) and not too small (to avoid adhesive interactions at a collision bril2007 ), the viscoelastic contact model may be applied bril96 . This yields the velocity-dependent restitution coefficient schwager1998 ; ramirez :

(3) |

Here , and is the dimensionless temperature, expressed in terms of current granular temperature,

(4) |

with being the velocity distribution function of grains, is a gas number density, is the particle mass and is the thermal velocity; is, correspondingly, the dimensionless impact velocity.

The restitution coefficient (3) depends on the small dissipation parameter

(5) |

proportional to the dissipative constant ,

(6) |

which depends on the Young modulus of the particle’s material , its Poisson ratio and the viscous constants and , see e.g. book . The parameter reads:

(7) |

Due to inelastic nature of inter-particle collisions, the velocity distribution function of grains deviates from the Maxwellian, so that the dimensionless distribution function esipov

(8) |

may be represented in terms of the Sonine polynomials expansion, e.g. goldshtein ; vannoije ; book ; pf :

(9) |

where is the dimensionless Maxwellian distribution and the first few Sonine polynomials read book ,

(10) | |||

(11) | |||

(12) |

Provided that the expansion (9) converges, the Sonine coefficients completely determine the form of . According to the definition of temperature book , that is, the first non-trivial coefficient in the expansion (9) is . For the case of a constant the coefficient has been found analytically goldshtein ; vannoije and the coefficient analytically breakdown and numerically breakdown ; Brey1996 ; Santos1999 ; Nakanishi2003 . The expansion (9) quantifies the deviation of the distribution function from the Maxwellian for the main part of the velocity distribution, that is, for ; the high-velocity tail of for is exponentially over-populated esipov ; vannoije ; book ; pf and requires a separate analysis.

The impact-velocity dependence of the restitution coefficient, as it follows from the realistic visco-elastic collision model, has a drastic impact on the granular gas properties. Namely, the form of the velocity distribution function and its time dependence significantly change vdve , similarly changes the time dependence of the kinetic coefficients Bril_PRE2003 ; Bril_in_NLP . Moreover, the global behavior of the system qualitatively alters: Instead of evolving to a highly non-uniform final state of a rarified gas and dense clusters, as predicted for mcnamara ; goldhirsch , the clustering goldhirsch and the vortex formation BritoErnst98 occurs in a gas of viscoelastic particles only as a transient phenomena BrilPRL2004 .

Based on the restitution coefficient given by Eq. (3) the theory of granular gases of viscoelastic particles has been developed, e.g. book , where as in the case of a constant , only the first non-trivial Sonine coefficient has been taken into account. Although the evolution of the high-velocity tail of has been analyzed vdve , neither the location of the tail, nor its amplitude was quantified.

Recently, however, it has been shown for the case of a constant restitution coefficient, that the next Sonine coefficient is of the same order of magnitude as the main coefficient breakdown ; moreover, it was also shown that the amplitude of the high-velocity tail of and its contribution to the kinetic coefficients may be described quantitatively vdeps ; tail2007 .

Finally, a new expression for the velocity-dependent restitution coefficient has been derived delayed , which takes into account the effect of ”delayed recovery” in a collision. The delayed recovery implies, that at the very end of an impact, when the colliding particles have already lost their contact, their material remains deformed delayed . This affects the total dissipation at a collision, so that the revised reads delayed :

(13) |

Here , , , and the other numerical coefficients up to are given in delayed . Fig. 1 (inset) illustrates the dependence of the restitution coefficient (13) on the dissipative parameter for for collisions with the characteristic thermal velocity, that is for .

In the present study we use the revised expression (13) for and investigate the evolution of the velocity distribution function in a gas of viscoelastic particles in a homogeneous cooling state. With the new restitution coefficient we are able to describe collisions with significantly larger dissipation, than it was possible before with the previously available expression for . For the larger inelasticity one expects the increasing importance of the next-order terms in the Sonine polynomials expansion and of the high-energy tail of the velocity distribution. In what follows we study analytically and numerically time evolution of , using the Sonine expansion up to the third-order term and analyze the amplitude and slope of the high-energy tail of . In addition, we consider self-diffusion – the only non-trivial transport process in the homogeneous cooling state and compute the respective coefficient.

The rest of the paper is organized as follows. In the next Sec. II we address the Sonine polynomial expansion and calculate the time-dependent Sonine coefficients, along with the granular temperature. In Sec. III the high-energy tail is analyzed, while Sec. IV is devoted to the self-diffusion coefficient. Finally, in Sec. V we summarize our findings.

## 2 Sonine polynomial expansion: evolution of the expansion coefficients

Evolution of the velocity distribution function of a granular gas of spherical particles of diameter in a homogeneous cooling state obeys the Boltzmann-Enskog equation, e.g. book ; Goldhirsh_AnRev2003 :

(14) |

where is the collision integral. Generally, depends on the two-particle distribution function . Within the hypothesis of molecular chaos, (see e.g. book ), the closed-form equation (14) for is obtained; denotes here a pair correlation function at a contact, which accounts for the increasing collision frequency due to the effect of excluded volume book .

Eq. (14) yields the equation for the dimensionless distribution function book ,

(15) |

where and

(16) |

is the -th moment of the dimensionless collision integral:

(17) |

in the above equation is the Heaviside step-function and the dimensionless velocities and are the pre-collision velocities in the so-called inverse collision, which results with and as the after-collision velocities, e.g. book .

Eq. (15) is coupled to the equation for the granular temperature, e.g. book :

(18) |

which defines the cooling coefficient .

Multiplying Eq. (15) with , integrating over and writing Eq. (18) in the dimensionless form, we obtain for book (similar equations has been used in huthmann ):

(19) |

Here is the dimensionless time, measured in the mean collision units, at initial time . The equations for the Sonine coefficients and in (19) are the first two of the infinite system of equations for huthmann ; book . The moments in these equations depend on all the Sonine coefficients , with , that is, only the infinite set of equation is closed. To make the system tractable one has to truncate it. Following Refs. breakdown ; vestnik we truncate the Sonine series at the third term, i.e. we approximate,

(20) |

Within this approximation in Eqs. (19) depend only on , and , which makes the system (19) closed.

The moments of the collisional integral were calculated analytically up to , using the formula manipulation program, explained in detail in book . The complete expressions are rather cumbersome, therefore, we present here for illustration only linear approximations with respect to , and :

(21) | |||

where .

Since the system of equations (19) is strongly non-linear, only a numerical solution is possible. Still, one can find a perturbative solution in terms of small dissipative parameter ,

(22) |

The zeroth-order solution (with ) reads,

(23) |

where . For the first-order solution only its asymptotics for may be obtained analytically:

(24) | |||

(25) | |||

(26) | |||

Interestingly, the third Sonine coefficient is of the same order of magnitude with respect to the small parameter , as the second Sonine coefficient , albeit five times smaller. This conclusion is in a sharp contrast with the hypothesis suggested in Ref. huthmann , that the Sonine coefficients are of an ascending order of some small parameter , that is, .

The numerical solution of Eqs. (19) confirms the obtained asymptotic dependence, Eqs. (24)-(26). This is seen in Fig. 1, where the time dependence of the reduced temperature is plotted and in Figs. 2 and 3 (insets). Fig. 1 demonstrates that the larger the dissipative parameter , the earlier the asymptotic behavior of is achieved.

The numerical solution for the Sonine coefficients and , shown in Fig. 2 and 3, corresponds to the initial Maxwellian distribution, . As it follows from the figures, the absolute values of the both coefficients initially increase, reach their maxima and and, eventually, decay to zero. In other words, the velocity distribution function for viscoelastic particles evolves towards the Maxwellian. It is interesting to note, that the maxima and first increase with the increasing dissipative parameter , then saturate at and do not anymore grow. The location of the maxima, however, shifts with increasing to later time, Figs. 2 and 3. This illustrates the general tendency – the larger the dissipation parameter, the slower the gas evolution. Again we see that the third Sonine coefficient is of the same order of magnitude as (although a few times smaller) for all evolution stages.

To understand the observed behavior of the Sonine coefficients, consider the dependence of these coefficients on for a constant restitution coefficient (Fig. 4) and the respective dependence of the restitution coefficient on the dissipative parameter (Fig. 1, inset). For [which corresponds to , Fig. 1 (inset)] one observes a fast relaxation, on a collision time scale, of the Sonine coefficients , to their maximal values , , roughly corresponding to the respective values for the constant . In a course of time the granular temperature decreases and the effective restitution coefficient alters in accordance with decreasing dissipative parameter , that is, the effective increases with time. For this interval of () the increasing restitution coefficient implies the decrease of the absolute values of the Sonine coefficients, Fig. 4, which is indeed observed in the evolution of and . The larger the (i.e. the smaller the effective restitution coefficient), the larger the maxima and , in accordance with the dependencies and depicted in Fig. 4.

Similarly, for [which corresponds to , Fig. 1 (inset)] the initial fast relaxation to the related values of and first takes place. Then the decreases of and the respective increase of the effective causes the increase of the absolute value of , until it reaches the maximum, corresponding to (or ), Fig. 4. The further decrease of leads to the according decay of , in agreement with the dependence of shown in Fig. 2. This qualitatively explains the evolution of and the saturation of its maximum for . The qualitative behavior of the third Sonine coefficient may be explained analogously.

It is interesting to note that the observed dependence of with the new restitution coefficient (13) differs qualitatively for from that obtained previously for the old restitution coefficient (3). While in the latter case a positive bump at initial time, was detected for vdve , in the former case the positive bump appears at much larger , Fig. 5. This is again in agreement with the behavior, expected from the dependence of for a constant : The coefficient becomes positive for , which corresponds to , Fig. 1 (inset). Note, however, that with increasing more and more terms in the expansion (13) are to be kept. While for it suffice to keep terms up to , for the Sonine coefficients computed with the accuracy , and noticeably differ, Fig. 5. Therefore we conclude that the revised restitution coefficient for a viscoelastic impact (13) may be accurately used up to . The loss of the accuracy in the computation of the Sonine coefficients for larger values of may be a manifestation of the breakdown of the Sonine polynomials expansion, as it has been found previously for a constant restitution coefficient breakdown . Similarly as in the case of constant breakdown , we expect that in the domain of convergence of the Sonine expansion, , the magnitude of the next-order Sonine coefficients , is very small, so that an acceptable accuracy may be achieved with the use of the two coefficients, and only.

## 3 High-velocity tail of the velocity distribution function

The expansion (9) refers to the main part of the velocity distribution, , that is, to the velocities, close to the thermal one, . The high-velocity tail is however exponentially overpopulated esipov . It develops in a course of time, during a first few tens of collisions vdeps ; tail2007 . For viscoelastic particles the velocity distribution reads for vdve :

(27) |

where the function satisfies the equation vdve ,

(28) |

with and defined previously. Using , obtained by the formula manipulation program book (see the discussion above), we solve numerically Eq. (28) to obtain . In the linear approximation, , the function has the form:

(29) |

with book .

Following vdeps we neglect the transition region between the main part of the velocity distribution function, , and its high-energy part, and write the distribution function as

(30) | |||||

The coefficients , and the threshold velocity , which separates the main and the tail part of can be obtained, using the normalization condition:

(31) |

and the continuity condition for the function itself and its first derivative vdeps :

(32) | |||

(33) |

Substituting Eq. (30) into (31), (32) and (33) we arrive at,

(34) |

where is the solution of Eq. (28).

To find the amplitudes and together with the threshold velocity the system (34) was solved numerically.

The asymptotic dependence of on may be easily found if we take into account that and are of the same order of magnitude for , while and . Keeping in the first equation in (34) only the largest terms, yields , which implies that

(35) |

The typical velocity distribution function, computed for and is shown in Fig. 6. The threshold velocity, separating the main and the tail part of the velocity distribution reads in this case, . The threshold velocity increases with decreasing inelasticity and shifts to larger values at later time, Fig. 7, which means that the high-energy tail becomes less pronounced. Again, we see that in contrast to a gas of particles with a constant restitution coefficient, where the tail of persists after its relaxation, in a gas of viscoelastic particles the velocity distribution function tends to a Maxwellian.

## 4 Self-diffusion

Self-diffusion is the only transport process which takes place in a granular gas in a homogeneous cooling state: In spite of the lack of macroscopic currents, a current of tagged particles, identical to the particles of the surrounding gas, but somehow marked, may exist. Moreover, the diffusion coefficient is directly related to the mobility coefficient of the tagged particles via the Einstein relation, , which approximately holds true for granular gases, e.g. Puglisi2004 ; Garzo2004 . The mean-square displacement of the tagged particles reads book

(36) |

where the time-dependent diffusion coefficient (diffusivity) is the solution of the following equation book :

(37) |

Here is, as previously, the cooling coefficient and is the velocity correlation time, which characterizes the time after which the memory about the initial particle velocity is lost; it reads in terms of the distribution function book :

(38) | |||

With the approximation (20) for was calculated up to , using the formula manipulation program, described in book . Here we present for illustration only its linear with respect to , and part:

(39) |

Here is the Enskog velocity correlation time of elastic particles at initial time :

(40) |

Using the obtained expressions for and up to , we solve Eq. (37) numerically and compute the diffusion coefficient . In Fig. 8 the ratio of is plotted, where

(41) |

is the Enskog self-diffusion coefficient for a gas of elastic particles at initial temperature . As it is clearly seen from the figure, the diffusion coefficient decreases with time in a way, that the ratio approaches – the ratio of the two diffusion coefficients for elastic particles at the current temperature and the initial temperature . Hence in a course of time the self-diffusion in a gas of viscoelastic particles tends to that in a gas of elastic particles.

In the linear approximation with respect to the small dissipative parameter one can obtain an analytical expression for . Keeping only first-order terms in the expressions for , , and and substituting them into Eq. (37) yields the diffusion coefficient:

(42) |

where, as previously, and . Note that due to the account of the third Sonine coefficient in and the linear approximation (42) for differs from the previously obtained result book .

The linear approximation (42) is compared in Fig. 8 with the complete numerical solution. As it follows from the figure, after about ten collision per particle the two solution become practically indistinguishable. As it may be seen from the inset of Fig. 8, the impact of the third Sonine coefficient on the behavior of self-diffusion coefficient is rather small even for the large value of dissipation parameter .

## 5 Conclusion

We study evolution of a granular gas of viscoelastic particles in a homogeneous cooling state. We use a new expression for the restitution coefficient , which accounts for the delayed recovery of the particle material at a collision and allows to model collisions with much larger dissipation, as compared to previously available result for . We analyze the velocity distribution function and the self-diffusion coefficient. To describe the deviation of the velocity distribution function from the Maxwellian we use the Sonine polynomial expansion. In contrast to the commonly used approximation, which neglects all terms in the Sonine expansion beyond the second one, we consider explicitly the third Sonine coefficient. We detect a complicated evolution of this coefficient and observe that it is of the same order of magnitude, with respect to the (small) dissipative parameter , as the second coefficient. This contradicts the existing hypothesis huthmann , that the subsequent Sonine coefficients , , are of an ascending order of some small parameter, characterizing particles inelasticity. Similarly as for the case of a constant restitution coefficient, we obtain an indication of divergence of the Sonine expansion for large dissipation, . For the asymptotic long-time behavior we derive analytical expressions for both Sonine coefficients, which agree well with the numerical data.

Using the obtained third Sonine coefficient we compute the self-diffusion coefficient and derive an analytical expression for in the linear, with respect to , approximation. We show that the complete solution approaches the approximate analytical solution after a transient time of about ten collisions per particle. We observe, that in spite of the importance of for an accurate description of the velocity distribution function, its impact on is rather small.

We also study the evolution of the high-energy tail of the velocity distribution. Using the equation for the time-dependent slope of the tail and the obtained Sonine coefficients we find the amplitude of the tail and the threshold velocity, which demarcates the main part of the velocity distribution and the high-energy tail. We find the analytical expression for the asymptotic behavior of the threshold velocity, which agrees well with numerics, and observe, that in a course of time it shifts to larger values; this implies that the high-velocity tail becomes less pronounced.

Such behavior of the threshold velocity, of the Sonine coefficients and of the coefficient of self-diffusion, naturally, manifests that the properties of the system tend to those of a gas of elastic particles; our theory quantifies evolution towards this limit.

## References

- (1) Brilliantov N V and Pöschel T 2004 Kinetic theory of Granular Gases (Oxford University Press)
- (2) Goldhirsch I 2003 Annual Review of Fluid Mechanics 35 267
- (3) Hinrichsen H, Wolf D E 2004 The Physics of Granular Media Berlin: Wiley
- (4) Pöschel T and Luding S 2001 Granular Gases 564 of Lecture Notes in Physics Berlin: Springer
- (5) Pöschel T and Brilliantov N V 2003 Granular Gas Dynamics 624 of Lecture Notes in Physics Berlin: Springer
- (6) Goldsmit W 1960 The Theory and Physical Behavior of Colliding Solids (London: Arnold)
- (7) Bridges F G, Hatzes A and Lin D N C 1984 Nature 309 333
- (8) Kuwabara G, Kono K 1987 J. Appl. Phys. Part 1 26 1230
- (9) Tanaka T, Ishida T and Tsuji Y 1991 it Trans. Jap. Soc. Mech. Eng. 57 456
- (10) Ramirez R, Pöschel T, Brilliantov N V and Schwager T 1999 Phys. Rev. E 60 4465
- (11) Brilliantov N V, Spahn F, Hertzsch J-M, Poeschel T 1996 Phys. Rev. E. 53 5382
- (12) Morgado W A M and Oppenheim I 1997 Phys. Rev. E 55 1940
- (13) Brilliantov N V, Albers N, Spahn F and Poeschel T 2007 Phys. Rev. E 76 051302
- (14) Schwager T and Pöschel T 1998 Phys. Rev. E 57 650
- (15) Esipov S E and Pöschel T 1997 J. Stat. Phys. 86 1385
- (16) Goldshtein A and Shapiro M 1995 J. Fluid. Mech. 282 75
- (17) van Noije T P C and Ernst M H 1998 Granular Matter 1 57
- (18) Brilliantov N V and Pöschel T 2006 Europhys. Lett. 74 424
- (19) Brey J J, Cubero D and Ruiz-Montero M J 1996 Phys. Rev. E 54 3664
- (20) Montanero J M and Santos A 1999 Granular Matter 2 53
- (21) Nakanishi H 2003 Phys. Rev. E 67 010301R
- (22) Noskowicz S H, Bar-Lev O, Serero D, Goldhirsch I 2007 Europhys. Lett. 79 60001
- (23) Brilliantov N V and Pöschel T 2000 Phys. Rev. E 61 5573
- (24) Brilliantov N V and Pöschel T 2003 Phys. Rev. E 67 061304
- (25) Pöschel T and Brilliantov N V 2003 in book2 p.131
- (26) McNamara S 1993 Physics of Fluids A 5 3056
- (27) Goldhirsch I and Zanetti G 1993 Phys. Rev. Lett. 70 1619
- (28) Brito R and Ernst M H 1998 Europhys. Lett. 43 497
- (29) Brilliantov N V, Saluena C, Schwager T and Poeschel T 2004 Phys. Rev. Lett. 93 134301
- (30) Pöschel T, Brilliantov N V and Formella A 2006 Phys. Rev. E 74 041302
- (31) Pöschel T, Brilliantov N V and Formella A 2007 Int. J. Mod. Phys. C 18 701
- (32) Schwager T and Pöschel T 2008 Phys. Rev. E 78 051304
- (33) Bodrova A S and Brilliantov N V 2009 Moscow University Physics Bulletin 64 131
- (34) Huthmann M, Orza J A and Brito R 2000 Granular Matter 2 189
- (35) Barrat A , Loreto V and Puglisi A 2004 Physica A 334 513
- (36) Garzo V 2004 Physica A 343 105