Granular Gas Cooling and Relaxation to the Steady State in Regard to the Overpopulated Tail of the Velocity Distribution

0.57 MB·English

by  Thorsten Poeschel

February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
7
0
0
2
International JournalofModernPhysicsC
n (cid:13)c WorldScientificPublishingCompany
a
J

4
2

]
h GRANULAR GAS COOLING AND RELAXATION TO THE
c STEADY STATE IN REGARD TO THE OVERPOPULATED TAIL
e
m OF THE VELOCITY DISTRIBUTION
-
t
a
t THORSTENPO¨SCHEL
s
. Charit´e, Augustenburger Platz, 10439 Berlin, Germany
t
a [email protected]
m
- NIKOLAIV.BRILLIANTOV
d
Institute of Physics, Universityof Potsdam, Am NeuenPalais 10, 14469 Potsdam, Germany
n
[email protected]
o
c
[ ARNOFORMELLA

Universidad de Vigo, Department of Computer Science, Edificio Polit´ecnico,
1
32004 Ourense, Spain
v
[email protected]
0
8
5 ReceivedFebruary6,2008
1 RevisedDayMonthYear
0
7 Wepresentauniversaldescriptionofthevelocitydistributionfunctionofgranulargases,
0 f(v), valid for both, small and intermediate velocities where v is close to the thermal
/ velocity and also for large v where the distribution function reveals an exponentially
t
a decaying tail. By means of large-scale Monte Carlo simulations and by kinetic theory
m we show that the deviation fromthe Maxwelldistribution inthe high-energy tail leads
to small but detectable variation of the cooling coefficient and to extraordinary large
-
d relaxationtime.
n
o Keywords:granulargases;velocitydistributionfunction;overpopulatedhigh-energytail.
c
PACSNos.:47.70.-n,51.10.+y
:
v
i
X 1. Introduction
r
Granulargasesarecharacterizedbyacertainlossofkineticenergyaccordingtothe
a
dissipativepropertiesofparticlecollisions.The particlevelocitiesbeforea collision,
′
~v , and after,~v , are related by the collision rule
1/2 1/2
1+ε
′
~v =~v (~v ~e)~e
1 1− 2 12·
(1)
1+ε
′
~v =~v + (~v ~e)~e
2 2 2 12·
with ~v ~v ~v and the unit vector ~e (~r ~r )/ ~r ~r at the moment of
12 1 2 1 2 1 2
≡ − ≡ − | − |
the collision.
1
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
2 T. P¨oschel, N.V. Brilliantov, and A. Formella
In the absence of external excitation, a homogeneously initialized granular gas
stayshomogeneousduringthefirstpartofitsevolution,calledhomogeneous cooling
state (HCS). In later stages of its evolution, hydrodynamics instabilities develop
leading to vortex formation, that is, spacial correlations of the vectorial velocity
1 2
field, and finally cluster formation due to a pressure instability. In this paper we
restrict ourselves to the homogeneous cooling state. It should be mentioned that
3
heated granular gases with a Gaussian thermostat are equivalent to the HCS ,
therefore, the presented results apply also for such systems.
Accordingtothesteadylossofenergy,thevelocitydistributionfunction,f(~v,τ),
depends on time τ (measured in units of collisions per particle). Its shape can be
described by the time-independent function f˜(~c) via4
n
f(~v,τ)= f˜(~c) (2)
v3(τ)
T
where
~v
~c and v (τ) 2T(τ). (3)
T
≡ v (τ) ≡
T
p
5
TheevolutionofthegranulartemperatureT(τ)isgivenbyHaff’slaw (withparticle
mass m =1),
i
dT
= 2γT, i.e., T(τ)=T(0)exp( 2γτ). (4)
dτ − −
Here γ is the cooling coefficient, addressed in detail in the next section.
2. Universal velocity distribution function
For small and intermediate velocities, ~v v , that is, ~c 1, the scaled velocity
T
| |≈ | |≈
distribution function f˜(~c) is close to a Gaussian with small deviations that can
be characterized by the first non-trivial coefficient, a , of a Sonine polynomials
2
expansion around the leading Gauss distribution,6,7
1 5 15
f˜(c)=π−3/2exp c2 1+a S c2 with S c2 = c4 c2+ (5)
2 2 2
− 1 − 2 8
and (cid:0) (cid:1)(cid:2) (cid:0) (cid:1)(cid:3) (cid:0) (cid:1)
16(1 ε) 1 2ε2
a (ε)= − − . (6)
2 81 17ε+30ε+302(1 ε)
− (cid:0) (cid:1)−
For large kinetic energy, c 1, the distribution function decays exponentially
≫
slow, that is, the distribution function is overpopulated with respect to the Gauss
distribution one would expect for a molecular gas,4,7
3π
f˜(c)=Be−bc with b= (7)
µ
2
with the second moment of the dimensionless collision integral
3
µ d~c c2I˜ f˜,f˜ =√2π 1 ε2 1+ a (8)
2 ≡− 1 1 − 16 2
Z (cid:16) (cid:17) (cid:0) (cid:1)(cid:20) (cid:21)
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
Granular gas cooling and relaxation to the steady state 3
and with
1
I˜ f˜,f˜ = d~c d~e Θ( ~c ~e) ~c ~e f˜(~c′′)f˜(~c′′) f˜(~c )f˜(~c ) (9)
2 − 12· | 12· | ε2 1 2 − 1 2
(cid:16) (cid:17) Z Z (cid:20) (cid:21)
′′
where~c are the pre-collisionalscaled velocities.
1/2
There is a transition region between the ranges of validity of the Sonine expan-
sion (c 1) and the expression for the tail (c 1) where the distribution function
≈ ≫ 8
was quantified in terms of complicated implicit expressions. For practical compu-
tations,however,a muchsimpler approachmaybe sufficient:Itwasshownrecently
that the simple combination of the shown Sonine expansion and the exponential
tail,
f˜(c)= Ac2e−c2 1+a2S2(c2) for c<c∗ (10)
( (cid:2) Bc2e−b(cid:3)c for c c∗
≥
with the parameters A, B, and the transition velocity c∗ agrees with large-scale
9
DSMCsimulationuptoaverygoodaccuracyintheentirerangeofvelocities. The
parameters A, B, and c∗ follow from the conditions of smoothness
f˜(c∗+0)=f˜(c∗ 0)
−
df˜ df˜ (11)
=
dc(cid:12) dc(cid:12)
(cid:12)c∗+0 (cid:12)c∗−0
(cid:12) (cid:12)
and normalization: (cid:12) (cid:12)
(cid:12) (cid:12)
b a2 2c3∗ 5c∗
c∗ = + − (12)
2 2[1(cid:0)+a2S2(c2∗(cid:1))]
A1 = k(bc3∗)[2+bc∗(2+bc∗)]e−bc∗+√4πErf(c∗)−c8∗ 4+a2c2∗ 2c2∗−5 e−c2∗ (13)
(cid:2) (cid:0) (cid:1)(cid:3)
B = Ak(c∗) (14)
with
k(c∗) e−c2∗+bc∗ 1+a2S2 c2∗ . (15)
≡
Solving the fifth order equation (12) nu(cid:2)merically f(cid:0)or(cid:1)c(cid:3)∗, we obtain A, k and finally
9
B.
3. Impact of the high-velocity tail on the coefficient of cooling
Bymeansofthedistributionfunction,Eq.(10),anditsparametersA,B,c∗ depend-
ingonthe coefficientofrestitutionε,wecanquantifytheimpactofthe exponential
tail on the cooling coefficient γ, which is one of the most important characteristics
10
ofagranulargas.Usingthestandardanalysis, onecancalculatethecoolingcoef-
ficient for the homogeneouscooling state,when the stationaryvelocity distribution
is achieved:
J
γ = 1 ε2 2 , (16)
− 12J
0
(cid:0) (cid:1)
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
4 T. P¨oschel, N.V. Brilliantov, and A. Formella
where the coefficients J depend on the velocity distribution function,
k
J = d~c d~c d~eΘ( ~c ~e)~c ~e k+1f˜(c )f˜(c ). (17)
k 1 2 12 12 1 2
− · | · |
Z Z
7
If we neglect the exponential tail, we obtain the energy decay rate
1 ε21+ 3 a (ε)
γ (ε)= − 16 2 . (18)
0 6 1 1 a (ε)
− 16 2
Complementary, we determine the cooling coefficient numerically, using Direct
Simulation Monte Carlo (DSMC)11,12. We initialized a granular gas of N = 108
particleswithvelocitiesdrawnfromaGaussdistributionatT(0)=1andsimulated
until the particle velocities reachedthe limit of the double precision number repre-
sentation, i.e., until T 10−23. For ε = 0.9 this corresponds to a total of 5 1010
≈ ×
collisions or 1,000 collisions per particle.
For different ε we recorded the temperature T (τ). Then γ was
DSMC DSMC
determined by fitting T (τ) for τ 1 to its asymptotic law, T
DSMC DSMC
≫ ∝
exp( 2γ τ), Eq. (4).
DSMC
−
Thedependenceofγ onthecoefficientofrestitutionεisshowninFig.1together
0
with the numerical result. The numerical results (see below) and the theoretical
resultsdisregardingthetail,Eq.(18),virtuallycoincide.Whenweplotthedifference
ofbothcurves,however,wefindasystematicdeviationbetweenthenumericalvalues
and the analytical result (Fig. 2).
Toestimatethis rathersmalldifference quantitatively,wewritethe distribution
function as a sum of two terms,
f˜(c)=f˜(c)+f˜(c), (19)
0 t
0.16 Dγ0S, EMqC. (17)
0.14
0.12
γ 0.1
2
0.08
0.06
0.04
0.2 0.4 0.6 0.8
ε
Fig. 1. The cooling coefficient γ(ε) as a function of the coefficient of restitution as obtained by
DirectSimulationMonteCarlo(DSMC)(points)andduetoEq.(18),γ0.
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
Granular gas cooling and relaxation to the steady state 5
10-3 DbeSsMt eCxponential fit
−2γDSMC10-4 0.003625 e-5.9899 ε
0
γ
2
10-5
0 0.2 0.4 0.6 0.8
ε
Fig. 2. Difference of the cooling coefficient due to Eq. (18) and the numerically values, (γ0−
γDSMC), over the coefficient of restitution. Since Eq. (18) disregards the exponential tail of the
distribution, the shown difference is due to the presence of the tail. The dashed line shows the
bestexponential fit.
that is, the main partf˜(c) (with c extended to infinity), and the pure overpopula-
0
tion of the tail, f˜(c),
t
f˜(c)=Ac2e−c2 1+a S c2 (20)
0 2 2
f˜t(c)= Bc2exp(cid:2)( bc) f˜(cid:0)0(c)(cid:1)(cid:3)Θ(c c∗) . (21)
− − −
h i
The product f˜(c )f˜(c ) in Eq. (17) reads then
1 2
f˜(c )f˜(c )=f˜(c )f˜(c )+f˜(c )f˜(c )+f˜(c )f˜(c )+f˜(c )f˜(c ). (22)
1 2 0 1 0 2 0 1 t 2 t 1 0 2 t 1 t 2
For large c∗ one can neglect the last term in the above equation. Moreover, in this
case we approximate f˜(c) Bexp( bc) and~c ~e ~c ~e, assuming c c . A
t 12 1 1 2
≈ − · ≈ · ≫
similar approximation may be applied for the opposite case, c c . Taking then
2 1
≫
into account the symmetry of the integrand in Eq. (17) with respect to ~c and~c ,
1 2
we finally obtain an estimate for J :
k
π
J = A2J(0)+2 d~c f˜(c ) d~c d~e Θ( ~c ~e) ~c ~e k+1f˜(c )
k 16 k 2 0 2 1 − 1· | 1· | t 1
π (cid:20)Z√π ∞(cid:21)Z Z π (23)
= A2J(0)+2A 4π ck+3Be−bc1dc 2π sinθ cosθ k+1 ,
16 k 4 Zc∗ 1 1 Zπ2 | |
where J(0) is the value of the coefficient J for the distribution function f˜(c) with
k k 0
A in Eq. (20) equal to 4/√π, that is, for the case when the exponential tail is
disregarded. In particular, we need
1
J(0) =2√2π 1 a
0 − 16 2
(cid:18) (cid:19) (24)
3
J(0) =4√2π 1+ a .
2 16 2
(cid:18) (cid:19)
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
6 T. P¨oschel, N.V. Brilliantov, and A. Formella
14 2(γ0-γDSMC), numerics
2(γ-γ), theory, Eqs. (17,25)
0
12
]
0
010
0
0,
1 8
x
[
γ 6
∆
2
4
2
0 0.2 0.4
ε
Fig.3. Influenceofthetailonthecoolingcoefficient.ThesymbolsshowthesamedataasinFig.
2(butinlinearscale).Thefulllineshows2(γ0−γ)asgivenbyEqs.(18)and(26).
Performing the integration in Eq. (23) we obtain
π √π 2AB 8π2
Jk = 16A2Jk(0)+ 4 (k+2)bk+4Γ(k+4,bc∗), (25)
where Γ(x,y) is the incomplete gamma-function. Inserting into Eq. (16) yields the
cooling rate,
γ = (1−ε2) 1+ 136a2 Ab6+2π√2BΓ(6,bc∗) . (26)
6 (cid:0)1− 116a2 (cid:1)Ab6+8π√2Bb2Γ(4,bc∗)
Naturally, for B = 0, corres(cid:0)ponding t(cid:1)o absence of the tail (see Eq. (10)), Eq. (26)
reduces to the previous relation Eq. (18). In Fig. 3 the difference (γ γ) which
0
−
quantifies the contribution to the cooling coefficient from the tail is compared with
the numerical results.
Thescalinglaw,γ γ exp( 6ε),showninFig.2asanumericalresultis,
0 DSMC
− ∝ −
however,difficulttoconfirmanalytically.The approximateanalyticaltheoryshown
in Fig. 3 agrees with the numerical results only qualitatively.
4. Slow relaxation of the velocity distribution
In the above analysis we assumed that the velocity distribution function f˜(c) had
already relaxed to its steady state scaling shape. In our numerical experiments
we observed, however, that the relaxation to the stationary form occurs extremely
slowly, as compared to the relaxation of molecular gases to equilibrium distribu-
tion. To study this retardedrelaxationquantitatively,we use the coefficient γ
DSMC
described above, and define the temperature T (τ) exp( 2γ τ). By defini-
fit DSMC
∝ −
tion, for τ 1 one obtains T T since γ was determined as the best
DSMC fit DSMC
≫ ≈
exponential fit to T (τ) for τ 1. Hence, the quantity 1 T (τ)/T (τ)
DSMC DSMC fit
≫ −
characterizes the relaxation of the distribution function to its stationary form. We
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
Granular gas cooling and relaxation to the steady state 7
observethattherelaxationtime dependssensitivelyonthecoefficientofrestitution
ε.
Figures4-7illustratetherelaxationkineticsfordifferentvaluesofthecoefficient
ofrestitutionε.InFigs.4and5 the initialrelaxationis shown.We observeainitial
quick relaxation with characteristic time of 3-5 collisions per particles, similar as
for molecular gases.
For values ε<0.7 the simulatedtemperature approachesHaff’s law frombelow
(Fig.4),whereasforε>0.7itapproachesHaff’s lawfromabove(notethe different
labeling of the vertical axis in these figures). For ε=0.7 the initial relaxation rate
vanishes. This may be explained by the fact that at ε 0.7 the value of the second
≈
100
Tfit exp(-8/15 τ)
/ MC ε=0.1
S
D 10
T
-
1 ε=0.6 ε=0.3
ε=0.4
1
0 1 2 3 4
τ
Fig. 4. Short-time relaxation of the temperature decay to Haff’s law, Eq. (4) for a system of
N = 108 particles for different coefficients of restitution, ε= 0.1, ε= 0.3, ε = 0.4, and ε= 0.6.
Theinitialrelaxationoccursonthetimescale1−TDSMC/Tfit∝exp(−8/15τ).
10
ε=0.8
/ TDSMCfit ε=0.9 ε=0.7~exp(-0.384 τ)
T
+ const
1
- ~exp(-1/2 τ
)
1
0 1 2 3 4
τ
Fig. 5. Same as Fig. 4 but for ε = 0.7, ε = 0.8, and ε = 0.9. Note that the vertical axis has
oppositesignastheaxisinFig.4.
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
8 T. P¨oschel, N.V. Brilliantov, and A. Formella
Sonine coefficient changes its sign, see Eq. (6). That is, for ε . 0.7 the stationary
velocitydistributionf˜isbenttowardslowervelocitiesascomparedwiththeMaxwell
distribution while for ε&0.7 the distribution is bent towardshigher velocities. For
ε=0.7 the second Sonine coefficient is very small, a (0.7) 0, therefore, for c 1
2
≈ ∼
the distribution function is very close to the Maxwell distribution and, thus, there
is no initial relaxation. These arguments prove that the initial relaxation shown
in Figs. 4 and 5 corresponds to the relaxation of the main part of the velocity
distribution, c 1, whose deviation from the Maxwell distribution is described by
∼
the second Sonine coefficient a .
2
For small coefficient of restitution, ε 0.6, the initial slope of the relaxation
≤
curves is almostindependent of ε. Its value is surprisinglyclose to the slope 8/15
−
ofthe relaxationcurveofthe secondSonine coefficientin thelimit ofalmost elastic
10
particles, ε. 1. Also surprisingly, for larger values of the restitution coefficients
ε > 0.7, which correspond to negative values of the second Sonine coefficient, the
slope becomes smaller and depends noticeably on ε. For almost elastic particles,
ε = 0.9 we find the initial slope 1/2. Presently we do not have an explanation for
this behavior of the initial relaxation.
In Figs. 6 and 7 we present the complete relaxation to the steady state. The
completerelaxationrequiresmuchlongertimeascomparedtomoleculargases.This
isrelatedtotherelativelyslowformationoftheexponentialtail.Therelaxationtime
dependsonthecoefficientofrestitutionasshowninFig.6.Thelargerthecoefficient
ofrestitutionthelongerittakestoformthecompletevelocitydistributionincluding
the exponential tail.
To explain this effect we use the equation for the relaxation of the velocity
distribution function to its scaling form,10,13
µ ∂ ∂
2 3+c f˜(c,τ)+J f˜(c,τ)=I˜ f˜,f˜ (27)
0
3 ∂c ∂τ
(cid:18) (cid:19)
(cid:16) (cid:17)
ε=0.1
0
ε=0.1
Tfit -2 ε=0.3
/ C
M
-4
S
TD ε=0.3
-
1 -6 ε=0.4
-8
0 10 20 30 40 50
τ
Fig.6. Long-timerelaxationofthetemperaturedecayrateforε=0.1,ε=0.3,andε=0.4.
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
Granular gas cooling and relaxation to the steady state 9
where I˜(f˜,f˜) is the reduced collision integral, Eq. (9), and J is given in Eq. (17).
0
As usual we neglect the incoming term for c 14,7, leading to the approximation
≫
I˜ f˜,f˜ πcf˜(c) for c 1. (28)
≈− ≫
(cid:16) (cid:17)
With the Ansatz f˜(c,τ)=Bexp[ w(τ)c] we recast Eq. (27) into
−
dw µ π
2
+ w = for c 1 (29)
dτ 3J J ≫
0 0
with the solution
τ
w(τ)=b+(1 b)exp , (30)
− −τ (ε)
(cid:20) 0 (cid:21)
where b=3π/µ coincideswith Eq.(7)andτ−1(ε)=µ /3J .Neglecting a ,which
2 0 2 0 2
quantifies(small)deviationsofthemainpartofthedistributionwithrespecttothe
Maxwellian, and the contribution from the tail, Eq. (24) yields J = 2√2π, which
0
together with Eq. (8) leads to the relaxation time
1 ε2
τ−1(ε)= − . (31)
0 6
For example, for ε = 0.4 we obtain the theoretical relaxation time, τ 7.1. As
0
≈
showninFig.6forε=0.4inthesimulationthequantity(1 T /T )decreases
DSMC fit
−
by the factor of 10 in the time span ∆τ =25, ranging from τ =10 to τ =35. This
gives the estimate for the relaxation time τ = 25/log(10) 10.8, in agreement
0
≈
with the above theoretical value τ =7.1.
0
From the above theory we expect that the relaxation time increases with ε.
While thistendencyisconfirmedforsmallcoefficientsofrestitutionfromε=0.1to
ε=0.4,Fig.6,therelaxationtime seemsto saturateforlargerε 0.6,Fig.7.This
≥
is, presumably, a finite size effect: For large values of ε the number of particles in
ε=0.6
0
Tfit ε=0.9
/ C
M
S
TD-2 ε=0.6
-
1
ε=0.7
ε=0.8
-4
0 10 20
τ
Fig.7. SameasFig.6butforε=0.6,ε=0.7,ε=0.8,andε=0.9.
February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068˙poeschel
10 T. P¨oschel, N.V. Brilliantov, and A. Formella
the tail, which is determined by the threshold velocity c∗, increasing with ε, is not
large enough to develop a significant tail. Therefore for such systems the apparent
relaxation occurs much faster then it would be in a sufficiently large system. The
relaxation time is mainly determined by the number of particles in the tail, rather
then bythe coefficients ofrestitutionε.Nevertheless,evenforthese systems,which
are not sufficiently large for a quantitative study of the tail relaxation, the latter
occurs anomalously slow.
5. Conclusion
We studied analytically and numerically the impact of the high-energy tails of the
velocitydistributionfunctioningranulargasesinthehomogeneouscoolingstateon
the cooling coefficient and on the relaxationtime towards the steady state velocity
distribution.
InouranalyticaltheoryweusedanuniversalAnsatzforthevelocitydistribution
functionfortheentirerangeofvelocities.ThisAnsatzcomprisesthemainpartofthe
distribution function, where the velocities are of the order of the thermal velocity,
v v as well as the tail, where v v .
T T
∼ ≫
We derived the coefficients of the proposed Ansatz, which allows to estimate
the cooling coefficient and to characterizethe impact of the high-energytail on the
cooling. Our analytical results are in qualitative agreement with numerical data
obtained by Direct Simulation Monte Carlo of 108 particles.
In the simulations we found anomalously slow relaxation of the velocity distri-
bution function to its stationary form as compared to the relaxation of molecular
gases, where the relaxation occurs during few collision times. Instead, for granular
gases, we observed much longer relaxation, in the order of 20-30 collisions per par-
ticle.Toexplainthisslowrelaxationwedevelopedatheory,whichpredictsthatthe
relaxationtimeincreaseswithdecreasinginelasticity,inagreementwiththenumer-
ical observations.For a large coefficient of restitution the relaxationtime saturates
with increasing ε, which may be attributed to finite size effects, that is, the system
of 108 particles is not large enough for the quantitative numerical analysis of the
tail relaxation.
References
1. R. Brito, M. H. Ernst, Extension of Haff’s cooling law in granular flows, Europhys.
Lett. 43 (1998) 497.
2. I. Goldhirsch, G. Zanetti, Clustering instability in dissipative gases, Phys. Rev.Lett.
70 (1993) 1619.
3. J.M.Montanero,A.Santos,Computersimulationofuniformlyheatedgranularfluids,
Granular Matter 2(1999) 53–64.
4. S.E. Esipov,T. P¨oschel, Thegranular phasediagram, J. Stat.Phys.86 (1997) 1385.
5. P. K. Haff, Grain flow as a fluid-mechanical phenomenon, J. Fluid Mech. 134 (1983)
401.
6. A.Goldshtein,M.Shapiro,Mechanicsofcollisionalmotionofgranularmaterials.Part
1: General hydrodynamicequations, J. Fluid Mech. 282 (1995) 75.