- Email: [email protected]
Stationary state of thermostated inelastic hard spheres
Physica A 310 (2002) 308 – 324
www.elsevier.com/locate/physa
Stationary state of thermostated inelastic hard spheres Thierry Bibena; ∗ , Ph.A. Martinb , J. Piaseckic a Laboratoire
de spectrometrie physique, Universite Joseph Fourier, BP87 F-38402 Saint Martin d’H&eres cedex, France b Institut de physique th eorique, Ecole Polytechnique Federale, CH-1015, Lausanne, Switzerland c Institute of Theoretical Physics, University of Warsaw, PL-00 681 Warsaw, Poland Received 10 July 2001; received in revised form 11 February 2002
Abstract We consider a -uid composed of inelastic hard spheres moving in a thermostat modelled by a hard sphere gas. The losses of energy due to inelastic collisions are balanced by the energy transfer via elastic collisions from the thermostat particles. The resulting stationary state is analysed within the Boltzmann kinetic theory. A numerical iterative method permits to study the nature of deviations from the Gaussian state. Some analytic results are obtained for a one-dimensional c 2002 Elsevier Science B.V. All rights reserved. system. PACS: 05.20.−y; 51.10.+y; 81.−05.Rm
1. Introduction The dynamics of -uidized granular matter is often modelled by that of a system of inelastic hard spheres. The spheres have a mass M and a diameter . During binary collisions, the magnitude of the component of the relative velocity normal to the surface of spheres at contact is instantaneously reduced by a factor , 0 6 ¡ 1. The restitution coe7cient is supposed to be independent of the velocities (for general background and references see for instance Refs. [1–3]). More precisely, if (v1 ; v2 ) and (:v1 ; v: 2 ) denote the velocities of two inelastic spheres of mass M before and after the collision, they are related by M v: 1 + M v: 2 = M v1 + M v2 ; ∗
Corresponding author. Tel.: +33-4-76-51-48-94; fax: +33-4-76-63-54-95. E-mail address: [email protected] (T. Biben).
c 2002 Elsevier Science B.V. All rights reserved. 0378-4371/02/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 2 ) 0 0 7 7 9 - 3
(1)
T. Biben et al. / Physica A 310 (2002) 308 – 324
ˆ · v: 12 = −ˆ · v12 ; ˆ ⊥ · v: 12 = ˆ ⊥ · v12 ;
309
06¡1 ; (2)
where v12 = (v1 − v2 ), ˆ is a unit vector normal to the surface of the spheres at the point of impact, and ˆ ⊥ points in the tangent direction: ˆ ⊥ · ˆ = 0. The Grst relation is the conservation of the centre of mass momentum, whereas the second one says that the normal component of the relative velocity reverses its direction with a magnitude reduced by the factor (the tangent component remains unchanged). Solving (1), (2) for v: 1 and v: 2 gives 1+ 1+ v: 1 = v1 − (3) ˆ v: 2 = v2 + (ˆ · v12 )ˆ : (ˆ · v12 ); 2 2 The inverse relation is obtained by exchanging the roles of initial and Gnal velocities and replacing by −1 in (3). In a collision the amount M (1 − 2 )|ˆ · v12 |2 (4) 4 of kinetic energy is dissipated. Inelastic hard spheres evolving in isolation cool down because of the constant loss of kinetic energy at collisions. A much studied case is the homogeneous cooling state determined by the solution of the Boltzmann equation with collision law (2) (for kinetic equations with dissipative collisions, see for instance Refs. [4,5] and references in Ref. [3], the homogeneous cooling state is studied in Refs. [6 –8], see also the study of binary mixtures in Ref. [9]). However, if the spheres are allowed to interact with an external agency (a heat source) providing an energy -ux compensating for the internal dissipative loss, a non-equilibrium steady state may be created. A heat source can be introduced by means of an external random force (a white noise force) of intensity f0 , acting on the spheres. In the Boltzmann description, this amounts to supplementing the kinetic equation by a diJusion term. The resulting (homogeneous) steady state, determined for all values of f0 ¿ 0 and ¡ 1, has been analysed in Ref. [3, Chapter 5]. One of its characteristic features is an enhanced tail of the velocity distribution of the spheres (compared to equilibrium), namely the non-Gaussian behaviour B exp (−Av3=2 ) as v → ∞. As f0 → 0 or → 1, this state looses its stability since either the external source or the dissipation mechanism disappears. 1 Several other thermostats have been studied since, both stochastic and deterministic (see Ref. [10]) with various results for the asymptotic velocity distribution tail. While the Gaussian deterministic thermostat presented in Ref. [10] leads to an asymptotic tail of the form B exp (−Av), the non-Gaussian deterministic thermostat presented in the same reference leads to B exp (−Av2 ), a Gaussian tail. It is then clear on the basis of these results that the thermal equilibrium of granular matter strongly depends on the thermalizing process. It is then important to consider the precise physical interaction mechanism between a granular particle and the thermostat to determine the proper velocity distribution. The random force thermostat studied in Ref. [3] is a model for a stochastic thermostat where the particles of the -uid undergo collisions with the particles of the thermostat. 1
A limiting state can, however, be obtained if one lets f0 → 0 and → 1 jointly in an appropriate way.
310
T. Biben et al. / Physica A 310 (2002) 308 – 324
The collisions are supposed to be accounted for by the random force acting on the particles of the -uid. Although for purely elastic particles the random force thermostat is equivalent to the collisional thermostat, from the above discussion it is not clear that this property still holds for dissipative -uids. We show in this article that it is indeed not fully the case. In this work, the system of inelastic hard spheres (called hereafter the colloidal particles) is immersed in a heat bath maintained at a Gxed temperature T . This thermostat, introduced now at a microscopic level, consists itself of hard spheres with mass m and diameter undergoing purely elastic collisions with the colloidal particles. The dynamics of these colloidal particles together with their interactions with the bath particles is assumed to obey a Boltzmann equation. In Section 2, we write the equation determining the (homogeneous) steady state in a form amenable to numerical resolution by successive iterations. The stationary (non-Gaussian) velocity distribution of the colloidal particles is obtained numerically and discussed in Section 3 for various values of the restitution parameter and of the mass ratio m=M . Our Gndings diJer from those obtained for the above-described system of randomly driven inelastic spheres on two important points. Nothing singular occurs in the elastic limit since then the colloidal particles get simply thermalized at the same temperature as that of the bath; hence perturbation theory at weak inelasticity is possible. Secondly, the tail of the velocity distribution remains Gaussian for 0 6 ¡ 1, so that deviations from equilibrium occur only locally (i.e., deviations only appear in a Gnite region of the velocity space). In Section 4, we give analytic argument for it in the simpler setting of the one-dimensional system. Some concluding remarks are presented in Section 5.
2. Analysis of the Boltzmann equation for the stationary state Consider a homogeneous binary mixture of colloidal particles of mass M and diameter moving in an inGnite three-dimensional thermostat modelled by a hard sphere -uid at temperature T . The thermostat particles have mass m and their state is described by the one particle density nT T (v), where T is the Maxwell velocity distribution given by 3=2 m mv2 ; (5) exp − T (v) = 2kB T 2kB T where nT represents the number density of the thermostat particles, and kB is the Boltzmann constant. It is important to note that the system we consider does not represent a true binary mixture since we assume that the velocity distribution of the thermostat particles is not modiGed by the collisions with the colloidal particles (i.e., it is a true thermostat). The temperature of the thermostat particles has then to be maintained by external means. The collisions between the thermostat particles and the colloidal particles are supposed to be elastic. The velocity v1 of a sphere of mass M and diameter encountering a sphere of mass m and diameter moving with relative velocity (v12 = v1 − v2 ) is
T. Biben et al. / Physica A 310 (2002) 308 – 324
311
instantaneously changed into v1 − (v12 · ˆ )ˆ ;
(6)
where 2m : m+M Correspondingly, the postcollisional velocity of the thermostat particle is given by =
v2 + (2 − )(v12 · ˆ )ˆ :
(7)
Here the meaning of the unit vector ˆ is the same as in Eq. (2). We denote by nc and F(v) the number density and the stationary velocity distribution of the colloidal particles, respectively. Within the Boltzmann kinetic theory the stationary state nc F(v) satisGes the equation G(v|F; T ) = L(v|F; T )
(8)
representing the balance between the gain and the loss in the density of colloidal particles with velocity v per unit time resulting from binary collisions. Putting 2 + = nc nT (9) 2 we write the gain term in the form G(v1 |F; T ) = dv2 d ˆ (v12 · ˆ )(v12 · ˆ ) × F(v1 − (v12 · ˆ )ˆ )T (v2 + (2 − )(v12 · ) ˆ ) ˆ ! 1+ 1+ + 2 F v1 − ; (v12 · ˆ )ˆ F v2 + (v12 · ˆ )ˆ 2 2
(10)
where !=
4nc ; nT (1 + )2
=
(11)
and is the unit Heaviside step function. Physically ! represents the ratio of the two mean free paths: for “colloid–thermostat” and for “colloid–colloid” collisions. The term in (10) linear in F describes the eJect of thermalizing elastic collisions, whereas the nonlinear term represents the dissipative encounters between the colloidal particles. In accordance with the collision laws (1),(2), and (6), (7) the arguments of the distributions F and T in the gain term (10) are such that the corresponding postcollisional velocity of the colloidal particle equals v1 . In the loss term the precollisional velocity of the colloidal particle is v1 . We thus have L(v1 |F; T ) = LcT (v1 |F; T ) + !Lcc (v1 |F)
(12)
312
T. Biben et al. / Physica A 310 (2002) 308 – 324
with
LcT (v1 |F; T ) =
and
d ˆ (v12 · ˆ )(v12 · ˆ )F(v1 )T (v2 )
(13)
Lcc (v1 |F) =
dv2
d ˆ (v12 · ˆ )(v12 · ˆ )F(v1 )F(v2 ) :
dv2
(14)
The structure of the loss term indicates that the Boltzmann equation (8) can be rewritten in the form appropriate for developing an iterative method for constructing numerically its solution. Indeed, several schemes are possible. A natural choice is to extract F(v1 ) from the loss term (12), and to divide both sides of (8) by the term (15) D1 (v1 |T ) = dv2 d ˆ (v12 · ˆ )(v12 · ˆ ){T (v2 ) + !F(v2 )} which leads to F(v1 ) =
G(v1 |F; T ) : D1 (v1 |T )
(16)
An alternative formulation is obtained by suppressing F(v2 ) from the deGnition (15). Dividing both sides of (8) by the term (17) D2 (v1 |T ) = dv2 d ˆ (v12 · ˆ )(v12 · ˆ )T (v2 ) we get the relation G(v1 |F; T ) − !Lcc (v1 |F) F(v1 ) = : D2 (v1 |T )
(18)
It turns out, that starting with Maxwell’s distribution with temperature T the iterations of equation (18) yield numerically a well converging series of approximations which permit to determine the stationary state F(v1 ) with good accuracy. However, in order to make the numerical problem tractable, we have to simplify the right-hand side of (18) by reducing the number of integrations. To this end we consider the orthonormal basis composed of unit vectors ˆ ;
ˆ ⊥ 1 ;
ˆ ⊥ 2
(19)
and we introduce for integration purposes the components u; w1 ; w2 of the relative velocity v12 deGned by ⊥ v12 = uˆ + w1 ˆ ⊥ 1 + w2 ˆ 2 :
(20)
Moreover, using the spherical symmetry of the stationary state we put F(v) = (v2 );
v = |v| :
(21)
The gain term (10) takes then the form ∞ ∞ ∞ dw1 dw2 du u d ˆ (v12 + 2 u2 − 2(v1 · ˆ )u) −∞
−∞
0
⊥ ⊥ ⊥ ×T ([(v1 · ˆ ) + (1 − )u]ˆ + [(v1 · ˆ ⊥ 1 ) − w1 ]ˆ 1 + [(v1 · ˆ 2 ) − w2 ]ˆ 2 )
T. Biben et al. / Physica A 310 (2002) 308 – 324
! + 2 ×
v12
+
1+ u 2
2
1+ − u(v1 · ˆ )
313
2 1− 2 2 ⊥ ⊥ u + [(v1 · ˆ 1 ) − w1 ] + [(v1 · ˆ 2 ) − w2 ] : (v1 · ˆ ) + 2 (22)
In the term linear in F, owing to the factorized structure of the Maxwell distribution, one can readily perform the integration over variables (w1 ; w2 ). In the nonlinear term we put u1 = (v1 · ˆ ⊥ 1 ) − w1 ;
u2 = (v1 · ˆ ⊥ 2 ) − w2
and we perform the angular integration using polar coordinates in the (u1 ; u2 ) plane, setting w = u12 + u22 . Finally, introducing the integration variable ! equal to the cosine of the angle between the velocity v1 and the unit vector ˆ we eventually arrive at a representation of the gain term in the form (10) +1 ∞ d! du u (v12 + 2 u2 − 2v1 u!)T [v1 ! + (1 − )u] G = 2 −1
! + 2 ×
0
∞
0
v12
+
1+ u 2
2
1+ uv1 ! −
2 1− v1 ! + u +w 2
dw
containing the one-dimensional Maxwell distribution 1=2
m m T [v1 ! + (1 − )u] = exp − (v1 ! + (1 − )u)2 : 2kB T 2kB T The most complicated nonlinear term involves thus just three integrations. D2 (v1 |T ) and Lcc (v1 |F) can be reduced to a single integral using the relation: d (v ˆ 12 · ˆ )(v12 · ˆ ) = |v12 | from which integrals of the form I(v1 |f) ≡ dv2 d ˆ (v12 · ˆ )(v12 · ˆ )f(v2 ) ; where f is an arbitrary function, can be rewritten: +∞ +1 I(v1 |f) = 22 v22 dv2 f(v2 ) dy v12 + v22 − 2v1 v2 y 0
−1
which can be integrated with respect to y to give +∞ (v1 + v2 )3 − |v1 − v2 |3 2 I(v1 |f) = 22 dv2 v2 f(v2 ) : 3v1 v2 0
(23)
314
T. Biben et al. / Physica A 310 (2002) 308 – 324
Since from expressions (14) and (17) Lcc (v1 |F) = F(v1 )I(v1 |F) and D2 (v1 |T ) = I(v1 |T ), these functions take the form +∞ (v1 + v2 )3 − |v1 − v2 |3 2 2 v2 T (v2 ) (24) D2 = 2 dv2 3v1 v2 0 while Lcc = 22 F(v1 )
0
+∞
dv2
(v1 + v2 )3 − |v1 − v2 |3 2 v2 (v22 ) : 3v1 v2
(25)
The ratio of the derived expressions (23) – (25) and (24) is equal to the stationary distribution, and represents a very useful starting point for numerical solution of the Boltzmann equation. In the stationary state the cooling down through dissipative collisions is exactly balanced by the -ow of energy arriving via elastic collisions from the thermostat particles. The energy density qc absorbed per unit time by the colloidal particles is thus given by Mv12 F(v1 )F(v2 ) ˆ 12 · ˆ )(v12 · ˆ ) qc = n2c 2 dv1 dv2 d (v 2 1 1+ 1+ (v12 · ˆ )ˆ F v2 + (v12 · ˆ )ˆ − 2 F v1 − : (26) 2 2 By standard transformations of the Boltzmann collision operator one Gnds from (26) the relation M 2 2 nc (1 − 2 ) dv1 dv2 |v12 |3 F(v1 )F(v2 ) : qc = (27) 16 It represents the collisional average of the energy loss (4) in a single encounter. In the low inelasticity limit ( → 1), owing to the presence of the factor (1 − 2 ), one can evaluate qc by approximating F(v) with the thermostat Maxwell distribution T (v). Finally, we give here the detailed form which takes the Boltzmann equation (8) in a one-dimensional space. The gain and the loss terms are then given by 2m 2M G(v1 |F; T ) = dv2 |v12 | nc nT F v1 − v12 T v2 + v12 m+M m+M
n 2 1+ 1+ c + ; (28) F v1 − v12 F v2 + v12 2 2 L(v1 |F; T ) = dv2 |v12 |{nc nT F(v1 )T (v2 ) + n2c F(v1 )F(v2 )} : (29) The stationary state satisGes the constraint G(v1 |F; T ) = L(v1 |F; T ). 3. Numerical solution Although a direct iteration of (18), using expressions (23) – (25), converges within some few tens of iterations provided the mass ratio is not too small, using a simple
T. Biben et al. / Physica A 310 (2002) 308 – 324
315
mixing scheme allows to compute the solutions for any value of the parameters. If Fnin (v) is the input function at step n and Fnout (v) is the corresponding output using expression (18), the new input function can be constructed as a linear interpolation in between Fnin and Fnout : Fn+1 (v) = (1 − ')Fnin (v) + 'Fnout (v). With ' = 0:1 the convergence becomes very regular. The normalization of the velocity distribution is forced at each step. Integrals appearing in expressions (23) – (25) have been estimated using Simpson’s trapezoidal method. The number of grid points for velocity integrations has been Gxed to 200, with dv=0:1 kB T=M . Angular integrations (variable !, see Eq. (23)) have been performed using 81 points, which corresponds to d! = 0:025. The velocity distribution (v2 ) has been sampled on the velocity lattice, and values of for velocities which do not correspond to a lattice point have been estimated using a quadratic interpolation scheme. The three relevant parameters in this problem are: the restitution coe7cient which has been deGned in relation (2) and is directly related to the energy loss during a collision; the “density ratio” ! which has been deGned in relation (11) and contains a weak dependence in the size ratio of the particles; and Gnally the mass ratio m=M . Since we are interested in dissipative systems, varies between 0 and 1. The variation domain of ! is more subjective, since it can be varied from 0 to inGnity, but we decided to restrict the variation domain to the interval [0,1], which corresponds to the physical situation 6 1 and nc =nT ¡ 1 suitable for an application to colloidal suspensions. For the same reason the mass ratio will also be varied in the interval [0,1]. A measure of the non-Gaussian nature of the velocity distribution function is provided by moments +∞ 4v2 dvv n F(v) : (30) Mn = v n = 0
For a Gaussian distribution even moments are simply proportional to a power of the second moment. One can easily check that 5 M4 = v4 = (M2 )2 : (31) 3 A very sensitive test is then to plot ( = 3M4 =(5(M2 )2 ) − 1 which should be zero for a Gaussian distribution. Such a plot is presented in Fig. 1, for ! = 1, which corresponds to equal densities and equal diameters of the particles, and for various values of the restitution coe7cient and the mass ratio. First of all, we can notice that the distribution function deviates from a Gaussian as soon as = 1. Next, the amplitude of this deviation can reach more than 10% when = 0, m=M = 1 and ! = 1. It is interesting to comment on the amazing behaviour of ( when the mass ratio goes to zero. We can see in Fig. 1 that the deviation from a Gaussian distribution can be of positive or of negative sign, that the sign can change when the restitution coe7cient decreases. This kind of behaviour seems to be quite generic and has indeed been observed for other thermostats, as can be seen in Refs. [2,10]. We can see moreover that curves corresponding to diJerent values of the mass ratio can cross. Indeed, a vanishing mass ratio while ! remains Gxed is not of great relevance from the point of view of physics. Rather, we should
316
T. Biben et al. / Physica A 310 (2002) 308 – 324 0.14
m/M=1 m/M=0.1 m/M=0.01 m/M=0.001
0.12 0.10
ζ
0.08
ω=1
0.06 0.04 0.02 0.00
−0.02
0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 1. Deviation from a Gaussian distribution as a function of for ! = 1 and various mass ratios. 0.15
m/M=1, ω=1 m/M=0.1, ω=0.1 m/M=0.01, ω=0.01
ζ
0.10
0.05
0.00 0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 2. Deviation from a Gaussian distribution as a function of in the colloidal limit.
introduce the “colloidal” limit, in which the mass ratio m=M tends to zero while the mass densities of the individual particles, m=3 and M=3 , are held Gxed, together with their volume fractions. In this limit, m=M ∼ 3 and nc =nT ∼ 3 (see Eq. (11)), which implies that ! → 0 like the mass ratio. The corresponding variation of ( is reported in Fig. 2, which shows the deviation from a Gaussian distribution in the colloidal limit. We can now see that in this limit the deviation becomes vanishingly small, which simply illustrates the fact that a colloidal particle is mostly subject to collisions with the solvent molecules (the thermostat), and occasionally with other colloidal particles. Fig. 3 shows the variation of ( with the “density ratio” !, and again, not surprisingly, the deviation from a Gaussian distribution increases (roughly linearly) with !. Although the velocity distribution is not strictly Gaussian, it is usual to deGne an eJective temperature from the average kinetic energy of the dissipative particles: M 2 2 Te; = : (32) v 3kB 2
T. Biben et al. / Physica A 310 (2002) 308 – 324
317
0.15
α=0 α=0.5 α=0.9
ζ
0.10
m/M=1 0.05
0.00 0.0
0.2
0.4
ω
0.6
0.8
1.0
Fig. 3. Deviation from a Gaussian distribution as a function of ! for equal masses. 1.0
m/M=1 m/M=0.1 m/M=0.01
0.9
ω=1
Teff / T
0.8 0.7 0.6 0.5 0.4 0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 4. EJective temperature versus for ! = 1.
We can see in Figs. 4 and 5 the variations of this eJective temperature corresponding to the parameters of Figs. 1 and 2, namely for ! = 1, and in the colloidal limit. At Gxed value of ! (Fig. 4) the eJective temperature decreases with the mass ratio. On the contrary, in the colloidal limit the eJective temperature approaches the temperature of the thermostat. An interesting comparison is to plot the velocity distribution function F(v) together with the “eJective” Gaussian distribution (i.e., the Maxwell distribution corresponding to Te; ). Such a comparison is presented in Fig. 6 for equal masses and two very diJerent values of the restitution coe7cient =0 and 0:9. The two top curves are for = 0, the dashed line corresponds to the Gaussian distribution. We clearly see in this case the inaccuracy of an eJective Gaussian description. For = 0:9 (the two nearly superimposed bottom curves) such a description looks fairly reasonable. To emphasize this point, we have plotted in the insert the diJerence between the actual distribution function F(v) and the so-called eJective Gaussian distribution for both cases = 0 and 0:9. One should not be confused by the apparent discrepancy between the normalization of F(v) and the eJective Gaussian distribution for the case = 0: since these are three-dimensional distributions, normalization is Gxed by the condition
318
T. Biben et al. / Physica A 310 (2002) 308 – 324 1.0
Teff / T
0.9
0.8
0.7
m/M=1, ω=1 m/M=0.1, ω=0.1 m/M=0.01, ω=0.01
0.6
0.5 0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 5. EJective temperature versus in the colloidal limit. ω=1
0.15
α=0 α=0.9
Distribution function
0.04
α=0
0.02
0.10 0
−0.02
0.05
0
α=0.9
1
2
F(v) φT (v) eff
0.00 0.0
1.0
2.0
3.0
4.0
v* Fig. 6. Velocity distribution function for equal masses and two values of . v∗ = v M=kB T is the reduced velocity, and ! = 1. Te; (v) is the Maxwellian at temperature Te; . For = 0:9 F(v) and Te; (v) are hardly distinguishable. The insert shows the diJerence F(v) − Te; (v) for both values of .
4v2 dvF(v) = 1, and low velocities contributes only weakly. The corresponding velocity modulus distribution functions 4v2 F(v) and 4v2 Te; (v) have been plotted in Fig. 7, and still we can observe easily the non-Gaussian behaviour for = 0, and the quasi-Gaussian behaviour for = 0:9, for which the two curves are indistinguishable. To investigate further the Gaussian versus non-Gaussian behaviour, it is interesting to consider the high velocity asymptotic behaviour. If the distribution function behaves like F(v) ∼ A exp (−Bv n )
when v → ∞ ;
(33) n(v)= d lnd v
then the exponent n should be given by n = limv→∞ n(v) where ln{−ln [F(v)= F(0)]}. A plot of this function is given in Fig. 8 for m=M = 1 and various values of the restitution coe7cient , and we can see that although n(v) seems to go to n = 32 , which would correspond to the result reported in Refs. [2,3], it Gnally goes back to n = 2 at very high velocity. However, this high-velocity region is irrelevant from the physical point of view because the Gaussian distribution is vanishingly small therein
Modulus distribution function
T. Biben et al. / Physica A 310 (2002) 308 – 324
319
α=0
0.7
2
4πv F(v) 2 4πv φT (v)
0.6
eff
0.5 0.4
ω=1
0.3
α=0.9
0.2 0.1 0.0 0.0
1.0
2.0
3.0
4.0
v* Fig. 7. Velocity modulus distribution functions 4v2 F(v) (full line) and 4v2 Te; (v) (dashed line) for equal masses and = 0 and 0:9. v∗ = v M=kB T is the reduced velocity, and ! = 1 (see Fig. 6).
n(v)
2.0
1.8
α=0 DSMC α=0.5 α=0.9 α=0 iter. α=0.5 α=0.9 α=1.0
1.6
1.4 −0.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
ln(v)
Fig. 8. Asymptotic exponent for various values of ; m=M = 1 and ! = 1. v∗ = v M=kB T . Lines are obtained with the iterative method presented in this article while symbols correspond to the DSMC method.
(∼10−30 ). But, strictly speaking, the non-Gaussian nature of the velocity distribution is only local in the velocity space (i.e., it only appears in a Gnite region of the velocity space). At this level it is important to comment on the numerical procedure since Fig. 8 requires a very high accuracy in the computation of the distribution function. If one uses a Gaussian distribution function as an initial guess for F(v), the result presented in Fig. 8 could be an eJect of insu7cient convergence of the algorithm, the tail of the distribution being reminiscent of the initial Gaussian distribution. To prevent this eJect we used an initial function of the form: 3 F(v) = for v ¡ 1 and 0 otherwise : (34) 4 With this prescription, the Gaussian-like behaviour observed at very high velocity is a result of the converging scheme. Moreover, to test the reliability of the numerical procedure we have computed the Gaussian distribution function for =1 using the same
320
T. Biben et al. / Physica A 310 (2002) 308 – 324
algorithm. We can see in Fig. 8 that starting from the singular velocity distribution function (34) the numerical scheme nicely converged to the Gaussian distribution for all velocities. To emphasize the advantage of the iterative procedure we present in Fig. 8 a comparison with the direct simulation Monte Carlo method (DSMC) which is now a standard resolution method for Boltzmann problems [11]. The DSMC takes advantage of the mapping between the Boltzmann equation and a stochastic process, which is here simply deGned by the successive collisions of the colloidal particles. If we consider N colloidal particles with initial velocities chosen according to an arbitrary velocity distribution, the stationary velocity distribution can be constructed by generating the following sequence of collisions: at each step we choose at random whether there is a collision between two colloidal particles (with a probability !=(1 + !)) or between a colloidal particle and a particle from the thermostat (with a probability 1=(1 + !)). Next the colliding partners are chosen at random, and the collision is performed with a probability proportional to (v12 · ˆ )(v12 · ˆ ) where v12 is their relative velocity and ˆ their centre to centre direction. The velocity of a particle from the thermostat is chosen according to the distribution function T . This method allows the arbitrary initial velocity distribution to relax to the solution of the Boltzmann equation (8). Although this stochastic method is much simpler to implement than the iterative scheme presented in this article, it is plagued by an insu7cient sampling in the asymptotic high-velocity region, as can be seen in Fig. 8, where the velocity distribution reaches values as small as 10−6 . But we can see in this Ggure that both methods agree quite well. A similar study can be done for a one-dimensional system. In this case, when masses are equal, the equation for the stationary state can be put in the form (see Eqs. (28) and (29)) T (v)J (v) − F(v)I (v) 1+ 1− = !F(v)J (v) − ! du |u|F v + u F v− u ; 2 2
(35)
where F(v) is the unknown velocity distribution for colloidal particles, T (v) is the thermostat Maxwell distribution. J (v) and I (v) are integrals of the form +∞ J (v) = dw |v − w|F(w) ; (36) −∞
I (v) =
+∞
−∞
dw |v − w|T (w) :
(37)
One can see from these last two expressions that J (v) and I (v) behave like |v| when |v| → ∞. Dividing equation (35) by T (v)J (v) we get F(v) I (v) + ! = 1 + !,(v) ; (38) T (v) J (v)
T. Biben et al. / Physica A 310 (2002) 308 – 324
321
2
Γ(v)
α=1 α=0.99 α=0
1
ω=0.01 ω=0.5 ω=1 0
0
10
20
30
v* Fig. 9.Asymptotic behaviour of ,(v) when |v| → ∞ for various values of the dissipation and !. v∗ = v M=kB T (m=M = 1). For each value of ! decreases from bottom to top. We can note the dramatic change when changes from 0:99 to unity.
where
,(v) =
du |u|F(v + [(1 + )=2]u)F(v − [(1 − )=2]u) : T (v)J (v)
(39)
Since both J (v) and I (v) behave like |v| at high velocities, their ratio approaches unity in this limit. To go further we need to know the asymptotic behaviour of ,(v). When = 1 the solution of the problem is simply F(v) = T (v) so that from the deGnitions (36) and (39) we immediately deduce ,(v) ≡ 1. However, as soon as = 1 (and ! = 0), we can see in Fig. 9 that limv→∞ ,(v) = 0. As a consequence, we can conjecture from Eq. (38) that the asymptotic behaviour of the velocity distribution of the colloidal particles in one-dimension for equal masses is given by F(v) ∼
1 T (v) 1+!
when |v| → ∞
(40)
which implies that the asymptotic behaviour remains Gaussian.
4. Asymptotic Gaussian behaviour For simplicity, the study is restricted to the one-dimensional system with equal masses described by Eq. (35). We supplement the numerical result (40) by analytic argument which prove that under reasonable assumptions the asymptotic behaviour of the stationary distribution is necessarily Gaussian. For this it is now convenient to divide Eq. (35) by F(v)J (v): du |u|F(v + [(1 + )=2]u)F(v − [(1 − )=2]u) T (v) I (v) = +!−! : (41) F(v)J (v) F(v) J (v)
322
T. Biben et al. / Physica A 310 (2002) 308 – 324
We specify a large class of distributions by the conditions (i) – (iii) stated below and show that within this class, Eq. (41) imposes an asymptotic Gaussian character to F(v). Since F(v) = F(−v) we only treat the case of positive v. (i) dv|v|F(v) ¡ ∞. (ii) F(v) is monotonously decreasing for v large: there is v0 ¿ 0 such that F(v1 ) 6 F(v2 );
v1 ¿ v2 ¿ v0 :
The large velocity behaviour is precisely qualiGed by one or the other of the two following conditions: (iiia) For ' ¿ 1 and any !, limv→∞ F('v + !)=F(v) = 0. (iiib) limv→∞ v F(v) = C ¿ 0; ¿ 3. It is easily veriGed that condition (iiia) includes for instance all functions F(v) that behave like F(v) ∼ c1 exp (−c2 v );
v → ∞; ¿ 0; c1 ; c2 ¿ 0 :
(42)
Condition (iiib) covers the case of functions having an algebraic decay at inGnity. Under condition (i), both J (v) and I (v) are asymptotic to v as v → ∞ so that equation (41) implies
J (v) + J− (v) T (v) : (43) lim = 1 + ! − ! lim v→∞ F(v) v→∞ vF(v) In (43) the u-integral occurring in (41) has been splitted into the range of positive and negative u setting ∞ 1+ 1− duuF v + J (v) = u F v− u : (44) 2 2 0 Introducing the new integration variable w = [(1 − )=2]u − v ( ¡ 1), one gets 2 ∞
w J (v) 2 dw 1 + = vF(v) 1− v −v
F([2=(1 − )]v + [(1 + )=(1 − )]w) F(w) : (45) × F(v) In the integration domain, one has always [2=(1−)]v+[(1+)=(1−)]w ¿ v implying F([2=(1 − )]v + [(1 + )=(1 − )]w) 6 1; F(v)
v ¿ v0
(46)
by the monotonicity assumption (ii). Hence the integrand in Eq. (45) is bounded uniformly with respect to v by the integrable function (1 + |w|=v0 )F(w); v ¿ v0 . Using the condition (iiia) with ' = 2=(1 − ) ¿ 1 and ! = [(1 + )=(1 − )]w, the integrand
T. Biben et al. / Physica A 310 (2002) 308 – 324
323
in Eq. (45) tends pointwise to zero as v → ∞, so that by dominated convergence J (v) =0: v→∞ vF(v) lim
(47)
The same arguments apply to J− (v) with '=2=(1+) ¿ 1 ( ¡ 1). We thus conclude from (43) that for the functions belonging to the class (iiia) lim
v→∞
T (v) =1+! F(v)
(48)
namely relation (40) is proven when ¡ 1. We now show that an algebraic decay characterized by condition (iiib) is excluded. Using (ii) and (iiib) we obtain again from (45) by dominated convergence that −2 1− J (v) = ; (49) lim v→∞ vF(v) 2 where the normalization dvF(v) = 1 has been used. Thus, for class (iiib) Eq. (43) leads to the equality −2 −2 1− 1 1+ (50) =1+ : + 2 ! 2 But when ¿ 3 and 0 6 ¡ 1 the left-hand side of (50) does not exceed 1 so that the above equality cannot be satisGed, and functions in the class (iiib) cannot obey Eq. (41).
5. Concluding remarks We have investigated in this article the competition between a dissipative medium, modelled by a dissipative hard sphere gas, and a thermostat. The stationary state has been analysed within the Boltzmann kinetic theory, and we have used an iterative procedure to compute numerically the velocity distribution of the dissipative particles. From this analysis we can Grst conclude that the stationary velocity distribution is non-Gaussian. But we have shown that the non-Gaussian nature of the distribution is only local in the velocity space, and the asymptotic high-velocity tail remains Gaussian. To emphasize this point, a rigorous analysis of the one-dimensional case has been presented in the Gnal section. This result diJers from the analysis presented in Refs. [2,3] where the authors observe an overpopulation of the high-velocity region F(v) ∼ exp (−Av3=2 ). Whereas in the collision induced stationary state studied here deviations from the Gaussian distribution appear in the intermediate velocity region (in this sense they are local) the high energy part having a Gaussian shape, the DSMC studies of the homogeneous cooling state yielded an opposite situation: a Gaussian distribution for small velocities and deviations from it for large velocities in the form of a high-energy tail [8].
324
T. Biben et al. / Physica A 310 (2002) 308 – 324
It is important to note that the true (Gaussian) asymptotic regime appears at velocities larger than 10 kB T=M , which in practice corresponds to very low values of the distribution function. A Grst consequence is that this regime cannot be investigated using usual stochastic resolution methods such as DSMC. Another consequence is that it will be probably very di7cult to reach this regime in real experiments, even if the elastic component could be maintained at a constant temperature and thus represent a thermostat. Experimentalists could hopefully have access to the intermediate regime only, where the apparent asymptotic exponent n deGned by F(v) ∼ exp(−Av n ) will vary between 32 and 2 depending on the value of the restitution coe7cient and the mass ratio m=M . It must be noted that this behaviour is not speciGc to our thermostat, similar features have been observed [12], using a variant of the iterative method presented here, with the thermostat of Refs. [2,3]. References [1] H.M. Jaeger, S.R. Nagel, Rev. Mod. Phys. 68 (1966) 1259. [2] T.P.C. Van Noije, M.H. Ernst, Granular Matter 1 (1998) 57. [3] T.P.C. Van Noije, Kinetic and mesoscopic theories of granular -uids, Ph.D. Thesis, Institute for Theoretical Physics, University of Utrecht, 1999. [4] M.H. Ernst, J.R. Dorfman, W.R. Hoegy, J.M.J. van Leeuwen, Physica 45 (1969) 127. [5] J.J. Brey, J.W. Dufty, A. Santos, J. Stat. Phys. 87 (1997) 1051. [6] A. Goldshtein, M. Shapiro, J. Fluid Mech. 282 (1995) 75. [7] J.J. Brey, M.J. Ruiz-Montero, D. Cubero, Phys. Rev. E 54 (1996) 3664. [8] J.J. Brey, D. Cubero, M.J. Ruiz-Montero, Phys. Rev. E 59 (1999) 1256. [9] V. GarzTo, J.W. Dufty, Phys. Rev. E 60 (1999) 5706. [10] J.M. Montanero, A. Santos, Granular Matter 2 (2000) 53, cond-mat=0002323. [11] G. Bird, Molecular Gas Dynamics, Oxford University Press, New York, 1976; G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas -ows, Clarendon Press, Oxford, 1994. [12] A. Barrat, T. Biben, Z. Rcz, E. Trizac, F. van Wijland, J. Phys. A 35 (2002) 463.
www.elsevier.com/locate/physa
Stationary state of thermostated inelastic hard spheres Thierry Bibena; ∗ , Ph.A. Martinb , J. Piaseckic a Laboratoire
de spectrometrie physique, Universite Joseph Fourier, BP87 F-38402 Saint Martin d’H&eres cedex, France b Institut de physique th eorique, Ecole Polytechnique Federale, CH-1015, Lausanne, Switzerland c Institute of Theoretical Physics, University of Warsaw, PL-00 681 Warsaw, Poland Received 10 July 2001; received in revised form 11 February 2002
Abstract We consider a -uid composed of inelastic hard spheres moving in a thermostat modelled by a hard sphere gas. The losses of energy due to inelastic collisions are balanced by the energy transfer via elastic collisions from the thermostat particles. The resulting stationary state is analysed within the Boltzmann kinetic theory. A numerical iterative method permits to study the nature of deviations from the Gaussian state. Some analytic results are obtained for a one-dimensional c 2002 Elsevier Science B.V. All rights reserved. system. PACS: 05.20.−y; 51.10.+y; 81.−05.Rm
1. Introduction The dynamics of -uidized granular matter is often modelled by that of a system of inelastic hard spheres. The spheres have a mass M and a diameter . During binary collisions, the magnitude of the component of the relative velocity normal to the surface of spheres at contact is instantaneously reduced by a factor , 0 6 ¡ 1. The restitution coe7cient is supposed to be independent of the velocities (for general background and references see for instance Refs. [1–3]). More precisely, if (v1 ; v2 ) and (:v1 ; v: 2 ) denote the velocities of two inelastic spheres of mass M before and after the collision, they are related by M v: 1 + M v: 2 = M v1 + M v2 ; ∗
Corresponding author. Tel.: +33-4-76-51-48-94; fax: +33-4-76-63-54-95. E-mail address: [email protected] (T. Biben).
c 2002 Elsevier Science B.V. All rights reserved. 0378-4371/02/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 2 ) 0 0 7 7 9 - 3
(1)
T. Biben et al. / Physica A 310 (2002) 308 – 324
ˆ · v: 12 = −ˆ · v12 ; ˆ ⊥ · v: 12 = ˆ ⊥ · v12 ;
309
06¡1 ; (2)
where v12 = (v1 − v2 ), ˆ is a unit vector normal to the surface of the spheres at the point of impact, and ˆ ⊥ points in the tangent direction: ˆ ⊥ · ˆ = 0. The Grst relation is the conservation of the centre of mass momentum, whereas the second one says that the normal component of the relative velocity reverses its direction with a magnitude reduced by the factor (the tangent component remains unchanged). Solving (1), (2) for v: 1 and v: 2 gives 1+ 1+ v: 1 = v1 − (3) ˆ v: 2 = v2 + (ˆ · v12 )ˆ : (ˆ · v12 ); 2 2 The inverse relation is obtained by exchanging the roles of initial and Gnal velocities and replacing by −1 in (3). In a collision the amount M (1 − 2 )|ˆ · v12 |2 (4) 4 of kinetic energy is dissipated. Inelastic hard spheres evolving in isolation cool down because of the constant loss of kinetic energy at collisions. A much studied case is the homogeneous cooling state determined by the solution of the Boltzmann equation with collision law (2) (for kinetic equations with dissipative collisions, see for instance Refs. [4,5] and references in Ref. [3], the homogeneous cooling state is studied in Refs. [6 –8], see also the study of binary mixtures in Ref. [9]). However, if the spheres are allowed to interact with an external agency (a heat source) providing an energy -ux compensating for the internal dissipative loss, a non-equilibrium steady state may be created. A heat source can be introduced by means of an external random force (a white noise force) of intensity f0 , acting on the spheres. In the Boltzmann description, this amounts to supplementing the kinetic equation by a diJusion term. The resulting (homogeneous) steady state, determined for all values of f0 ¿ 0 and ¡ 1, has been analysed in Ref. [3, Chapter 5]. One of its characteristic features is an enhanced tail of the velocity distribution of the spheres (compared to equilibrium), namely the non-Gaussian behaviour B exp (−Av3=2 ) as v → ∞. As f0 → 0 or → 1, this state looses its stability since either the external source or the dissipation mechanism disappears. 1 Several other thermostats have been studied since, both stochastic and deterministic (see Ref. [10]) with various results for the asymptotic velocity distribution tail. While the Gaussian deterministic thermostat presented in Ref. [10] leads to an asymptotic tail of the form B exp (−Av), the non-Gaussian deterministic thermostat presented in the same reference leads to B exp (−Av2 ), a Gaussian tail. It is then clear on the basis of these results that the thermal equilibrium of granular matter strongly depends on the thermalizing process. It is then important to consider the precise physical interaction mechanism between a granular particle and the thermostat to determine the proper velocity distribution. The random force thermostat studied in Ref. [3] is a model for a stochastic thermostat where the particles of the -uid undergo collisions with the particles of the thermostat. 1
A limiting state can, however, be obtained if one lets f0 → 0 and → 1 jointly in an appropriate way.
310
T. Biben et al. / Physica A 310 (2002) 308 – 324
The collisions are supposed to be accounted for by the random force acting on the particles of the -uid. Although for purely elastic particles the random force thermostat is equivalent to the collisional thermostat, from the above discussion it is not clear that this property still holds for dissipative -uids. We show in this article that it is indeed not fully the case. In this work, the system of inelastic hard spheres (called hereafter the colloidal particles) is immersed in a heat bath maintained at a Gxed temperature T . This thermostat, introduced now at a microscopic level, consists itself of hard spheres with mass m and diameter undergoing purely elastic collisions with the colloidal particles. The dynamics of these colloidal particles together with their interactions with the bath particles is assumed to obey a Boltzmann equation. In Section 2, we write the equation determining the (homogeneous) steady state in a form amenable to numerical resolution by successive iterations. The stationary (non-Gaussian) velocity distribution of the colloidal particles is obtained numerically and discussed in Section 3 for various values of the restitution parameter and of the mass ratio m=M . Our Gndings diJer from those obtained for the above-described system of randomly driven inelastic spheres on two important points. Nothing singular occurs in the elastic limit since then the colloidal particles get simply thermalized at the same temperature as that of the bath; hence perturbation theory at weak inelasticity is possible. Secondly, the tail of the velocity distribution remains Gaussian for 0 6 ¡ 1, so that deviations from equilibrium occur only locally (i.e., deviations only appear in a Gnite region of the velocity space). In Section 4, we give analytic argument for it in the simpler setting of the one-dimensional system. Some concluding remarks are presented in Section 5.
2. Analysis of the Boltzmann equation for the stationary state Consider a homogeneous binary mixture of colloidal particles of mass M and diameter moving in an inGnite three-dimensional thermostat modelled by a hard sphere -uid at temperature T . The thermostat particles have mass m and their state is described by the one particle density nT T (v), where T is the Maxwell velocity distribution given by 3=2 m mv2 ; (5) exp − T (v) = 2kB T 2kB T where nT represents the number density of the thermostat particles, and kB is the Boltzmann constant. It is important to note that the system we consider does not represent a true binary mixture since we assume that the velocity distribution of the thermostat particles is not modiGed by the collisions with the colloidal particles (i.e., it is a true thermostat). The temperature of the thermostat particles has then to be maintained by external means. The collisions between the thermostat particles and the colloidal particles are supposed to be elastic. The velocity v1 of a sphere of mass M and diameter encountering a sphere of mass m and diameter moving with relative velocity (v12 = v1 − v2 ) is
T. Biben et al. / Physica A 310 (2002) 308 – 324
311
instantaneously changed into v1 − (v12 · ˆ )ˆ ;
(6)
where 2m : m+M Correspondingly, the postcollisional velocity of the thermostat particle is given by =
v2 + (2 − )(v12 · ˆ )ˆ :
(7)
Here the meaning of the unit vector ˆ is the same as in Eq. (2). We denote by nc and F(v) the number density and the stationary velocity distribution of the colloidal particles, respectively. Within the Boltzmann kinetic theory the stationary state nc F(v) satisGes the equation G(v|F; T ) = L(v|F; T )
(8)
representing the balance between the gain and the loss in the density of colloidal particles with velocity v per unit time resulting from binary collisions. Putting 2 + = nc nT (9) 2 we write the gain term in the form G(v1 |F; T ) = dv2 d ˆ (v12 · ˆ )(v12 · ˆ ) × F(v1 − (v12 · ˆ )ˆ )T (v2 + (2 − )(v12 · ) ˆ ) ˆ ! 1+ 1+ + 2 F v1 − ; (v12 · ˆ )ˆ F v2 + (v12 · ˆ )ˆ 2 2
(10)
where !=
4nc ; nT (1 + )2
=
(11)
and is the unit Heaviside step function. Physically ! represents the ratio of the two mean free paths: for “colloid–thermostat” and for “colloid–colloid” collisions. The term in (10) linear in F describes the eJect of thermalizing elastic collisions, whereas the nonlinear term represents the dissipative encounters between the colloidal particles. In accordance with the collision laws (1),(2), and (6), (7) the arguments of the distributions F and T in the gain term (10) are such that the corresponding postcollisional velocity of the colloidal particle equals v1 . In the loss term the precollisional velocity of the colloidal particle is v1 . We thus have L(v1 |F; T ) = LcT (v1 |F; T ) + !Lcc (v1 |F)
(12)
312
T. Biben et al. / Physica A 310 (2002) 308 – 324
with
LcT (v1 |F; T ) =
and
d ˆ (v12 · ˆ )(v12 · ˆ )F(v1 )T (v2 )
(13)
Lcc (v1 |F) =
dv2
d ˆ (v12 · ˆ )(v12 · ˆ )F(v1 )F(v2 ) :
dv2
(14)
The structure of the loss term indicates that the Boltzmann equation (8) can be rewritten in the form appropriate for developing an iterative method for constructing numerically its solution. Indeed, several schemes are possible. A natural choice is to extract F(v1 ) from the loss term (12), and to divide both sides of (8) by the term (15) D1 (v1 |T ) = dv2 d ˆ (v12 · ˆ )(v12 · ˆ ){T (v2 ) + !F(v2 )} which leads to F(v1 ) =
G(v1 |F; T ) : D1 (v1 |T )
(16)
An alternative formulation is obtained by suppressing F(v2 ) from the deGnition (15). Dividing both sides of (8) by the term (17) D2 (v1 |T ) = dv2 d ˆ (v12 · ˆ )(v12 · ˆ )T (v2 ) we get the relation G(v1 |F; T ) − !Lcc (v1 |F) F(v1 ) = : D2 (v1 |T )
(18)
It turns out, that starting with Maxwell’s distribution with temperature T the iterations of equation (18) yield numerically a well converging series of approximations which permit to determine the stationary state F(v1 ) with good accuracy. However, in order to make the numerical problem tractable, we have to simplify the right-hand side of (18) by reducing the number of integrations. To this end we consider the orthonormal basis composed of unit vectors ˆ ;
ˆ ⊥ 1 ;
ˆ ⊥ 2
(19)
and we introduce for integration purposes the components u; w1 ; w2 of the relative velocity v12 deGned by ⊥ v12 = uˆ + w1 ˆ ⊥ 1 + w2 ˆ 2 :
(20)
Moreover, using the spherical symmetry of the stationary state we put F(v) = (v2 );
v = |v| :
(21)
The gain term (10) takes then the form ∞ ∞ ∞ dw1 dw2 du u d ˆ (v12 + 2 u2 − 2(v1 · ˆ )u) −∞
−∞
0
⊥ ⊥ ⊥ ×T ([(v1 · ˆ ) + (1 − )u]ˆ + [(v1 · ˆ ⊥ 1 ) − w1 ]ˆ 1 + [(v1 · ˆ 2 ) − w2 ]ˆ 2 )
T. Biben et al. / Physica A 310 (2002) 308 – 324
! + 2 ×
v12
+
1+ u 2
2
1+ − u(v1 · ˆ )
313
2 1− 2 2 ⊥ ⊥ u + [(v1 · ˆ 1 ) − w1 ] + [(v1 · ˆ 2 ) − w2 ] : (v1 · ˆ ) + 2 (22)
In the term linear in F, owing to the factorized structure of the Maxwell distribution, one can readily perform the integration over variables (w1 ; w2 ). In the nonlinear term we put u1 = (v1 · ˆ ⊥ 1 ) − w1 ;
u2 = (v1 · ˆ ⊥ 2 ) − w2
and we perform the angular integration using polar coordinates in the (u1 ; u2 ) plane, setting w = u12 + u22 . Finally, introducing the integration variable ! equal to the cosine of the angle between the velocity v1 and the unit vector ˆ we eventually arrive at a representation of the gain term in the form (10) +1 ∞ d! du u (v12 + 2 u2 − 2v1 u!)T [v1 ! + (1 − )u] G = 2 −1
! + 2 ×
0
∞
0
v12
+
1+ u 2
2
1+ uv1 ! −
2 1− v1 ! + u +w 2
dw
containing the one-dimensional Maxwell distribution 1=2
m m T [v1 ! + (1 − )u] = exp − (v1 ! + (1 − )u)2 : 2kB T 2kB T The most complicated nonlinear term involves thus just three integrations. D2 (v1 |T ) and Lcc (v1 |F) can be reduced to a single integral using the relation: d (v ˆ 12 · ˆ )(v12 · ˆ ) = |v12 | from which integrals of the form I(v1 |f) ≡ dv2 d ˆ (v12 · ˆ )(v12 · ˆ )f(v2 ) ; where f is an arbitrary function, can be rewritten: +∞ +1 I(v1 |f) = 22 v22 dv2 f(v2 ) dy v12 + v22 − 2v1 v2 y 0
−1
which can be integrated with respect to y to give +∞ (v1 + v2 )3 − |v1 − v2 |3 2 I(v1 |f) = 22 dv2 v2 f(v2 ) : 3v1 v2 0
(23)
314
T. Biben et al. / Physica A 310 (2002) 308 – 324
Since from expressions (14) and (17) Lcc (v1 |F) = F(v1 )I(v1 |F) and D2 (v1 |T ) = I(v1 |T ), these functions take the form +∞ (v1 + v2 )3 − |v1 − v2 |3 2 2 v2 T (v2 ) (24) D2 = 2 dv2 3v1 v2 0 while Lcc = 22 F(v1 )
0
+∞
dv2
(v1 + v2 )3 − |v1 − v2 |3 2 v2 (v22 ) : 3v1 v2
(25)
The ratio of the derived expressions (23) – (25) and (24) is equal to the stationary distribution, and represents a very useful starting point for numerical solution of the Boltzmann equation. In the stationary state the cooling down through dissipative collisions is exactly balanced by the -ow of energy arriving via elastic collisions from the thermostat particles. The energy density qc absorbed per unit time by the colloidal particles is thus given by Mv12 F(v1 )F(v2 ) ˆ 12 · ˆ )(v12 · ˆ ) qc = n2c 2 dv1 dv2 d (v 2 1 1+ 1+ (v12 · ˆ )ˆ F v2 + (v12 · ˆ )ˆ − 2 F v1 − : (26) 2 2 By standard transformations of the Boltzmann collision operator one Gnds from (26) the relation M 2 2 nc (1 − 2 ) dv1 dv2 |v12 |3 F(v1 )F(v2 ) : qc = (27) 16 It represents the collisional average of the energy loss (4) in a single encounter. In the low inelasticity limit ( → 1), owing to the presence of the factor (1 − 2 ), one can evaluate qc by approximating F(v) with the thermostat Maxwell distribution T (v). Finally, we give here the detailed form which takes the Boltzmann equation (8) in a one-dimensional space. The gain and the loss terms are then given by 2m 2M G(v1 |F; T ) = dv2 |v12 | nc nT F v1 − v12 T v2 + v12 m+M m+M
n 2 1+ 1+ c + ; (28) F v1 − v12 F v2 + v12 2 2 L(v1 |F; T ) = dv2 |v12 |{nc nT F(v1 )T (v2 ) + n2c F(v1 )F(v2 )} : (29) The stationary state satisGes the constraint G(v1 |F; T ) = L(v1 |F; T ). 3. Numerical solution Although a direct iteration of (18), using expressions (23) – (25), converges within some few tens of iterations provided the mass ratio is not too small, using a simple
T. Biben et al. / Physica A 310 (2002) 308 – 324
315
mixing scheme allows to compute the solutions for any value of the parameters. If Fnin (v) is the input function at step n and Fnout (v) is the corresponding output using expression (18), the new input function can be constructed as a linear interpolation in between Fnin and Fnout : Fn+1 (v) = (1 − ')Fnin (v) + 'Fnout (v). With ' = 0:1 the convergence becomes very regular. The normalization of the velocity distribution is forced at each step. Integrals appearing in expressions (23) – (25) have been estimated using Simpson’s trapezoidal method. The number of grid points for velocity integrations has been Gxed to 200, with dv=0:1 kB T=M . Angular integrations (variable !, see Eq. (23)) have been performed using 81 points, which corresponds to d! = 0:025. The velocity distribution (v2 ) has been sampled on the velocity lattice, and values of for velocities which do not correspond to a lattice point have been estimated using a quadratic interpolation scheme. The three relevant parameters in this problem are: the restitution coe7cient which has been deGned in relation (2) and is directly related to the energy loss during a collision; the “density ratio” ! which has been deGned in relation (11) and contains a weak dependence in the size ratio of the particles; and Gnally the mass ratio m=M . Since we are interested in dissipative systems, varies between 0 and 1. The variation domain of ! is more subjective, since it can be varied from 0 to inGnity, but we decided to restrict the variation domain to the interval [0,1], which corresponds to the physical situation 6 1 and nc =nT ¡ 1 suitable for an application to colloidal suspensions. For the same reason the mass ratio will also be varied in the interval [0,1]. A measure of the non-Gaussian nature of the velocity distribution function is provided by moments +∞ 4v2 dvv n F(v) : (30) Mn = v n = 0
For a Gaussian distribution even moments are simply proportional to a power of the second moment. One can easily check that 5 M4 = v4 = (M2 )2 : (31) 3 A very sensitive test is then to plot ( = 3M4 =(5(M2 )2 ) − 1 which should be zero for a Gaussian distribution. Such a plot is presented in Fig. 1, for ! = 1, which corresponds to equal densities and equal diameters of the particles, and for various values of the restitution coe7cient and the mass ratio. First of all, we can notice that the distribution function deviates from a Gaussian as soon as = 1. Next, the amplitude of this deviation can reach more than 10% when = 0, m=M = 1 and ! = 1. It is interesting to comment on the amazing behaviour of ( when the mass ratio goes to zero. We can see in Fig. 1 that the deviation from a Gaussian distribution can be of positive or of negative sign, that the sign can change when the restitution coe7cient decreases. This kind of behaviour seems to be quite generic and has indeed been observed for other thermostats, as can be seen in Refs. [2,10]. We can see moreover that curves corresponding to diJerent values of the mass ratio can cross. Indeed, a vanishing mass ratio while ! remains Gxed is not of great relevance from the point of view of physics. Rather, we should
316
T. Biben et al. / Physica A 310 (2002) 308 – 324 0.14
m/M=1 m/M=0.1 m/M=0.01 m/M=0.001
0.12 0.10
ζ
0.08
ω=1
0.06 0.04 0.02 0.00
−0.02
0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 1. Deviation from a Gaussian distribution as a function of for ! = 1 and various mass ratios. 0.15
m/M=1, ω=1 m/M=0.1, ω=0.1 m/M=0.01, ω=0.01
ζ
0.10
0.05
0.00 0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 2. Deviation from a Gaussian distribution as a function of in the colloidal limit.
introduce the “colloidal” limit, in which the mass ratio m=M tends to zero while the mass densities of the individual particles, m=3 and M=3 , are held Gxed, together with their volume fractions. In this limit, m=M ∼ 3 and nc =nT ∼ 3 (see Eq. (11)), which implies that ! → 0 like the mass ratio. The corresponding variation of ( is reported in Fig. 2, which shows the deviation from a Gaussian distribution in the colloidal limit. We can now see that in this limit the deviation becomes vanishingly small, which simply illustrates the fact that a colloidal particle is mostly subject to collisions with the solvent molecules (the thermostat), and occasionally with other colloidal particles. Fig. 3 shows the variation of ( with the “density ratio” !, and again, not surprisingly, the deviation from a Gaussian distribution increases (roughly linearly) with !. Although the velocity distribution is not strictly Gaussian, it is usual to deGne an eJective temperature from the average kinetic energy of the dissipative particles: M 2 2 Te; = : (32) v 3kB 2
T. Biben et al. / Physica A 310 (2002) 308 – 324
317
0.15
α=0 α=0.5 α=0.9
ζ
0.10
m/M=1 0.05
0.00 0.0
0.2
0.4
ω
0.6
0.8
1.0
Fig. 3. Deviation from a Gaussian distribution as a function of ! for equal masses. 1.0
m/M=1 m/M=0.1 m/M=0.01
0.9
ω=1
Teff / T
0.8 0.7 0.6 0.5 0.4 0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 4. EJective temperature versus for ! = 1.
We can see in Figs. 4 and 5 the variations of this eJective temperature corresponding to the parameters of Figs. 1 and 2, namely for ! = 1, and in the colloidal limit. At Gxed value of ! (Fig. 4) the eJective temperature decreases with the mass ratio. On the contrary, in the colloidal limit the eJective temperature approaches the temperature of the thermostat. An interesting comparison is to plot the velocity distribution function F(v) together with the “eJective” Gaussian distribution (i.e., the Maxwell distribution corresponding to Te; ). Such a comparison is presented in Fig. 6 for equal masses and two very diJerent values of the restitution coe7cient =0 and 0:9. The two top curves are for = 0, the dashed line corresponds to the Gaussian distribution. We clearly see in this case the inaccuracy of an eJective Gaussian description. For = 0:9 (the two nearly superimposed bottom curves) such a description looks fairly reasonable. To emphasize this point, we have plotted in the insert the diJerence between the actual distribution function F(v) and the so-called eJective Gaussian distribution for both cases = 0 and 0:9. One should not be confused by the apparent discrepancy between the normalization of F(v) and the eJective Gaussian distribution for the case = 0: since these are three-dimensional distributions, normalization is Gxed by the condition
318
T. Biben et al. / Physica A 310 (2002) 308 – 324 1.0
Teff / T
0.9
0.8
0.7
m/M=1, ω=1 m/M=0.1, ω=0.1 m/M=0.01, ω=0.01
0.6
0.5 0.0
0.2
0.4
α
0.6
0.8
1.0
Fig. 5. EJective temperature versus in the colloidal limit. ω=1
0.15
α=0 α=0.9
Distribution function
0.04
α=0
0.02
0.10 0
−0.02
0.05
0
α=0.9
1
2
F(v) φT (v) eff
0.00 0.0
1.0
2.0
3.0
4.0
v* Fig. 6. Velocity distribution function for equal masses and two values of . v∗ = v M=kB T is the reduced velocity, and ! = 1. Te; (v) is the Maxwellian at temperature Te; . For = 0:9 F(v) and Te; (v) are hardly distinguishable. The insert shows the diJerence F(v) − Te; (v) for both values of .
4v2 dvF(v) = 1, and low velocities contributes only weakly. The corresponding velocity modulus distribution functions 4v2 F(v) and 4v2 Te; (v) have been plotted in Fig. 7, and still we can observe easily the non-Gaussian behaviour for = 0, and the quasi-Gaussian behaviour for = 0:9, for which the two curves are indistinguishable. To investigate further the Gaussian versus non-Gaussian behaviour, it is interesting to consider the high velocity asymptotic behaviour. If the distribution function behaves like F(v) ∼ A exp (−Bv n )
when v → ∞ ;
(33) n(v)= d lnd v
then the exponent n should be given by n = limv→∞ n(v) where ln{−ln [F(v)= F(0)]}. A plot of this function is given in Fig. 8 for m=M = 1 and various values of the restitution coe7cient , and we can see that although n(v) seems to go to n = 32 , which would correspond to the result reported in Refs. [2,3], it Gnally goes back to n = 2 at very high velocity. However, this high-velocity region is irrelevant from the physical point of view because the Gaussian distribution is vanishingly small therein
Modulus distribution function
T. Biben et al. / Physica A 310 (2002) 308 – 324
319
α=0
0.7
2
4πv F(v) 2 4πv φT (v)
0.6
eff
0.5 0.4
ω=1
0.3
α=0.9
0.2 0.1 0.0 0.0
1.0
2.0
3.0
4.0
v* Fig. 7. Velocity modulus distribution functions 4v2 F(v) (full line) and 4v2 Te; (v) (dashed line) for equal masses and = 0 and 0:9. v∗ = v M=kB T is the reduced velocity, and ! = 1 (see Fig. 6).
n(v)
2.0
1.8
α=0 DSMC α=0.5 α=0.9 α=0 iter. α=0.5 α=0.9 α=1.0
1.6
1.4 −0.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
ln(v)
Fig. 8. Asymptotic exponent for various values of ; m=M = 1 and ! = 1. v∗ = v M=kB T . Lines are obtained with the iterative method presented in this article while symbols correspond to the DSMC method.
(∼10−30 ). But, strictly speaking, the non-Gaussian nature of the velocity distribution is only local in the velocity space (i.e., it only appears in a Gnite region of the velocity space). At this level it is important to comment on the numerical procedure since Fig. 8 requires a very high accuracy in the computation of the distribution function. If one uses a Gaussian distribution function as an initial guess for F(v), the result presented in Fig. 8 could be an eJect of insu7cient convergence of the algorithm, the tail of the distribution being reminiscent of the initial Gaussian distribution. To prevent this eJect we used an initial function of the form: 3 F(v) = for v ¡ 1 and 0 otherwise : (34) 4 With this prescription, the Gaussian-like behaviour observed at very high velocity is a result of the converging scheme. Moreover, to test the reliability of the numerical procedure we have computed the Gaussian distribution function for =1 using the same
320
T. Biben et al. / Physica A 310 (2002) 308 – 324
algorithm. We can see in Fig. 8 that starting from the singular velocity distribution function (34) the numerical scheme nicely converged to the Gaussian distribution for all velocities. To emphasize the advantage of the iterative procedure we present in Fig. 8 a comparison with the direct simulation Monte Carlo method (DSMC) which is now a standard resolution method for Boltzmann problems [11]. The DSMC takes advantage of the mapping between the Boltzmann equation and a stochastic process, which is here simply deGned by the successive collisions of the colloidal particles. If we consider N colloidal particles with initial velocities chosen according to an arbitrary velocity distribution, the stationary velocity distribution can be constructed by generating the following sequence of collisions: at each step we choose at random whether there is a collision between two colloidal particles (with a probability !=(1 + !)) or between a colloidal particle and a particle from the thermostat (with a probability 1=(1 + !)). Next the colliding partners are chosen at random, and the collision is performed with a probability proportional to (v12 · ˆ )(v12 · ˆ ) where v12 is their relative velocity and ˆ their centre to centre direction. The velocity of a particle from the thermostat is chosen according to the distribution function T . This method allows the arbitrary initial velocity distribution to relax to the solution of the Boltzmann equation (8). Although this stochastic method is much simpler to implement than the iterative scheme presented in this article, it is plagued by an insu7cient sampling in the asymptotic high-velocity region, as can be seen in Fig. 8, where the velocity distribution reaches values as small as 10−6 . But we can see in this Ggure that both methods agree quite well. A similar study can be done for a one-dimensional system. In this case, when masses are equal, the equation for the stationary state can be put in the form (see Eqs. (28) and (29)) T (v)J (v) − F(v)I (v) 1+ 1− = !F(v)J (v) − ! du |u|F v + u F v− u ; 2 2
(35)
where F(v) is the unknown velocity distribution for colloidal particles, T (v) is the thermostat Maxwell distribution. J (v) and I (v) are integrals of the form +∞ J (v) = dw |v − w|F(w) ; (36) −∞
I (v) =
+∞
−∞
dw |v − w|T (w) :
(37)
One can see from these last two expressions that J (v) and I (v) behave like |v| when |v| → ∞. Dividing equation (35) by T (v)J (v) we get F(v) I (v) + ! = 1 + !,(v) ; (38) T (v) J (v)
T. Biben et al. / Physica A 310 (2002) 308 – 324
321
2
Γ(v)
α=1 α=0.99 α=0
1
ω=0.01 ω=0.5 ω=1 0
0
10
20
30
v* Fig. 9.Asymptotic behaviour of ,(v) when |v| → ∞ for various values of the dissipation and !. v∗ = v M=kB T (m=M = 1). For each value of ! decreases from bottom to top. We can note the dramatic change when changes from 0:99 to unity.
where
,(v) =
du |u|F(v + [(1 + )=2]u)F(v − [(1 − )=2]u) : T (v)J (v)
(39)
Since both J (v) and I (v) behave like |v| at high velocities, their ratio approaches unity in this limit. To go further we need to know the asymptotic behaviour of ,(v). When = 1 the solution of the problem is simply F(v) = T (v) so that from the deGnitions (36) and (39) we immediately deduce ,(v) ≡ 1. However, as soon as = 1 (and ! = 0), we can see in Fig. 9 that limv→∞ ,(v) = 0. As a consequence, we can conjecture from Eq. (38) that the asymptotic behaviour of the velocity distribution of the colloidal particles in one-dimension for equal masses is given by F(v) ∼
1 T (v) 1+!
when |v| → ∞
(40)
which implies that the asymptotic behaviour remains Gaussian.
4. Asymptotic Gaussian behaviour For simplicity, the study is restricted to the one-dimensional system with equal masses described by Eq. (35). We supplement the numerical result (40) by analytic argument which prove that under reasonable assumptions the asymptotic behaviour of the stationary distribution is necessarily Gaussian. For this it is now convenient to divide Eq. (35) by F(v)J (v): du |u|F(v + [(1 + )=2]u)F(v − [(1 − )=2]u) T (v) I (v) = +!−! : (41) F(v)J (v) F(v) J (v)
322
T. Biben et al. / Physica A 310 (2002) 308 – 324
We specify a large class of distributions by the conditions (i) – (iii) stated below and show that within this class, Eq. (41) imposes an asymptotic Gaussian character to F(v). Since F(v) = F(−v) we only treat the case of positive v. (i) dv|v|F(v) ¡ ∞. (ii) F(v) is monotonously decreasing for v large: there is v0 ¿ 0 such that F(v1 ) 6 F(v2 );
v1 ¿ v2 ¿ v0 :
The large velocity behaviour is precisely qualiGed by one or the other of the two following conditions: (iiia) For ' ¿ 1 and any !, limv→∞ F('v + !)=F(v) = 0. (iiib) limv→∞ v F(v) = C ¿ 0; ¿ 3. It is easily veriGed that condition (iiia) includes for instance all functions F(v) that behave like F(v) ∼ c1 exp (−c2 v );
v → ∞; ¿ 0; c1 ; c2 ¿ 0 :
(42)
Condition (iiib) covers the case of functions having an algebraic decay at inGnity. Under condition (i), both J (v) and I (v) are asymptotic to v as v → ∞ so that equation (41) implies
J (v) + J− (v) T (v) : (43) lim = 1 + ! − ! lim v→∞ F(v) v→∞ vF(v) In (43) the u-integral occurring in (41) has been splitted into the range of positive and negative u setting ∞ 1+ 1− duuF v + J (v) = u F v− u : (44) 2 2 0 Introducing the new integration variable w = [(1 − )=2]u − v ( ¡ 1), one gets 2 ∞
w J (v) 2 dw 1 + = vF(v) 1− v −v
F([2=(1 − )]v + [(1 + )=(1 − )]w) F(w) : (45) × F(v) In the integration domain, one has always [2=(1−)]v+[(1+)=(1−)]w ¿ v implying F([2=(1 − )]v + [(1 + )=(1 − )]w) 6 1; F(v)
v ¿ v0
(46)
by the monotonicity assumption (ii). Hence the integrand in Eq. (45) is bounded uniformly with respect to v by the integrable function (1 + |w|=v0 )F(w); v ¿ v0 . Using the condition (iiia) with ' = 2=(1 − ) ¿ 1 and ! = [(1 + )=(1 − )]w, the integrand
T. Biben et al. / Physica A 310 (2002) 308 – 324
323
in Eq. (45) tends pointwise to zero as v → ∞, so that by dominated convergence J (v) =0: v→∞ vF(v) lim
(47)
The same arguments apply to J− (v) with '=2=(1+) ¿ 1 ( ¡ 1). We thus conclude from (43) that for the functions belonging to the class (iiia) lim
v→∞
T (v) =1+! F(v)
(48)
namely relation (40) is proven when ¡ 1. We now show that an algebraic decay characterized by condition (iiib) is excluded. Using (ii) and (iiib) we obtain again from (45) by dominated convergence that −2 1− J (v) = ; (49) lim v→∞ vF(v) 2 where the normalization dvF(v) = 1 has been used. Thus, for class (iiib) Eq. (43) leads to the equality −2 −2 1− 1 1+ (50) =1+ : + 2 ! 2 But when ¿ 3 and 0 6 ¡ 1 the left-hand side of (50) does not exceed 1 so that the above equality cannot be satisGed, and functions in the class (iiib) cannot obey Eq. (41).
5. Concluding remarks We have investigated in this article the competition between a dissipative medium, modelled by a dissipative hard sphere gas, and a thermostat. The stationary state has been analysed within the Boltzmann kinetic theory, and we have used an iterative procedure to compute numerically the velocity distribution of the dissipative particles. From this analysis we can Grst conclude that the stationary velocity distribution is non-Gaussian. But we have shown that the non-Gaussian nature of the distribution is only local in the velocity space, and the asymptotic high-velocity tail remains Gaussian. To emphasize this point, a rigorous analysis of the one-dimensional case has been presented in the Gnal section. This result diJers from the analysis presented in Refs. [2,3] where the authors observe an overpopulation of the high-velocity region F(v) ∼ exp (−Av3=2 ). Whereas in the collision induced stationary state studied here deviations from the Gaussian distribution appear in the intermediate velocity region (in this sense they are local) the high energy part having a Gaussian shape, the DSMC studies of the homogeneous cooling state yielded an opposite situation: a Gaussian distribution for small velocities and deviations from it for large velocities in the form of a high-energy tail [8].
324
T. Biben et al. / Physica A 310 (2002) 308 – 324
It is important to note that the true (Gaussian) asymptotic regime appears at velocities larger than 10 kB T=M , which in practice corresponds to very low values of the distribution function. A Grst consequence is that this regime cannot be investigated using usual stochastic resolution methods such as DSMC. Another consequence is that it will be probably very di7cult to reach this regime in real experiments, even if the elastic component could be maintained at a constant temperature and thus represent a thermostat. Experimentalists could hopefully have access to the intermediate regime only, where the apparent asymptotic exponent n deGned by F(v) ∼ exp(−Av n ) will vary between 32 and 2 depending on the value of the restitution coe7cient and the mass ratio m=M . It must be noted that this behaviour is not speciGc to our thermostat, similar features have been observed [12], using a variant of the iterative method presented here, with the thermostat of Refs. [2,3]. References [1] H.M. Jaeger, S.R. Nagel, Rev. Mod. Phys. 68 (1966) 1259. [2] T.P.C. Van Noije, M.H. Ernst, Granular Matter 1 (1998) 57. [3] T.P.C. Van Noije, Kinetic and mesoscopic theories of granular -uids, Ph.D. Thesis, Institute for Theoretical Physics, University of Utrecht, 1999. [4] M.H. Ernst, J.R. Dorfman, W.R. Hoegy, J.M.J. van Leeuwen, Physica 45 (1969) 127. [5] J.J. Brey, J.W. Dufty, A. Santos, J. Stat. Phys. 87 (1997) 1051. [6] A. Goldshtein, M. Shapiro, J. Fluid Mech. 282 (1995) 75. [7] J.J. Brey, M.J. Ruiz-Montero, D. Cubero, Phys. Rev. E 54 (1996) 3664. [8] J.J. Brey, D. Cubero, M.J. Ruiz-Montero, Phys. Rev. E 59 (1999) 1256. [9] V. GarzTo, J.W. Dufty, Phys. Rev. E 60 (1999) 5706. [10] J.M. Montanero, A. Santos, Granular Matter 2 (2000) 53, cond-mat=0002323. [11] G. Bird, Molecular Gas Dynamics, Oxford University Press, New York, 1976; G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas -ows, Clarendon Press, Oxford, 1994. [12] A. Barrat, T. Biben, Z. Rcz, E. Trizac, F. van Wijland, J. Phys. A 35 (2002) 463.