Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://poly.phys.msu.ru/~bodrova/PhysA.pdf Äàòà èçìåíåíèÿ: Mon Dec 14 16:28:32 2009 Äàòà èíäåêñèðîâàíèÿ: Sun Apr 10 22:42:38 2016 Êîäèðîâêà:

This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier's archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright


Author's personal copy
Physica A 388 (2009) 3315­3324

Contents lists available at ScienceDirect

Physica A
journal homepage: www.elsevier.com/locate/physa

Granular gas of viscoelastic particles in a homogeneous cooling state
A.S. Bodrova a, , N.V. Brilliantov
a b

b,a

Faculty of Physics, M.V.Lomonosov Moscow State University, Leninskie Gory, Moscow, 119991, Russia Department of Mathematics, University of Leicester University Road, Leicester LE1 7RH, UK

article

info

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 us to describe systems with large inelasticity. In contrast to previous studies, the third coefficient a3 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 a3 is always of the same order of magnitude as the leading second Sonine coefficient a2 ; this contradicts the existing hypothesis that the subsequent Sonine coefficients a2 , a3 . . ., 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 D 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. © 2009 Elsevier B.V. All rights reserved.

Article history: Received 10 December 2008 Received in revised form 17 March 2009 Available online 8 May 2009 Keywords: Granular gas Viscoelastic particles Velocity dependent coefficient of restitution Velocity distribution function Self-diffusion

1. Introduction An ensemble of macroscopic particles, which move ballistically between dissipative collisions, is usually termed as a granular gas, e.g. Refs. [1­5] and the loss of energy at impacts is quantified by a restitution coefficient ,

=

v12 · e v12 · e

.

(1)

Here v12 and v12 are relative velocities of two particles before and after an impact and e 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) 2 2 The simplest model for the restitution coefficient, = const, facilitates significantly the theoretical analysis and often leads to qualitatively valid results, e.g. Refs. [1,2]. This assumption, however, contradicts the experimental observations,

v1 = v1 - (1 + ) v12 · e e,

1

v2 = v2 + (1 + ) v12 · e e.

1



Corresponding author. Tel.: +7 495 452 14 54. E-mail address: bodrova@polly.phys.msu.ru (A.S. Bodrova).

0378-4371/$ ­ see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2009.04.040


Author's personal copy
3316 A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324

e.g. Refs. [6­8] and basic mechanical laws, e.g. Refs. [9,10]. 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 [8,11,12]. 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 [13]), the viscoelastic contact model may be applied [11]. This yields the velocity-dependent restitution coefficient [14,10]:

= 1 - C1 (2u(t ))1
Here C1 1.15, C temperature, 3 2 nT (t ) =
2

/10

c12 e

1/5

+ C2 2 (2u(t ))1/

5

c12 e

2/5

± ....

(3)

0.798 and u(t ) = T (t )/T (0) is the dimensionless temperature, expressed in terms of current granular mv 2
2

dv

f (v , t ) =

mvT (t ) 2

,

(4)

with f (v , t ) being the velocity distribution function of grains, n is a gas number density, m is the particle mass and vT (t ) is the thermal velocity; c12 = v12 /vT is, correspondingly, the dimensionless impact velocity. The restitution coefficient (3) depends on the small dissipation parameter

= A

2 5

T (0) m

1 10

,

(5)

proportional to the dissipative constant A, A= 1 (32 - 1 )2 (1 - 2 )(1 - 2 ) 3 (32 + 21 ) Y
2

,

(6)

which depends on the Young modulus of the particle's material Y , its Poisson ratio and the viscous constants 1 and 2 , see e.g. Ref. [1]. The parameter reads:

=

3 2

3 2

. m 1 - 2
Y

(7)

Due to the inelastic nature of inter-particle collisions, the velocity distribution function of grains deviates from the Maxwellian, so that the dimensionless distribution function [15] f (v , t ) = n
3 vT

~ f (c , t ),

(8)

may be represented in terms of the Sonine polynomials expansion, e.g. Refs. [16,17,1,22]:


~ f (c , t ) = (c ) 1 +
p=1

ap (t )Sp (c 2 ) ,

(9)

where (c ) = -3/2 exp(-c 2 ) is the dimensionless Maxwellian distribution and the first few Sonine polynomials read [1], S0 (x) = 1, S2 (x) = S3 (x) = 15 8 35 16 S1 (x) = 5 1 2 3 2

- x,

(10) (11)

- x + x2 , -
2 35 8 x+ 7 4 x2 - 1 6 x3 .

(12)

~ Provided that the expansion (9) converges, the Sonine coefficients ap completely determine the form of f (c , t ). According to the definition of temperature a1 (t ) = 0 [1], that is, the first non-trivial coefficient in the expansion (9) is a2 (t ). For the case of a constant the coefficient a2 has been found analytically [16,17] and the coefficient a3 analytically [18] and numerically [18­21]. The expansion (9) quantifies the deviation of the distribution function from the Maxwellian for the main part of ~ the velocity distribution, that is, for c 1; the high-velocity tail of f (c ) for c 1 is exponentially over-populated [15,17, 1,22] 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 [23], similarly the time dependence of the kinetic coefficients also changes [24,25]. 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 = const. [26,27], the clustering [27] and the vortex formation [28] occurs in a gas of viscoelastic particles only as a transient phenomena [29].


Author's personal copy
A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324 3317

Fig. 1. Evolution of the dimensionless granular temperature u( ) = T ( )/T (0), as it follows from the numerical solution of Eqs. (19), for different values of the dissipative parameter . Time is measured in the collision units 0 (see the text). For 1 all curves demonstrate the same slope, u( ) -5/3 (shown by the dotted line), in accordance with the theoretical prediction, Eq. (26). The inset illustrates the dependence (13) of the restitution coefficient on the dissipative parameter for u = 1 for collisions with the characteristic thermal velocity, that is for c12 · e = 1. The particular value of = 0.7, discussed below, corresponds to 0.3.

Based on the restitution coefficient given by Eq. (3) the theory of granular gases of viscoelastic particles has been developed, e.g. Ref. [1], where as in the case of a constant , only the first non-trivial Sonine coefficient a2 (t ) has been ~ taken into account. Although the evolution of the high-velocity tail of f (c ) has been analyzed [23], 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 a3 is of the same order of magnitude as the main coefficient a2 [18]; moreover, it was also shown that the amplitude of the ~ high-velocity tail of f (c ) and its contribution to the kinetic coefficients may be described quantitatively [30,31]. Finally, a new expression for the velocity-dependent restitution coefficient has been derived [32], 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 [32]. This affects the total dissipation at a collision, so that the revised reads [32]:


=1+
i=1

hi i/2 (2u(t ))i/20

c12 e

i/10

.

(13)

Here h1 = 0, h2 = -C1 , h3 = 0, h4 = C2 and the other numerical coefficients up to i = 20 are given in Ref. [32]. Fig. 1(inset) illustrates the dependence of the restitution coefficient (13) on the dissipative parameter for u = 1 for collisions with the characteristic thermal velocity, that is for c12 · e = 1. 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 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 f (c , t ), using the Sonine expansion up to the third-order term and analyze the amplitude and slope ~ of the high-energy tail of f (c , t ). 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 Section 2 we address the Sonine polynomial expansion and calculate the time-dependent Sonine coefficients, along with the granular temperature. In Section 3 the high-energy tail is analyzed, while Section 4 is devoted to the self-diffusion coefficient. Finally, in Section 5 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. Refs. [1,2]:

f (v , t ) = g2 ( )I (f , f ), (14) t where I (f , f ) is the collision integral. Generally, I depends on the two-particle distribution function f2 (v1 , v2 , r12 , t ). Within the hypothesis of molecular chaos, f2 (v1 , v2 , r12 , t ) = g2 ( )f (v1 , t )f (v2 , t ) (see e.g. Refs. [1]), the closed-form Eq. (14) for f (v , t ) is obtained; g2 ( ) denotes here a pair correlation function at a contact, which accounts for the increasing collision
frequency due to the effect of excluded volume [1].


Author's personal copy
3318 A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324

~ Eq. (14) yields the equation for the dimensionless distribution function f (c , t ) [1], µ
2

3

3+c

c

~ ~ f (c , t ) + B-1 f (c , t ) = ~(f , f ), I~~ t

(15)

where B = vT g2 ( ) 2 n and

µp = -

dc c p ~(f , f ), I~~

(16)

is the p-th moment of the dimensionless collision integral:

~(f , f ) = I~~

dc2

de (-c12 · e ) -c12 · e

1



2

~ ~ ~ ~ f (c1 , t )f (c2 , t ) - f (c1 , t )f (c2 , t ) .

(17)

(x) in the above equation is the Heaviside step-function and the dimensionless velocities c1 and c2 are the pre-collision
velocities in the so-called inverse collision, which results with c1 and c2 as the after-collision velocities, e.g. Ref. [1]. Eq. (15) is coupled to the equation for the granular temperature, e.g. Ref. [1]: dT dt

= - BT µ2 = - T ,
3

2

(18)

which defines the cooling coefficient = (2/3)Bµ2 . Multiplying Eq. (15) with c p , integrating over c and writing Eq. (18) in the dimensionless form, we obtain for p = 4, 6 [1] (similar equations have been used in Ref. [34]):

2 du = - µ2 u d 6
da
2

3 2



2 µ2 (1 + a2 ) - µ4 u d 3 15 d a3 u 2 2 2u = µ 2 ( 1 - a2 + a3 ) - µ 4 u + µ6 . d 5 105 2 =
2u



(19)

Here = t /c (0) is the dimensionless time, measured in the mean collision units, c-1 (0) = 4 g2 ( ) 2 n T (0)/m at initial time t = 0. The equations for the Sonine coefficients a2 and a3 in (19) are the first two of the infinite system of equations for p = 4, 6, 8, . . . [34,1]. The moments µ2 , µ4 , µ6 , . . . in these equations depend on all the Sonine coefficients ak , with k = 2, . . . , , that is, only the infinite set of equations are closed. To make the system tractable one has to truncate it. Following Refs. [18,33] we truncate the Sonine series at the third term, i.e. we approximate,





~ f (c , )

(c ) 1 + a2 ( )S2 (c 2 ) + a3 ( )S3 (c 2 ) .

(20)

Within this approximation µ2 , µ4 , µ6 in Eqs. (19) depend only on a2 , a3 and u, which makes the system (19) closed. The moments of the collisional integral µp were calculated analytically up to O ( 10 ), using the formula manipulation program, explained in detail in Ref. [1]. The complete expressions are rather cumbersome, therefore, we present here for illustration only linear approximations with respect to a2 , a3 and ( ) = (2u( ))1/10 :

µ2 = 0 1 + µ4 =

6 25

a2 +

2 125

a

3

( ) +
129 125 a2 - 179 1250 1079 100 a
3

2 (4a2 - a3 ) + 7
3

4 5

0 ( )
40717 1250 a2 - 39353 3125 a
3

(21)

45 µ6 = 3 2 15a2 - a
4

+ 30 ( )

+

,

where 0 = 2 2 21/10 (21/10) C1 6.485. 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 , a2 = a
20



+ a21 + · · · ,

a3 = a

30

+ a31 + · · · ,

u = u0 + u1 + · · · .

(22)

The zeroth-order solution (with a2 (0) = a3 (0) = 0) reads, a
20

= 0,

a

30

= 0,

u0 = (1 + /0 )-5/3 ,

(23)


Author's personal copy
A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324 3319

Fig. 2. Evolution of the second Sonine coefficient a2 ( ). At the first stage of the evolution |a2 ( )| increases, then reaches the maximum value |a2 max |, and eventually decays to zero. With the increasing dissipation parameter the maximum |a2 max | shifts to later time. The inset shows the asymptotic dependence of |a2 ( )| (full line) together with the analytical result (24) (dotted line) for = 0.01 in the log­log scale.

Fig. 3. Evolution of the third Sonine coefficient a3 ( ). The inset shows the asymptotic dependence of |a3 ( )| (full line) together with the analytical result (25) (dotted line) for = 0.01 in the log­log scale.

- where 0 1 = 26/5 (21/10) C1 /5 analytically:

0.55 . For the first-order solution only its asymptotics for may be obtained

a2 = -A2 ( /0 )-1/ a3 = -A3 ( /0 )-1/

6

6

157 A2 = 21/5 (21/10) C1 500 28 A3 = 21/5 (21/10) C1 500

0.44 0.08

(24) (25) (26)

u = ( /0 )-5/3 + q ( /0 )-11/6 , q = 2 C1
1 5

(16/5) (21/10) + 15625 (21/10)
2383

3.28.

Interestingly, the third Sonine coefficient a3 is of the same order of magnitude with respect to the small parameter , as the second Sonine coefficient a2 , albeit five times smaller. This conclusion is in a sharp contrast with the hypothesis suggested in Ref. [34], that the Sonine coefficients ak are of an ascending order of some small parameter , that is, ak k . 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 u( ) is plotted and in Figs. 2 and 3(insets). Fig. 1 demonstrates that the larger the dissipative parameter , the earlier the asymptotic behavior of u( ) is achieved. The numerical solution for the Sonine coefficients a2 ( ) and a3 ( ), shown in Figs. 2 and 3, corresponds to the initial Maxwellian distribution, a2 (0) = a3 (0) = 0. As it follows from the figures, the absolute values of both the coefficients initially increase, reach their maxima |a2 max | and |a3 max | 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 |a2 max | and |a3 max | first increase with the increasing dissipative parameter , then saturate at 0.15 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


Author's personal copy
3320 A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324

Fig. 4. Dependence of the second and third (inset) Sonine coefficients on the restitution coefficient for the case of a constant [1,17,18]. Note that corresponds to the vanishing a2 .

0.7

Fig. 5. The time dependence of the second Sonine coefficient a2 ( ), computed with a different number of terms in the expansion (13) for the restitution coefficient. The full lines correspond to the accuracy up to O ( 10 ), the dashed lines ­ up to O ( 19/2 ) and the dotted lines ­ up to O ( 9 ). Note that for = 0.33 the accuracy up to O ( 10 ) is insufficient to obtain a reliable convergence.

tendency -- the larger the dissipation parameter, the slower the gas evolution. Again we see that the third Sonine coefficient a3 is of the same order of magnitude as a2 (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 0 < < 0.15 [which corresponds to 0.85 < < 1, Fig. 1(inset)] one observes a fast relaxation, on a collision time scale, of the Sonine coefficients a2 , a3 to their maximal values |a2 max |, |a3 max |, roughly corresponding to the respective values for the constant . In a course of time, the granular temperature u( ) decreases and the effective restitution coefficient alters in accordance with decreasing dissipative parameter ( ) = (2u( ))1/10 , that is, the effective increases with time. For this interval of (0.85 < < 1) 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 a2 ( ) and a3 ( ). The larger the (i.e. the smaller the effective restitution coefficient), the larger the maxima |a2 max | and |a3 max |, in accordance with the dependencies a2 ( ) and a3 ( ) depicted in Fig. 4. Similarly, for 0.15 < < 0.30 [which corresponds to 0.7 < < 0.85, Fig. 1(inset)] the initial fast relaxation to the related values of |a2 max | and |a3 max | first takes place. Then the decreases of ( ) and the respective increase of the effective causes the increase of the absolute value of a2 ( ), until it reaches the maximum, corresponding to ( ) = 0.15 (or = 0.85), Fig. 4. The further decrease of ( ) leads to the according decay of a2 ( ), in agreement with the dependence of a2 ( ) shown in Fig. 2. This qualitatively explains the evolution of a2 ( ) and the saturation of its maximum |a2 max | for > 0.15. The qualitative behavior of the third Sonine coefficient may be explained analogously. It is interesting to note that the observed dependence of a2 ( ) with the new restitution coefficient (13) differs qualitatively for > 0.145 from that obtained previously for the old restitution coefficient (3). While in the latter case a positive bump at the initial time, 1 was detected for > 0.145 [23], in the former case the positive bump appears at the much larger > 0.3, Fig. 5. This is again in agreement with the behavior expected from the dependence of a2 ( ) for a constant : The coefficient a2 becomes positive for < 0.7, which corresponds to > 0.3, Fig. 1(inset). Note, however, that with increasing more and more terms in the expansion (13) are to be kept. While for = 0.27 it is suffice to keep terms up to 9 , for = 0.33 the Sonine coefficients a2 computed with the accuracy O ( 9 ), O ( 19/2 ) and O ( 10 ) noticeably differ, Fig. 5. Therefore we conclude that the revised restitution coefficient for a viscoelastic impact (13) may be accurately used up to


Author's personal copy
A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324 3321

= 0.3. 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 [18]. Similarly as in the case of constant [18], we expect that in the domain of convergence of the Sonine expansion, < 0.3, the magnitude of the next-order Sonine coefficients a4 , a5 , . . . is very small, so that an acceptable accuracy may be achieved with the use of the two coefficients, a2 and a3 only. 3. High-velocity tail of the velocity distribution function The expansion (9) refers to the main part of the velocity distribution, c 1, that is, to the velocities, close to the thermal one, vT . The high-velocity tail c 1 is however exponentially overpopulated [15]. It develops in a course of time, during a first few tens of collisions [30,31]. For viscoelastic particles the velocity distribution reads for c 1 [23]:

~ f (c , ) exp (- ( )c ) ,
where the function ( ) satisfies the equation [23], 1 3

(27)

+ µ2 B = B,

(28)

with B and µ2 defined previously. Using µ2 ( ), obtained by the formula manipulation program [1] (see the discussion above), we solve numerically Eq. (28) to obtain ( ). In the linear approximation, µ2 6.49 (2u( ))1/10 , the function ( ) has the form:

( ) = (b/ )(1 + /0 )

1 6

(29)

with b 1.129 [1]. Following [30] we neglect the transition region between the main part of the velocity distribution function, c 1, and its high-energy part, c 1 and write the distribution function as

~ f (c , ) = A( )c 2 exp(-c 2 ) 1 + a2 ( )S2 (c 2 ) + a3 ( )S3 (c 2 ) (c - c ) +B( )c 2 exp (- ( )c ) (c - c ).


(30)

~ The coefficients A( ), B( ) and the threshold velocity c , which separates the main and the tail part of f (c , ) can be obtained, using the normalization condition: ~ f (c )dc = 1
and the continuity condition for the function itself and its first derivative [30]: (31)

~ ~ f (c - 0, ) = f (c + 0, ) ~ ~ f (c - 0, ) f (c + 0, ) = . c c


(32) (33)

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

(2c - ) 1 + a2 S2 (c 2 ) + a3 S3 (c 2 ) = a2 2c 3 - 5c + a3 7c 3 - c 5 - 35 c 4 B = A exp(-c 2 + c ) 1 + a2 S2 (c 2 ) + a3 S3 (c 2 ) (34) 12 erf(c ) exp(c 2 ) - 12a2 c 5 + 30a2 c 3 + 35a3 c 3 - 28a3 c 5 + 4a3 c 7 - 24c A â exp(-c 2 ) + B 2 + c (2 + c ) exp(- c ) = 1 48 3 where ( ) is the solution of Eq. (28). To find the amplitudes A( ) and B( ) together with the threshold velocity c the system (34) was solved numerically. The asymptotic dependence of c on may be easily found if we take into account that a2 and a3 are of the same order of magnitude for 1, while ( ) 1 and c 1. Keeping in the first equation in (34) only the largest terms, yields 6 5 (2c - )a3 c -a3 c , which implies that
c ( )

( )/2 = (b/2 )(1 + /0 )1/6 .

(35)

The typical velocity distribution function, computed for = 50 and = 0.3 is shown in Fig. 6. The threshold velocity, c separating the main and the tail part of the velocity distribution reads in this case, c 4.31. The threshold velocity increases with decreasing inelasticity and shifts to larger values at a 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 f (c ) persists after its relaxation, in a gas of viscoelastic particles the velocity distribution function tends to a Maxwellian.


Author's personal copy
3322 A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324

Fig. 6. A typical velocity distribution function in a as a sum of velocity distribution function, obtained dash-dotted line). For comparison the Maxwellian respective quantity c = 3.77 for a gas of particles

granular gas of viscoelastic particles at = 50 for = 0.30 (Eq. (30), solid line), which is represented in the Sonine polynomial approximation (Eq. (20), dashed line), and the exponential function (Eq. (27), distribution is also shown (dotted line). The threshold velocity, c = 4.31, may be compared with the with a constant restitution coefficient 0.71.

Fig. 7. Evolution of the threshold velocity c in a gas of viscoelastic particles for = 0.3, 0.2, 0.15 (from top to bottom). The corresponding quantities for the case of a constant restitution coefficient are: c ( 0.71) = 3.77, c ( 0.79) = 3.65 and c ( 0.84) = 3.60 [30]. In a course of time the threshold velocity shifts to larger values, that is, the high-energy tail becomes less pronounced. The inset shows the asymptotic dependence of c (full line) together with the analytical result (35) (dotted line) for = 0.2 in the log­log scale.

4. Self-diffusion Self-diffusion is the only transport process which takes place in a granular the lack of macroscopic currents, a current of tagged particles, identical to the marked, may exist. Moreover, the diffusion coefficient D is directly related to via the Einstein relation, D/T , which approximately holds true for displacement of the tagged particles reads [1]
t

gas in a homogeneous cooling state: In spite of particles of the surrounding gas, but somehow the mobility coefficient of the tagged particles granular gases, e.g. [35,36]. The mean-square

[ r ( )]2 =

D(t )dt ,

(36)

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

D T -1 + Dv,ad = . T m Here = (2/3)Bµ2 is, as previously, the cooling coefficient and - T
-1 v,ad =

(37)

v ,ad 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 [1]:

1 6

vT (t )g2 ( ) 2 n â

dc1 c2

~ ~ de (-c12 · e) c12 · e f (c1 , t )f (c2 , t )(1 - )(c12 · e)2 .

(38)

~ With the approximation (20) for f (c , t )v ,ad was calculated up to 10 , using the formula manipulation program, described in Ref. [1]. Here we present for illustration purposes only its linear with respect to a2 , a3 and part: v,ad 3 1 = 1 + a2 + a E (0) 16 64
3

-

2 0

1 8

+

3 100

a2 +

1 500

a

3

( ).

(39)


Author's personal copy
A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324 3323

Fig. 8. Self-diffusion coefficient as a function of time, measured in the collision units 0 (see the text). The dissipative parameters, from the top to bottom are = 0.3 and = 0.1. D0 is the Enskog self-diffusion coefficient in a gas of elastic particles at the initial temperature T0 = T (0). The full lines correspond to the complete solution of Eq. (37), while the dashed lines depict the linear approximation to D, (42). Note that after about ten collisions, the complete and linear solutions become practically indistinguishable. With the dotted lines we show u1/2 ( ) = (T ( )/T0 )1/2 , which is equal to the ratio of two diffusion coefficients for elastic particles at temperatures T ( ) and T0 . The plotted ratio D( )/D0 clearly tends to u1/2 ( ), which manifests that the self-diffusion in ~ a gas of viscoelastic particles tends to that in a gas of elastic particles. The inset shows the relative deviation D = (D - D)/D of the diffusion coefficient ~ D computed with the use of both Sonine coefficients, a2 and a3 from the respective value, D, obtained with the use of a2 only ( = 0.3). Naturally, the location of maximum of D coincides with that of |a3 (t )|, Fig. 3.

Here E (0) is the Enskog velocity correlation time of elastic particles at initial time = 0:
- E 1 (0) =

8 3

T (0)
m

n 2 g2 .

(40)

Using the obtained expressions for µ2 and v ,ad up to O ( 10 ), we solve Eq. (37) numerically and compute the diffusion coefficient D( ). In Fig. 8 the ratio of D( )/D0 is plotted, where D0 = E (0)T (0)/m, (41) is the Enskog self-diffusion coefficient for a gas of elastic particles at initial temperature T (0). As it is clearly seen from the figure, the diffusion coefficient decreases with time in a way that the ratio D( )/D0 approaches u1/2 ( ) = (T ( )/T (0))1/2 -- the ratio of the two diffusion coefficients for elastic particles at the current temperature T ( ) and the initial temperature T (0). 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 D( ). Keeping only first-order terms in the expressions for µ2 , a2 , a3 and v ,ad and substituting them into Eq. (37) yields the diffusion coefficient: 4239 D = D0 u1/2 ( ) 1 + 16000 2
0





,

(42)

where, as previously, = (2u( ))1/10 and 0 6.485. Note that due to the account of the third Sonine coefficient in µ2 and v ,ad the linear approximation (42) for D( ) differs from the previously obtained result [1]. 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 solutions 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 the self-diffusion coefficient is rather small even for the large value of dissipation parameter = 0.3. 5. Conclusion We study the 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 us to model collisions with much larger dissipation, as compared to previously available results 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 [34], that the subsequent Sonine coefficients a2 , a3 , . . ., ak 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


Author's personal copy
3324 A.S. Bodrova, N.V. Brilliantov / Physica A 388 (2009) 3315­3324

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