Table Of ContentCharit´e, Biochemistry 03/01
Hydrodynamics and transport coefficients for Granular Gases
Nikolai Brilliantov1,2, Thorsten Po¨schel1
1Institut fu¨r Biochemie, Charit´e, Monbijoustraße 2, 10117 Berlin, Germany
2Moscow State University, 119899 Moscow, Russia
(Dated: February 23, 2009)
3
0 The hydrodynamics of granular gases of viscoelastic particles, whose collision is described by
0 an impact-velocity dependent coefficient of restitution, is developed using a modified Chapman-
2 Enskog approach. We derive the hydrodynamic equations and the according transport coefficients
n with the assumption that the shape of the velocity distribution function follows adiabatically the
a decayingtemperature. Weshownumericallythatthisapproximationisjustifieduptointermediate
J dissipation. The transport coefficients and the coefficient of cooling are expressed in terms of
0 the elastic and dissipative parameters of the particle material and by the gas parameters. The
1 dependence of these coefficients on temperature differs qualitatively from that obtained with the
simplifying assumption of a constant coefficient of restitution which was used in previous studies.
] The approach formulated for gases of viscoelastic particles may be applied also for other impact-
h
velocity dependencies of the restitution coefficient.
c
e
PACSnumbers: 81.05.Rm,47.10.+g,05.20.Dd,51.10.+y
m
-
t
a I. INTRODUCTION been derived [19]. Using this generalized Hertz law the
t coefficient of restitution of viscoelastic particles can be
s
. given as a function of the impact velocity and material
t Granular systems composed of a large number of dis-
a parameters [20]:
sipatively interacting particles behave in many respects
m
as a continuous medium and may be, in principle, de-
3
- scribed by a set of hydrodynamic equations with appro- ε=1 C Aα2/5g1/5+ C2A2α4/5g2/5 ... (1)
d − 1 5 1 ∓
n priate boundary conditions. Although this approach is
o successively used in various fields of engineering and soil with
c mechanics (e.g. [1, 2]) a first-principle theory for dense
[
granular media is still lacking. Hydrodynamics may be g ~e (~v1 ~v2) = ~e ~v12 . (2)
1 also applied to much simpler systems, such as rarefied ≡| · − | | · |
v granulargases. It hasbeen usedto describe many differ- Theunitvector~e=~r /r specifies thecollisiongeome-
12 12
2 ent processes,e.g. rapidgranularflows,structure forma- try,i.e.,therelativeposition~r =~r ~r oftheparticles
12 1 2
5 tion,etc. (see[3]foranoverview). Forthesesystemsthe at the collision instant. Their pre-col−lision velocities are
1
hydrodynamic equations are not postulated but derived given by~v and~v . The elastic constant
1 1 2
from the Boltzmann equation. The correspondingtrans-
0
3 port coefficients are not phenomenologicalconstants, in- 3 3/2 Y√Reff
0 stead they are obtained by regular methods, such as the α= 2 meff(1 ν2), (3)
/ Grad method [4] or the Chapman-Enskog method [5]. (cid:18) (cid:19) −
t
a In most of the studies, which address the derivation of depends on the effective mass and radius meff
m hydrodynamic equations and kinetic coefficients, it was m m /(m +m ),Reff R R /(R +R )ofthecollidin≡g
1 2 1 2 1 2 1 2
- assumedthatthe coefficientofrestitutionε is a material spheres, on the Young≡modulus Y, and on the Poisson
d
constant, e.g. [6, 7, 8, 9, 10, 11, 12, 13].
n ratio ν of the particle material. The dissipative coeffi-
o This assumption simplifies the analysis enormously, cient A is a function of dissipative and elastic constants
c however, it is neither in agreement with experimental (see [19] for details). Finally, C is a numerical constant
1
v: observations (e.g. [14, 15, 16]) nor with basic mechanics [17, 20]:
of particle collisions [17]. The coefficient of restitution
i
X depends on the impact velocity g and tends to unity for √π Γ(3/5)
C = 1.15344. (4)
r very small g. As a consequence, particles behave more 1 21/552/5Γ(21/10) ≈
a
and moreelastically as the averagevelocity of the grains
decreases. The simplest collision model which accounts Equation(1)describespureviscoelasticinteraction. The
for dissipative material deformation, is the model of vis- assumption of viscoelastic deformation is justified if the
coelastic particles. It is assumed in this model that the impact velocity is not too large to avoid plastic defor-
elasticstressinthe bulk ofthe particlematerialdepends mation of the particles and not too small to neglect sur-
linearly on the strain, whereas the dissipative stress de- face effects such as adhesion, van der Waals forces etc.
pends linearly on the strain rate. We alsoassumethatthe rotationaldegreesoffreedomof
Based on the fundamental work by Hertz [18] the in- particlesmaybeneglectedandconsideragranulargasof
teractionforce between colliding viscoelastic spheres has identicalparticlesintheabsenceofexternalforces. Then
2
the coefficient of restitution gives the velocities of parti- II. VELOCITY DISTRIBUTION AND
cles after a collision~v′,~v′ in terms oftheir values before TEMPERATURE IN THE HOMOGENEOUS
1 2
the collision: COOLING STATE
′ 1+ε A. Evolution equations for temperature and for
~v =~v (~e ~v )~e (5)
1/2 1/2∓ 2 · 12 the second Sonine coefficient
The Boltzmann equation for a granular gas of vis-
The impact velocity dependence of the coefficient of
coelasticparticlesinthehomogeneouscoolingstatereads
restitution implies serious consequences for the granu-
[24, 27]
lar gas dynamics: For the simplified case ε = const. the
form of the velocity distribution function f(~v,t) is char-
∂
acterizedby a time-independent scaled function f˜(~c). It f(~v ,t) = σ2g (σ) d~v d~eΘ( ~v ~e) ~v ~e
1 2 2 12 12
∂t − · | · |
depends only on the scaled velocity ~c = ~v/vT, where Z′′ Z ′′
[χf(~v ,t)f(~v ,t) f(~v ,t)f(~v ,t)]
vT(t) = 2T(t)/m is the thermal velocity and T(t) is × 1 2 − 1 2
the granular temperature [21, 22, 23]. The distribution g2(σ)I(f,f) , (6)
p ≡
function depends on time only via the time dependence
of the granular temperature, its shape is time indepen- where σ = 2R is the particle diameter. The velocities
′′ ′′
dent. The small deviations of the velocity distribution ~v1 and ~v2 denote the pre-collision velocities of the in-
function from the Maxwell distribution are determined verse collision, which leads to the after-collision veloc-
by the time-independent coefficient of restitution. For ities ~v1 and ~v2. The factor ~v12 ~e characterizes the
| · |
granular gases of viscoelastic particles, however, the ef- length of the collision cylinder of cross-section σ2 and
fective value of ε changes with time along with the ther- theHeavisidestep-functionΘ( ~v12 ~e)assuresthatonly
− ·
mal velocity v (t), which gives a typical velocity of the approaching particles collide. The contact value of pair
T
particles. Therefore, the shape of the velocity distribu- correlationfunction g2(σ) accountsforthe increasedcol-
tion function evolves in a rather complicated way [24]. lision frequency due to excluded volume effects. Finally,
thefactorχinthegaintermaccountsfortheJacobianof
Thetime dependenceofthe velocitydistributionfunc- ′′ ′′
the transformation (~v ,~v ) (~v ,~v ) and for the ratio
tion has to be taken into account when the hydrody- 1 2 → 1 2 ′′
of the lengths of the collision cylinders ~v ~e / ~v ~e
namic equations and the transport coefficients are de- | 12· | | 12· |
for the direct and the inverse collision. For the case of
rived. Since the simple scaling is violated, the stan-
spherescollidingwithaconstantcoefficientofrestitution
dard methods of kinetic theory of gases, developed for
χ = 1/ε2, while for viscoelastic spheres, with ε = ε(g)
ε=const. (e.g. [10, 12, 25]) must be revised.
given by Eq. (1), it reads [24, 27]
As it follows from our analysis the transport coeffi-
cients for gases of viscoelastic particles depend on tem-
11
peratureandtimeratherdifferentlyascomparedwiththe χ=1+ C1Aα2/5 ~v12 ~e 1/5
5 | · |
case ε = const. Correspondingly, the behavior of gases
66
of viscoelastic particles differs qualitatively from that + C2A2α4/5 ~v ~e 2/5+... . (7)
25 1 | 12· |
of gases of particles with the simplified collision model
ε=const. The dependence of χ on the impact velocity does not
We wish to remark that the impact velocity depen- allow to derive from the Boltzmann equation a time-
dence of ε has been already taken into account in Ref. independentequationforthescaleddistributionfunction
[26] for the hydrodynamic description of granular shear f˜. Contrary to the case of ε = const. [21, 22, 23] the
flow. In this study an empirical expression of the coef- scaled distribution function depends explicitly on time.
ficient of restitution was applied and a Maxwellian ve- Therefore, we write for a gas of viscoelastic particles
locity distribution was assumed. Moreover, the authors
n ~v
have used the standard Chapman-Enskog method with- f(~v,t)= f˜(~c,t) ~c= , (8)
out the modifications required for gases of dissipatively vT3(t) vT(t)
colliding particles. These modifications for the case of
with the number density of the granular gas n and the
ε=const. have been extensively elaborated in [10].
thermal velocity v (t) defined by the granular tempera-
T
Theaimofthepresentstudyistodevelopacontinuum
ture:
descriptionofgranulargasesofviscoelasticparticles. We
derivethehydrodynamicequationsalongwiththetrans- 3 mv2 3 mv2(t)
nT(t)= d~v f(~v,t)= n T . (9)
portcoefficientsandthecoefficientofcooling. Intherest 2 2 2 2
of the paper, in Sec. II-V, we discuss in detail the most Z
important case of the three-dimensional gases, while in Hence, the shape of the velocity distribution function,
Sec. VI we present the results for the two-dimensional characterized by the rescaled function f˜(~c,t), does not
systems,whicharefrequentlyaddressedinmoleculardy- persist but evolves along with temperature [24, 27]. We
namics studies. wish to stress that the time dependence of f˜(~c,t) is
3
caused by the dependence of the factor χ on the im- where
pact velocity. Contrary, for a gas of simplified particles
(ε=const.) we obtain χ=1/ε2 =const. and, therefore, ∆ψ(~c ) ψ(~c′) ψ(~c ) (16)
i ≡ i − i
the rescaled distribution function is time independent,
f˜(~c,t)=f˜(~c). denotes the change of some function ψ(~ci) according to
For slightly dissipative particles the velocity distribu- a collision. The coefficients µ depend on time via the
p
tion function is close to the Maxwell distribution. It time dependence of a2. For small enoughdissipation the
may be described by a Sonine polynomial expansion moments µ and µ may be obtained as expansions in
2 4
′
[8, 22, 23, 28]: the time-dependent dissipative parameter δ ,
∞
f˜(~c,t)=φ(c) 1+ ap(t)Sp c2 , (10) δ′(t) Aα2/5[2T(t)]1/10 δ 2T(t) 1/10 , (17)
! ≡ ≡ T
Xp=1 (cid:0) (cid:1) (cid:20) 0 (cid:21)
where φ(c) π−3/2exp c2 is the scaled Maxwell dis- with δ Aα2/5T1/10 and with the initial temperature
tribution, S≡(x) are the S−onine polynomials, ≡ 0
p T . These expansions read [24, 27]
(cid:0) (cid:1) 0
S (x)=1
0
2 2
S1(x)=−x2+ 23 (11) µ2 = Aknδ′kan2
k=0n=0
XX (18)
x2 5x 15 2 2
S2(x)= 2 − 2 + 8 , etc. µ4 = Bknδ′kan2
and ak(t) are the time-dependent Sonine coefficients, kX=0nX=0
which characterize the form of the velocity distribution
where and arenumericalcoefficients. Theymay
kn kn
[23, 24]. The first Sonine coefficient is trivial, a = 0, A B
1 be writteninthe compactmatrix notation(rowsreferto
due to the definition of temperature [22, 23], while the
the first index):
othercoefficientsquantifydeviationsofthemoments ck
of the velocity distribution from the moments of the
Maxwell distribution ck , e.g. (cid:10) (cid:11) 0 0 0
0
ˆ= ω 6 ω 21 ω
(cid:10) (cid:11)c4 c4 A 0 25 0 2500 0
a2 = (cid:10) (cid:11)hc−4i(cid:10)0 (cid:11)0 . (12) −ω1 −141090ω1 −6440604010ω1 (19)
Thus, the first nontrivialSonine coefficienta character-
2 0 4√2π 1√2π
izes the fourth moment of the distribution function. 8
For small enough inelasticity (ε & 0.6) the distribu- Bˆ= 258ω0 910235ω0 −12556070ω0
tion function is well approximated by the second Sonine
77ω 476973ω 4459833 ω
coefficient a [22, 23, 28], i.e. higher coefficients a = 0 −10 1 − 44000 1 70400000 1
2 k
for k 3 may be neglected. With this approximation
≥ with
the evolution a granular gas of viscoelastic particles in
the homogeneous cooling state is described by a set of 21
coupled equations for the granular temperature and for ω0 2√2π21/10Γ 6.485
≡ 10 ≈
the second Sonine coefficient [24, 27]: (cid:18) (cid:19) (20)
16
dT = 2BTµ ζT (13) ω1 ≡√2π21/5Γ 5 C12 ≈9.285.
2 (cid:18) (cid:19)
dt −3 ≡−
da2 4 4 The coupled equations (13) and (14) together with (17),
= Bµ (1+a ) Bµ . (14)
2 2 4
dt 3 − 15 (18) and (19) determine the evolution of a granular gas
ofviscoelasticparticlesinthehomogeneouscoolingstate.
The coefficient B B(t) = v (t)σ2g (σ)n is propor-
≡ T 2 In particular, they define the velocity distribution func-
tional to the mean collision frequency. The moments of
tion which is the starting point for the investigation of
the collision integral, µ and µ read with the approxi-
mation f˜=φ(c) 1+a 2(t)S c24 : inhomogeneous gases.
2 2
1 (cid:2) (cid:0) (cid:1)(cid:3)
µ = d~c d~c d~eΘ( ~c ~e) ~c ~e φ(c )
p 1 2 12 12 1 B. Adiabatic approximation for the second Sonine
−2 − · | · |
Z Z Z coefficient
φ(c ) 1+a S c2 +S c2 +a2S c2 S c2
× 2 2 2 1 2 2 2 2 1 2 2
In the limit of small dissipation, δ 1, the coupled
(cid:8) (cid:2) (cid:0) (cid:1) (cid:0) (cid:1)(cid:3) ∆(cp(cid:0)+(cid:1)cp) (cid:0)(1(cid:1)5(cid:9)) ≪
× 1 2 equations (13),(14) may be solved analytically. In linear
4
approximation with respect to δ the solution reads [24, Using (18) for µ , µ we find a as an expansion in the
2 4 2
′
27] small parameter δ :
T(t) t −5/3 a2 = a21δ′+a22δ′2+... (28)
= 1+ (21) 3ω
T0 (cid:18) τ0(cid:19) a21 = −20√02π ≈−0.388
where we introduce the characteristic time 12063 ω2 27 ω
a = 0 + 1 2.752.
22
τ−1 = 16q δτ (0)−1 = 48q δ4σ2n πT0 , (22) 640000 π 40√2π ≈
0 5 0 c 5 0 m Figures 1 and 2 show a in adiabatic approximationdue
r 2
to Eq. (28) together with the numerical solution of Eqs.
withtheinitialmeancollisiontimeτ (0)andtheconstant
c (13,14) and the result of the linear theory, Eqs.(24,25).
Theadiabaticapproximationisratheraccurateforsmall
21 C dissipation, δ < 0.05, after the initial relaxation (5-10
q =21/5Γ 1 0.173. (23)
0 10 8 ≈ collisions per particle) has passed. Even for the larger
(cid:18) (cid:19)
dissipation (δ = 0.2) it is in agreement with the numer-
Correspondingly in linear approximation the second So- ical result. This value of the dissipative parameter cor-
nine coefficient depends on time as [24] responds to the initial coefficient of restitution ε 0.75
≈
for the thermal velocity. The adiabatic approximation
12
a (t)= w(t)−1 Li[w(t)] Li[w(0)] , (24) becomes more and more accurate as the system evolves.
2
− 5 { − }
Contrary,thelineartheoryfailswithincreasingδandbe-
with comes qualitatively incorrectfor largerδ. In Figs. 1 and
2theevolutionofthesystemhasstartedfromaMaxwell
1/6
w(t) exp (q δ)−1 1+ t (25) distribution,i.e.,a2(0)=0. Fig. 3showsthe relaxation
0
≡ " (cid:18) τ0(cid:19) # of a2 fromcertaininitial conditions a2(0)6=0 which cor-
respond to non-Maxwellian distribution functions. The
and with the logarithmic integral relaxationoccurs during the first 5-10 collisions per par-
x 1 ticle. Then the adiabatic approximation becomes valid
Li(x) dt. (26) evenfor intermediate values ofthe dissipative parameter
≡ ln(t)
Z0 δ. The adiabatic approximation yields very accurate re-
For small dissipation δ the time dependence of a re- sultsfortheevolutionoftemperature. Itcoincidesalmost
2
veals two different regimes: (i) fast initial relaxation on perfectlywiththenumericalsolutionfortheshortandfor
the meancollision-time scale τ (0)and(ii) subsequent the long time scales and for all values of the dissipative
c
slowevolutiononthetimesca∼le τ τ (0),i.e. onthe parameter 0 < δ 0.2 (see Fig. 4). Thus we conclude
time scale of the temperature e∼volu0t≫ion.c Therefore the thatthe adiabatic≤approximationmay be appliedfor the
coefficient a (and hence the form of the velocity distri- hydrodynamic description of granular gases. By means
2
butionfunction)evolvesinaccordancewithtemperature. of Eq. (28) for the second Sonine coefficient we can de-
Forthehydrodynamicdescriptionofgranulargaseswe termine the moments µ2, µ4, and the cooling coefficient
assume that there exist well separated time and length ζ:
scales. The short time and length scales are given by
2 2 2T
the mean collisiontime and the mean free path, and the ζ = µ2B = nσ2g2(σ) ω0δ′ ω2δ′2+... ,
3 3 m −
longtime andlengthscalesarecharacterizedbythe evo- r
(cid:0) (cid:1)(29)
lution of the hydrodynamic fields (to be defined in the
where
next section) and their spatial inhomogeneities. The hy-
drodynamic approach corresponds to the coarse-grained ω ω + 9 ω2 2 . (30)
descriptionofthe system,whereallprocesseswhichtake 2 ≡ 1 500 0 π
r
place on the short time and length scales are neglected. Below we will need the derivatives of a and ζ with
2
Therefore the first stage of the relaxation of the velocity respect to temperature and density. Using
distribution function does not affect the hydrodynamic
′ ′
∂δ δ
description and only the second stage of its evolution on T = (31)
∂T 10
the time scale τ τ (0) is to be taken into account.
0 c
≫ ′
To this end we apply an adiabatic approximation: We according to the definition of of δ , Eq. (17), we obtain
omit the term da /dt in the left-hand side of Eq. (14)
2 ∂a 1 1
whichdescribesthefastrelaxationandassumethata2 is T ∂T2 = 10a21δ′+ 5a22δ′2 (32)
determinedbythecurrentvaluesofµ andµ duetothe
2 4 ∂ζ 2 3 7
presenttemperature. Hence,intheadiabaticapproxima- T = B ω δ′ ω δ′2+... (33)
0 2
∂T 3 5 − 10
tion a2 is determined by (cid:18) (cid:19)
∂ζ 1
5µ (1+a ) µ =0. (27) = ζ(0). (34)
2 2 4 ∂n n
−
5
0 0
δ=0.01
δ=0.01
-0,5 -0,5
] ]
0 0
0 0
1 -1 1 -1
x x
[ [
a 2 a 2 δ=0.05
-1,5 δ=0.05 -1,5
-2 -2
0 5 10 15 0 1 2 3 4 5
time [τ ] 10 10 10 tim1e0 [τ ] 10 10
c
c
0 0
δ=0.1 δ=0.1
-1 -1
] ]
0 0
0 0
x1 -2 x1 -2
[ [
2 2
a a
-3 -3
0 5 10 15 0 1 2 3 4 5
time [τ ] 10 10 10 tim1e0 [τ ] 10 10
c
c
4 4
2 2
δ=0.2
0] 0 0] 0 δ=0.2
0 0
1 1
a [x2-2 a [x2-2
-4 -4
-6 -6
0 5 10 15 0 1 2 3 4 5
time [τc] 10 10 10 tim1e0 [τ ] 10 10
c
FIG. 1: Evolution of the second Sonine coefficient a in the
2 FIG. 2: The same as in Fig. 1, but for longer time. The
homogeneouscooling stateontheshorttimescale. Solidline
adiabatic approximation improvesas thesystem evolves.
– numerical solution of Eqs.(13,14), dashed line – adiabatic
approximation (28),dot-dashedline–theresult ofthelinear
theory,Eqs.(24,25). Thetimeisgivenincollision unitsτc(0).
III. HYDRODYNAMIC EQUATIONS AND
TRANSPORT COEFFICIENTS
A. Derivation of hydrodynamic equations
Forthederivationofhydrodynamicequationsfromthe
Boltzmann equation it is assumed that there exist well
separated time and length scales. As already briefly dis-
6
1
4 δ=0.01 0,8 δ=0.01
2
]
0 0,6
0 0
1 0 T
[x2 T / 0,4 δ=0.2
a
-2
0,2
-4
0
0 2 4 6 8 10 12 0 5 10 15 20
time [τ ] time [τ ]
c c
δ=0.2
4 δ=0.2
-2
10
2
]
[x100 0 T / T010-4 ~t-5/3
a 2-2 10-6
-4 -8
10
0 1 2 3 4 5
0 2 4 6 8 10 12 10 10 10 10 10 10
time [τ ] time [τ ]
c c
FIG.3: Thesame asin Fig. 1,butfor differentinitial values FIG.4: Evolutionoftemperatureinthehomogeneouscooling
of a . The fast relaxation takes place during the first 5-10 state on the short time scale (top) and for longer time (bot-
2
collisions, independently on a (0). tom). The notations are the same as in Fig. 1. The dotted
2
line shows the asymptotic power-law at t→∞.
cussed in the preceding Section, in granular gases there
are (at least) two sets of scales: The microscopic scales and does not depend on time. Hence the evolution of
are characterized by the mean collision time and the temperature describes the evolution of the velocity dis-
mean free path. The macroscopic scales are given by tribution function exhaustively. For granular gases of
the characteristic time of the evolution of the hydrody- viscoelastic particles the shape of the velocity distribu-
namicfieldsandthesizeoftheirspatialinhomogeneities. tionfunction,characterizedbytheSoninecoefficientsa ,
k
Hence,scaleseparationmeansthatthemacroscopicfields is time dependent due to the time dependence of these
varyveryslowlyinspaceandtimeifmeasuredinthemi- coefficients. In the previous section we exemplified this
croscopic units. This assumption allows for a gradient property for the firstnon-trivialcoefficient a (t), assum-
2
expansion in space and time, i.e. the application of the ing that for small dissipation higher coefficients may be
Chapman-Enskog method [5]. disregarded. Therefore, the evolution of granular gases
The application of the Chapman-Enskog approach to of viscoelastic particles is described by the time depen-
granular gases is more sophisticated than its application dence of temperature, i.e., by the second moment of the
to moleculargases: Formoleculargasesthe unperturbed velocity distribution function, and by the time depen-
solution (the basic solution) for the velocity distribution dence of its higher-order even moments which are char-
function is the Maxwell distribution while for granular acterized by the Sonine coefficients a , a , etc. [see Eq.
2 3
gases the basic solution is the time-dependent velocity (12)]. Hence, the hydrodynamic description of granu-
distribution of the homogeneous cooling state. lar gases whose particles collide with ε = ε(g) requires
Forgasesofparticlesinthehomogeneouscoolingstate, an extended set of hydrodynamic fields which includes
which collide with ε = const., the velocity distribu- the higher-ordermoments. A regular approachto derive
tion function depends on time only through the time- hydrodynamic equations for this extended set of fields
dependent characteristic velocity. The shape of the ve- is the Grad method [4]. Alternatively, for small dissi-
locity distribution function, given by f˜, is fixed by ε pation an adiabatic approximation can be applied: It
7
is assumed that the shape of the velocity distribution Thestructureofthe hydrodynamicequations,exceptfor
function, albeitvaryingintime,followsadiabaticallythe the cooling term ζT, coincides with those for molecular
current temperature. Correspondingly, the higher-order gases.
moments are determined by the temperature too. With
this approximation a closed set of hydrodynamic equa-
tionsfordensityn(~r,t),velocity~u(~r,t),andtemperature B. Chapman-Enskog approach
T(~r,t) may be derived. These fields are defined respec-
tively by the zeroth, first, and second moments of the Thesystem(37,38,39)isclosedbyexpressingthepres-
velocity distribution function: sure tensor and the heat flux in terms of the hydrody-
namic fields and the fields gradients. To this end we
n(~r,t)= d~vf(~r,~v,t) apply the Chapman-Enskog approach [5]. This method
Z is based on two important assumptions: The evolution
n(~r,t) ~u(~r,t)= d~v~vf(~r,~v,t) , (35) of the distribution function is completely determined by
the evolutionof its firstfew moments, i.e., it depends on
Z
3 1 space and time only through the hydrodynamic fields:
n(~r,t) T (~r,t)= d~v mV2f(~r,~v,t) ,
2 2
Z f(~r,~v,t)=f[~v,n(~r,t),~u(~r,t),T (~r,t)] . (45)
where V~ ~v ~u(~r,t). For small dissipation the veloc-
≡ −
ity distribution is well approximated by the second So- As the second precondition it is assumed that the gas
nine coefficients a (t), i.e., higher-order coefficients are isonlyslightlyinhomogeneousonthemicroscopiclength
2
neglected. In adiabatic approximation a is determined scalewhichallowsforagradientexpansionofthevelocity
2
by temperature according to (28). distribution function:
Forsimplicityofthenotationweconsideradilutegran-
ular gas and approximate the Enskog factor g (σ) 1. f =f(0)+λf(1)+λ2f(2)+... , (46)
2
≈
Multiplying the Boltzmann equation for an inhomoge-
where each power k of the formal parameter λ corre-
neous gas
sponds to the order k in the spatial gradient. Thus, f(0)
∂ +~v ~ f(~r,~v ,t)=I(f,f) , (36) referstothehomogeneouscoolingstate,f(1) corresponds
∂t 1·∇ 1 tothelinearapproximationwithrespecttothefieldsgra-
(cid:18) (cid:19)
dients,f(2)isthesolutionwithrespecttoquadraticterms
correspondingly by v0, ~v , mv2/2, and integrating over
1 1 1 in the field gradients, etc. With these assumptions the
d~v we obtain the hydrodynamic equations (see e.g. [5])
1 Boltzmann equation may be solved iteratively for each
∂n orderinλ,togetherwiththehydrodynamicequationsfor
+ ~ (n~u)=0, (37)
the moments of the velocity distribution function. The
∂t ∇·
∂~u solution in zeroth order in λ yields the velocity distri-
+~u ~~u+(nm)−1~ Pˆ =0, (38) bution function for the homogeneous cooling state f(0)
∂t ·∇ ∇·
and the corresponding evolution of temperature. The
∂T 2
+~u ~T + Pˆ : ~~u+ ~ ~q +ζT =0. (39) function f(0) is then used to compute Pˆ and ~q, yielding
∂t ·∇ 3n ∇ ∇· Pˆ(0) = δ p, ~q(0) = 0 and the hydrodynamic equations
(cid:16) (cid:17) ij
The cooling coefficient ζ in the sink term ζT may be for the ideal fluid. These first-order equations contain
written as only linear gradient terms. Then f(1) may be found em-
ploying the first-order hydrodynamic equations and the
σ2m
ζ(~r,t)= d~v d~v d~eΘ( ~v ~e) ~v ~e distribution function f(0). The obtained f(1) as well as
1 2 12 12
12nT Z Z Z − · | · | thecorrespondingexpressionsforPˆ(1) and~q(1) arelinear
f(~r,~v1,t)f(~r,~v2,t)(~v12 ~e)2 1 ε2 . (40) in the field gradients:
× · −
The pressuretensorPˆ andthe heatflux(cid:0)~q arede(cid:1)fined by 2
P =pδ η u + u δ ~ ~u
ij ij i j j i ij
− ∇ ∇ − 3 ∇· (47)
P (~r,t) = D V~ f(~r,~v,t)d~v+pδ (41) (cid:18) (cid:19)
ij ij ij ~q = κ~T µ~n.
Z (cid:16) (cid:17) − ∇ − ∇
~q(~r,t) = S~ V~ f(~r,~v,t)d~v, (42) The transport coefficients η, κ and µ in these equa-
Z (cid:16) (cid:17) tions are expressed in terms of f(1). Hence, within the
where p = nT is the hydrostatic pressure. The velocity Chapman-Enskog approach for each order in the gradi-
tensor Dij and the vector S~ read ent expansion a closed set of equations may be derived.
Keeping only first-order field gradients for the distribu-
1
D (V~) m V V δ V2 (43) tionfunctionandrespectivelyforthepressuretensorand
ij i j ij
≡ − 3
(cid:18) (cid:19) fortheheatfluxasinEq. (47),theNavier-Stokeshydro-
mV2 5 dynamicsisobtained. Keepingnext-ordergradientterms
S~ V~ T V~ . (44)
≡ 2 − 2 correspondstotheBurnettorsuper-Burnettdescription.
(cid:16) (cid:17) (cid:18) (cid:19)
8
We will restrict ourselves to the Navier-Stokes level [32] The term ζ(1) is found by substituting f = f(0)+λf(1)
and skip for simplicity of the notation the superscript into Eq. (40) and collecting terms of the order (λ).
“(1)” for Pˆ and ~q in the above equations. Thesetermscontainthefactorf(0)(~v,~r,t)f(1)(~v,~rO,t)in
The Chapman-Enskog scheme also assumes a hierar- the integrand. They vanish upon integration according
chy of time scales and respectively a hierarchy of time to the different symmetry of the functions f(0) and f(1),
derivatives: thus, ζ(1) =0 [33].
Since the distribution function f(0) is known, we can
∂ ∂(0) ∂(1) ∂(2)
= +λ +λ2 +... , (48) evaluate those terms in (53) which depend only on f(0):
∂t ∂t ∂t ∂t
whereeachorderk in the time derivative,∂(k)/∂t≡∂t(k) ∂(1) +~v1 ~ f(0) = ∂f(0) ∂(1)n +~v1 ~n
corresponds to the related order in the space gradient. ∂t ·∇ ∂n ∂t ·∇
(cid:18) (cid:19) (cid:18) (cid:19)
Consequently,thehighertheorderinthespacegradient, ∂f(0) ∂(1)~u
the sloweris the accordingtime variation. Usingthe for- + +~v1 ~~u
∂~u · ∂t ·∇
mal expansion parameter λ we can write the Boltzmann (cid:18) (cid:19)
equation, + ∂f(0) ∂(1)T +~v ~T . (56)
1
∂T ∂t ·∇
(cid:18) (cid:19)
∂(0) ∂(1)
+λ +...λ~v ~ f(0)+λf(1)+...
∂t ∂t 1·∇ With~v1 =V~ +~u,thetimederivativesofn,~uandT given
(cid:18) (cid:19)(cid:16) (cid:17) in (55), and the relations
=I f(0)+λf(1)+... , f(0)+λf(1)+... (49)
h(cid:16) (cid:17) (cid:16) (cid:17)i ∂f(0) 1 ∂f(0) ∂f(0)
and collect terms of the same order in λ. The equation = f(0), = , (57)
∂n n ∂~u − ∂V~
in zeroth order
∂(0) which follow from Eq. (51), we recast Eq. (53) for f(1)
f(0) =I f(0),f(0) (50)
into the form:
∂t
(cid:16) (cid:17)
coincideswithEq. (6)forthehomogeneouscoolingstate. ∂(0)f(1)
AccordingtotheChapman-Enskogschemeweobtainthe +J(1) f(0),f(1) =f(0) ~ ~u V~ ~ logn
∂t ∇· − ·∇
velocity distribution function in zeroth order:
∂(cid:16)f(0) 2 (cid:17) (cid:16) (cid:17)
+ T~ ~u V~ ~T
f(0)(~v,~r,t)= vnT3((~r~r,,tt))(cid:2)1+a2S2(cid:0)c2(cid:1)(cid:3)φ(c), (51) +∂∂Tf(0(cid:18)) 3 V~∇·~ −ui ·∇1(cid:19) ip , (58)
where ~c = (~v−~u)/vT = V~/vT. The corresponding hy- ∂Vi (cid:18)(cid:16) ·∇(cid:17) − nm∇ (cid:19)
drodynamic equations in this order read
where we take into account ζ(1) =0.
∂(0) ∂(0) ∂(0) The right-hand side is known since it contains only
n=0, ~u=0, T = ζ(0)T , (52) f(0). It is convenient, however to rewrite it in a form
∂t ∂t ∂t −
whichshowsexplicitlythedependences onthe fieldsgra-
where ζ(0) is to be calculated using Eq. (40) with f = dients. Employing
f(0). In this way we reproduce the previous result (29).
Collecting terms of the order (λ) we obtain, 1 T T
O ip= ilogT + ilogn, (59)
nm∇ m∇ m∇
∂(0)f(1) ∂(1)
+ +~v ~ f(0)+J(1) f(0),f(1) =0,
∂t ∂t 1·∇ the right-hand side of (58) yields
(cid:18) (cid:19) (cid:16) (cid:17) (53)
where we introduce ∂(0)f(1)
+J(1) f(0),f(1)
∂t
−J(1) f(0),f(1) ≡I f(0),f(1) +I f(1),f(0) . (54) =A~ ~(cid:16)logT +B(cid:17)~ ~ logn+C u , (60)
ij j i
(cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17) ·∇ ·∇ ∇
The corresponding first-order hydrodynamic equations
with
read
∂(1)n = ~ (n~u) A~(V~) = V~T∂f(0) T ∂f(0) (61)
∂t −∇ − ∂T − m ∂V~
∂(∂1t)~u =−~u·∇~~u− n1m∇~p (55) B~(V~) = −V~f(0)− mT ∂∂V~ f(0) (62)
∂(1)T = ~u ~T 2T ~ ~u ζ(1)T . Cij(V~) = ∂ Vjf(0) + 2Tδij∂f(0) . (63)
∂t − ·∇ − 3 ∇· − ∂Vi 3 ∂T
(cid:16) (cid:17)
9
We will calculate these terms below. From the form of where J(1) is defined by Eq. (54).
the right-handside of Eq. (60) we expect the form of its
solution
f(1) =α~ ~ logT +β~ ~ logn+γ u , (64) C. Kinetic coefficients in terms of the velocity
·∇ ·∇ ij∇j i distribution function
whichisthemostgeneralformofascalarfunction,which
depends linearly on the vectorial gradients ~T, ~n and From the definition of the pressure tensor Eq. (41)
∇ ∇
onthetensorialgradients u . Thecoefficientsα~,β~ and anditsexpressionintermsofthefieldgradientsEq. (47)
j i
∇
γ are functions of V~ and of the hydrodynamic fields n, follows
ij
~u and T.
We derive now equations for the coefficients α~, β~ and D f(0)+α~ ~ logT +β~ ~ logn+γ u dV~
ij kl l k
γ by substituting f(1) as given by Eq. (64) into the ·∇ ·∇ ∇
ij Z (cid:16) (cid:17)
first-order equation (60) and equating the coefficients of 2
= η u + u δ ~ ~u , (72)
thecorrespondinggradients. Tothisendweneed∂(0)f(1) − ∇i j ∇j i− 3 ij∇·
t (cid:18) (cid:19)
and, therefore, the time derivatives of the coefficients α~,
β~, and γij: where the tensor Dij V~ has been defined above. The
∂(0)α~ ∂α~ ∂(0)T ∂α~ ∂(0)n ∂α~ ∂(0)T integrals (cid:16) (cid:17)
= + +
∂t ∂T ∂t ∂n ∂t ∂u ∂t
i D f(0)dV~ = D α~dV~ = D β~dV~ =0 (73)
∂α~ ij ij ij
= ζ(0)T , (65) Z Z Z
− ∂T
vanish since D is a traceless tensor and f(0) depends
where we use Eqs. (52) in zeroth order. Similarly we ij
obtain isotropically on V~. Moreover, as will be shown below,
the vectors α~ and β~ are directed along V~, hence, the
∂(0)β~ ∂β~ ∂(0)γ ∂γ
= ζ(0)T , ij = ζ(0)T ij (66) respective integrands are odd functions of V~. Therefore,
∂t − ∂T ∂t − ∂T
only the term with the factor γ u in the left-hand
kl l k
and, respectively, the time-derivatives of the gradients side of Eq. (72) is non-trivial. Eq∇uating the coefficients
∂(0)~ logn=0 of the gradient factor ∇luk we obtain
∂t ∇
2
∂(0) DijγkldV~ = η δliδkj +δljδki δijδkl . (74)
u =0 − − 3
∂t ∇j i Z (cid:18) (cid:19)
(67)
∂(0) For k =j, l=i the last equation turns into
~ logT = ~ζ(0)
∂t ∇ −∇
2
∂ζ(0) ∂ζ(0) D γ dV~ = η δ δ +δ δ δ δ = 10η,
= ~n ~T . ij ji − ii jj ij ij − 3 ij ij −
− ∂n ∇ − ∂T ∇ Z (cid:18) (cid:19)
(cid:18) (cid:19) (cid:18) (cid:19) (75)
The derivatives of ζ(0) are given by Eqs. (33,34). From (δiiδjj = 9 and δijδij = 3 according to the summation
(65), (66), and (67) we obtain: convention) and yields the coefficients of viscosity
∂(0)f(1) = T∂ζ(0)α~ ~ logT η = 1 Dij V~ γji V~ dV~ . (76)
∂t − ∂T ·∇ −10
(cid:18) (cid:19) Z (cid:16) (cid:17) (cid:16) (cid:17)
ζ(0)T ∂β~ +ζ(0)α~ ~ logn ζ(0)T∂γij u . UsingEq. (41)fortheheatfluxanditscorrespondingex-
− ∂T !·∇ − ∂T ∇j i pressionin terms of the field gradients,Eq. (47), we can
performacompletelyanalogouscalculationandarriveat
(68)
the kinetic coefficients κ and µ:
where we use Eq.(34). If we insert ∂(0)f(1) into Eq. (60)
t 1
and equate the coefficients of the gradients we arrive at κ = dV~S~ V~ α~ V~ (77)
a set of equations for the coefficients α~, β~ and γ : −3T ·
ij Z (cid:16) (cid:17) (cid:16) (cid:17)
1
∂ζ(0) µ = dV~S~ V~ β~ V~ , (78)
T α~ +J(1) f(0),α~ =A~ (69) −3n ·
− ∂T Z (cid:16) (cid:17) (cid:16) (cid:17)
∂β~ (cid:16) (cid:17) with S~ defined by Eq. (44). The coefficient of thermal
ζ(0)T ζ(0)α~ +J(1) f(0),β~ =B~ (70)
− ∂T − conductivityκhasthestandardinterpretation,whilethe
∂γ (cid:16) (cid:17) othercoefficientµ doesnothaveananalogformolecular
ζ(0)T ij +J(1) f(0),γij =Cij (71) gases.
− ∂T
(cid:16) (cid:17)
10
IV. COEFFICIENT OF VISCOSITY where we take into account that
2
Theviscositycoefficientisrelatedtothe coefficientγ D D = m2V4 (87)
ij ij ji
3
[see Eq. (76)] which, in its turn, is the solution of Eq.
(71) with the coefficient C in the right-hand side. Let
ij according the definition of D , Eq. (43), with the sum-
ij
us first find an explicit expression for C .
ij mation convention. Moreover, we have used the fourth
According to Eq. (51) the velocity distribution func-
moment of the Maxwell distribution. From Eq. (86) fol-
tion depends ontemperature throughthe thermalveloc-
lows
ityv andadditionallythroughthe secondSoninecoeffi-
T
cient a . Hence the temperature derivative of f(0) reads η
2 γ = . (88)
0
−nT
∂∂fT(0) = 21T ∂∂V~ ·V~f(0)+fM(V)S2 c2 ∂∂aT2 , (79) anMduulstiinpglyEinqg. (E7q6.) (y7ie1l)dsby Dij(V~1), integrating over V~1
(cid:0) (cid:1)
with f (V) being the Maxwell distribution
M
∂η
n 1 10ζ(0)T = dV~1Dij V~1 Cij V~1
f = φ(c) with φ(c)= exp c2 . (80) − ∂T −
M v3 π3/2 − Z (cid:16) (cid:17) (cid:16) (cid:17)
T (cid:0) (cid:1) + dV~1Dij V~1 J(1) f(0),γij . (89)
Using the relation
Z (cid:16) (cid:17) (cid:16) (cid:17)
∂f(0) V ∂f(0) To evaluate the first term in the right-hand side we use
i
= (81)
Eq. (63) for C , the relation
∂V V ∂V ij
i
and Eq. (79) the coefficient Cij reads ∂ D = ∂ m V V 1δ V2 =mV 1+ 1δ ,
ij i j ij j ij
∂V ∂V − 3 3
C V~ = V V 1δ V2 1 ∂f(0) i i (cid:18) (cid:19) (cid:18) ((cid:19)90)
ij i j − 3 ij V ∂V the definition of temperature, Eq. (35), and notice that
(cid:16) (cid:17) (cid:18) + 2δijS(cid:19)2 c2 fM(V)T∂a2 . (82) DInitjeδgirja=tio0nsbinycpeaDrtisjtihseantyriaecldelsess tensor [see Eq. (43)].
3 ∂T
(cid:0) (cid:1)
With
dV~ D C = (91)
1 ij ij
1 ∂f(0) m m 5 Z
= 1+a2S2 c2 fM+ c2 a2fM ∂ 2T ∂f(0)
V ∂V −T T − 2 = dV~ D V f(0)+ dV~ D δ
(cid:18) (cid:19) 1 ij 1j 1 ij ij
(cid:2) (cid:0) (cid:1)(cid:3) (83) Z ∂V1i 3 Z ∂T
we obtain 1 10
= dV~ f(0)mV V 1+ δ = dV~ f(0)mV2
1 1j 1j 3 ij 3 1 1
1 5 Z (cid:18) (cid:19) Z
Cij = Dij 1+a2 S2 c2 + c2 fM(V) =10nT.
−T 2 −
(cid:20) (cid:18) (cid:19)(cid:21)
2 (cid:0) (cid:1) ∂a
+ δ S c2 f (V)T 2 , (84) Forthesecondtermintheright-handsideof (89)weuse
ij 2 M
3 ∂T the definition of J(1), (54), and obtain
(cid:0) (cid:1)
where D V~ has been defined in Eq. (43).
ij
The exp(cid:16)ress(cid:17)ion for Cij determines the right-hand side dV~1DijJ(1) f(0),γij
of Eq. (71) for γij and hence it suggests the form for Z (cid:16) (cid:17)
γij. For small dissipation (when a2 is small) and for = dV~1DijI f(0),γij dV~1DijI γij,f(0) .
− −
smallfieldsgradientswekeeponlytheleadingtermswith Z (cid:16) (cid:17) Z (cid:16) (cid:17)(92)
respect to these variables. Therefore, we seek for γ in
ij
the form
We apply the property of the collision integral [5]
γ
γ V~ = 0D V~ f (V), (85)
ij ij M
T
(cid:16) (cid:17) (cid:16) (cid:17) dV~ D I f(0),γ = dV~ D I γ ,f(0)
where γ is a velocity-independent coefficient, i.e., we 1 ij ij 1 ij ij
0
neglect the dependence of γ on a [34]. The viscosity Z (cid:16) (cid:17) Z (cid:16) (cid:17)
ij 2 σ2
coefficient Eq. (76) reads then = dV~ dV~ f(0) V~ γ V~ d~eΘ V~ ~e
1 2 1 ij 2 12
2 − ·
Z Z (cid:16) (cid:17) (cid:16) (cid:17)Z (cid:16) (cid:17)
η =−1γ00 T1 dV~DijDijfM =−γ0nT , (86) × V~12·~e ∆ Dij V~1 +Dij V~2 , (93)
Z (cid:12) (cid:12) h (cid:16) (cid:17) (cid:16) (cid:17)i
(cid:12) (cid:12)
(cid:12) (cid:12)
