7 Kinetic Theory for Binary Granular Mixtures
at Low-Density
Vicente Garzó
Departamento de F́ısica, Universidad de Extremadura, E-06071 Badajoz, Spain
vicenteg@unex.es
Many features of granular media can be modelled as a fluid of hard spheres
with inelastic collisions. Under rapid flow conditions, the macroscopic be-
havior of grains can be described through hydrodynamic equations. At low-
density, a fundamental basis for the derivation of the hydrodynamic equations
and explicit expressions for the transport coefficients appearing in them is pro-
vided by the Boltzmann kinetic theory conveniently modified to account for
inelastic binary collisions. The goal of this chapter is to give an overview of the
recent advances made for binary granular gases by using kinetic theory tools.
Some of the results presented here cover aspects such as transport properties,
energy nonequipartition, instabilities, segregation or mixing, non-Newtonian
behavior, . . . . In addition, comparison of the analytical results with those ob-
tained from Monte Carlo and molecular dynamics simulations is also carried
out, showing the reliability of kinetic theory to describe granular flows even
for strong dissipation.
1 Introduction
Granular systems have attracted the attention of the physics community in
the past few years, in part because the behavior of these systems under many
conditions exhibit a great similarity to ordinary fluids [1]. These conditions
include rapid, dilute flows where the dominant transfer of momentum and
energy is through binary collisions of the grains. The main difference from
ordinary fluids is the absence of energy conservation, leading to both obvious
and subtle modifications of the usual macroscopic balance equations as well as
the constitutive equations for the irreversible fluxes. However, in spite of the
utility of hydrodynamics to describe rapid granular flows, there are still some
open questions about its domain of validity and the associated constitutive
equations appearing in the hydrodynamic equations [2].
To isolate the effects of such collisional dissipation from other important
properties of granular media, an idealized microscopic model system is usually
http://arxiv.org/abs/0704.1211v1
2 Vicente Garzó
considered: a system composed by smooth hard spheres with inelastic colli-
sions. As in the elastic case, the collisions are specified in terms of the change
in relative velocity at contact but with a decrease in the magnitude of the nor-
mal component measured by a positive coefficient of restitution α ≤ 1. This
parameter distinguishes the granular fluid (α < 1) from the ordinary fluid
(α = 1). Given that the hard sphere system with elastic collisions has been
widely studied for both equilibrium and non-equilibrium statistical mechanics
[3], it is tempting to apply the same methods for the case of inelastic collisions.
However, some care is warranted in translating properties of ordinary fluids to
granular fluids. In this presentation, a kinetic theory description based on the
Boltzmann kinetic equation (which applies at sufficiently low density) will be
considered as the appropriate tool to study granular flows from a microscopic
point of view.
Although many efforts have been devoted in the past few years to the
understanding of granular fluids, the derivation of the form of the constitu-
tive equations with explicit expressions for the transport coefficients is still a
subject of interest and controversy. The conditions to obtain a hydrodynamic
description are expected to be similar to those for normal fluids. For a given
initial state there are two stages of evolution. First, during the kinetic stage
there is rapid velocity relaxation to a “universal” velocity distribution that
depends on the average local density, temperature, and flow velocity. Subse-
quently, the hydrodynamic stage is described through a slower evolution of
these local hydrodynamic fields as they approach uniformity. The solution to
the Boltzmann equation in this second stage is said to be normal, where all
space and time dependence of the distribution function occurs through the
macroscopic hydrodynamic fields. The Chapman–Enskog method [4] provides
a constructive means to obtain an approximation to such a solution for states
whose spatial gradients are not too large. In this case, the explicit form of this
normal solution is given as a perturbation expansion in the spatial gradients of
the fields. This solution is then used to evaluate the fluxes in the macroscopic
balance equations in terms of these gradients. To lowest order the balance
equations become the granular Euler equations while to second order they are
the granular Navier–Stokes equations. In carrying out this analysis, explicit
forms for the transport coefficients are obtained as functions of the coefficient
of restitution and other parameters of the collision operator. In this general
context, the study of hydrodynamics for granular gases is the same as that
for ordinary fluids.
The derivation of hydrodynamics from the inelastic Boltzmann equation
has been widely covered in the case of a monodisperse gas where the particles
are of the same mass and size. As for elastic collisions, the transport coeffi-
cients are given in terms of the solutions of linear integral equations [5, 6],
which are approximately solved by using Sonine polynomial expansions. The
estimates for the transport coefficients provided by the Sonine solution com-
pare in general quite well with both direct Monte Carlo simulation (DSMC)
of the Boltzmann equation and molecular dynamics (MD) simulation of the
Kinetic Theory for Binary Granular Mixtures at Low-Density 3
gas, even for relatively strong degrees of dissipation [7, 8, 9, 10]. This good
agreement supports the formal theoretical analysis and the claim that hydro-
dynamics is not limited to nearly elastic particles [11, 12, 13].
Nevertheless, a real granular system is generally characterized by some
degrees of polydispersity in density and size, which leads to phenomena very
often observed in nature and experiments, such as separation or segregation.
Needless to say, the study of granular mixtures is much more complicated
than for a monodisperse gas since not only the number of transport coeffi-
cients in a multicomponent system is higher than that of a single gas, but
also they depend on parameters such as masses, sizes, composition as well
as several independent coefficients of restitution αij . Due to these difficulties,
studies for multicomponent gases are more scarce in the literature. Many of
the previous attempts [14] to derive hydrodynamics from kinetic theory were
carried out in the quasi-elastic limit where the equipartition of energy can be
considered as an acceptable assumption. In addition, according to this level of
approximation, the inelasticity is only accounted for by the presence of a sink
term in the energy balance equation, so that the expressions for the transport
coefficients are the same as those obtained for ordinary fluids. However, the
theoretical prediction of the failure of energy equipartition in multicomponent
granular gases [15] has been confirmed by computer simulations [16], and even
observed in real experiments [17].
Although the possibility of nonequipartition was already pointed out many
years ago [18], it has not been until recently that a systematic study of the
effect of nonequipartition on the Navier–Stokes hydrodynamic equations has
been carefuly carried out [19, 20]. These new equations and associated trans-
port coefficients provide a somewhat more stringent test of the analysis since
the parameter space is much larger. As in the monodisperse case, explicit ex-
pressions for the transport coefficients requires also to consider Sonine poly-
nomial expansions. The numerical accuracy of this Sonine expansion has been
confirmed by comparison with Monte Carlo simulations of the Boltzmann
equation in the cases of the shear viscosity [21] and the tracer diffusion [22]
coefficients. Exceptions to this agreement are extreme mass or size ratios and
strong dissipation, although these discrepancies between theory and simula-
tion diminish as one considers more terms in the Sonine polynomial approxi-
mation [22].
The explicit knowledge of the Navier-Stokes transport coefficients allows
quantitative application of the nonlinear hydrodynamic equations to a num-
ber of interesting problems for granular mixtures, such as to quantify the
violation of the Einstein relation [23, 24] or the Onsager reciprocal relations
[20], the stability analysis of the homogeneous cooling state [20], and segrega-
tion induced by a thermal gradient [25]. In all the cases, the analysis clearly
shows the important role played by the inelasticity in the different physical
situations.
The analogy between rapid granular flow and ordinary fluids can be also
extended to many other transport situations. A particularly simple case, al-
4 Vicente Garzó
lowing detailed analysis even in far from equilibrium conditions is the simple
or uniform shear flow (USF) problem. Macroscopically, it is characterized by
uniform density and temperature and a constant mean velocity profile. This
is a well-known non-equilibrium problem widely studied, for both granular
monodisperse [13, 26, 27] and ordinary gases [28]. Nevertheless, the nature of
this state is quite different in each system. While for elastic fluids the temper-
ature increases monotonically in time due to viscous heating, a steady state
is possible for granular media when the effect of the viscous heating is exactly
compensated by the dissipation in collisions. Thus, in the steady state, there
is an intrinsic connection between the shear field and dissipation so that the
collisional cooling sets the strength of the velocity gradient. As a consequence,
the USF state is inherently non-Newtonian and the rheological properties of
the system cannot be obtained from the Navier-Stokes description, at least
for finite dissipation [29].
The aim of this chapter is to offer a short review of recent results ob-
tained for binary granular mixtures from the Boltzmann kinetic theory. It is
structured as follows. The Boltzmann kinetic equation for a granular binary
mixture and its associated macroscopic balance equations are introduced in
Sec. 2. Section 3 deals with the solutions to the Boltzmann equation for homo-
geneous states in the free cooling case as well as when the mixture is heated
by an external thermostat. The Chapman–Enskog method around the local
version of the homogenous distributions obtained in Sec. 2 is applied in Sec. 3
to get the form of the Navier–Stokes hydrodynamic equations. Theoretical re-
sults for the diffusion and shear viscosity transport coefficients are compared
with simulation data in Sec. 5, while the Einstein and the Onsager relations for
granular mixtures are analyzed in Secs. 6 and 7, respectively. The dispersion
relations for the hydrodynamic equations linearized about the homogeneous
cooling state are obtained in Sec. 8, showing that the homogeneous reference
state is unstable to long wavelength perturbations. The conditions for sta-
bility are identified as functions of the wave vector, the dissipation, and the
parameters of the mixture. Segregation due to thermal diffusion is studied in
Sec. 9 by using the Navier–Stokes description. A new criterion for segregation
is found that is consistent with recent experimental results. Section 10 deals
with the USF problem for a granular mixture. Finally, the paper is closed in
Sec. 11 with a discussion of the results presented here.
Before ending this section, I want to remark that the present account is a
personal perspective based on the author’s work and that of his collaborators
so that no attempt is made to include the extensive related work of many
others in this field. The references given are selective and apologies are of-
fered at the outset to the many other important contributions not recognized
explicitly.
Kinetic Theory for Binary Granular Mixtures at Low-Density 5
2 Boltzmann Kinetic Equation for Binary Mixtures of
Inelastic Hard Spheres
Consider a binary granular mixture composed by smooth inelastic disks
(d = 2) or spheres (d = 3) of massesm1 andm2, and diameters σ1 and σ2. The
inelasticity of collisions among all pairs is characterized by three independent
constant coefficients of restitution α11, α22, and α12 = α21, where αij ≤ 1 is
the coefficient of restitution for collisions between particles of species i and j.
Since the spheres are assumed to be perfectly smooth, only the translational
degrees of freedom of grains are affected by dissipation. In the low density
regime, a simultaneous interaction of more than two particles is highly unlike
and so can be neglected. Consequently, in a dilute gas the interactions among
the particles reduce to a succession of binary collisions. At this level of de-
scription, all the relevant information on the state of the system is contained
in the one-body velocity distribution functions fi(r,v; t) (i = 1, 2) defined so
that fi(r,v; t)drdv is the most probable (or average) number of particles of
species i which at time t lie in the volume element dr centered at the point
r and moving with velocities in the range dv about v. For an inelastic gas,
the distributions fi(r,v; t) (i = 1, 2) for the two species satisfy the coupled
nonlinear Boltzmann equations [30, 31]
∂t + v · ∇+
fi(r,v, t) =
Jij [v|fi(t), fj(t)] , (1)
where the Boltzmann collision operator Jij [v|fi, fj ] is
Jij [v1|fi, fj] = σ
dσ̂Θ(σ̂ · g12)(σ̂ · v12)
α−2ij fi(r,v
1, t)fj(r,v
2, t)− fi(r,v1, t)fj(r,v2, t)
. (2)
In Eq. (2), d is the dimensionality of the system, σij = (σi + σj) /2, σ̂ is an
unit vector along the line of centers, Θ is the Heaviside step function, and
v12 = v1 −v2 is the relative velocity. The primes on the velocities denote the
initial values {v′1,v
2} that lead to {v1,v2} following a binary (restituting)
collision:
1 = v1 − µji
1 + α−1ij
(σ̂ · v12)σ̂,
2 = v2 + µij
1 + α−1ij
(σ̂ · v12)σ̂, (3)
where µij ≡ mi/ (mi +mj). In addition, Fi denotes an external conservative
force acting on species i (such as a gravity field) and Fi is an operator repre-
senting a possible effect of an external nonconservative forcing which injects
energy into the system to compensate for the energy dissipated by collisional
cooling. This type of force acts as a thermostat that tries to mimics a thermal
bath. Some explicit forms for the operator Fi will be chosen later.
6 Vicente Garzó
The relevant hydrodynamic fields for the mixture are the number densities
ni, the flow velocity u, and the temperature T . They are defined in terms of
moments of the velocity distribution functions fi as
dvfi(v) , (4)
dvvfi(v) , (5)
nT = p =
niTi =
dvV 2fi(v) , (6)
where V = v − u is the peculiar velocity, n = n1 + n2 is the total number
density, ρ = m1n1 + m2n2 is the total mass density, and p is the pressure.
Furthermore, the third equality of Eq. (6) defines the kinetic temperatures Ti
for each species, which measure their mean kinetic energies.
The collision operators conserve the particle number of each species and
the total momentum, but the total energy is not conserved:
dvJij [v|fi, fj ] = 0 , (7)
dvvJij [v|fi, fj] +mj
dvvJji[v|fj , fi] = 0 , (8)
dvV 2Jij [v|fi, fj] = −dnTζ . (9)
In Eq. (9), ζ is identified as the total “cooling rate” due to collisions among
all species. It measures the rate of energy loss due to dissipation. At a kinetic
level, it is also convenient to introduce the “cooling rates” ζi for the partial
temperatures Ti. They are defined as
ζij = −
dniTi
dvV 2Jij [v|fi, fj], (10)
where the second equality defines the quantities ζij . The total cooling rate ζ
can be written in terms of the partial cooling rates ζi as
ζ = T−1
xiTiζi, (11)
where xi = ni/n is the mole fraction of species i.
From Eqs. (7)–(9), the macroscopic balance equations for the number den-
sities ni, the total momentum density ρu and the energy density (d/2)nT can
be obtained. They are given, respectively, by [19]
Kinetic Theory for Binary Granular Mixtures at Low-Density 7
Dtni + ni∇ · u+
∇ · ji
= 0 , (12)
ρDtu+∇ · P =
niFi , (13)
DtT −
∇ · ji
∇ · q+ P : ∇u−
Fi · ji
= −(ζ − ξ)T . (14)
In the above equations, Dt = ∂t + u · ∇ is the material derivative,
ji = mi
dvV fi(v) (15)
is the mass flux for species i relative to the local flow,
dvVV fi(v) (16)
is the total pressure tensor, and
dv V 2V fi(v) (17)
is the total heat flux. On the right-hand side of the temperature equation
(14), the source term ξ (measuring the rate of heating due to the external
thermostat) is given by
ξ = −
dvV 2Fifi(v). (18)
In the balance equations (12)–(14) it is assumed that the external driving
thermostat does not change the number of particles of each species or the
total momentum, i.e., ∫
dvFifi(v) = 0, (19)
dv v Fifi(v) = 0. (20)
The utility of the balance equations (12)–(14) is limited without further
specification of the fluxes and the cooling rate, which in general have a com-
plex dependence on space and time. However, for sufficiently large space and
time scales, one expects that the system reaches a hydrodynamic regime in
which all the space and time dependence is given entirely through a functional
dependence on the six hydrodynamic fields ni, u, and T . The corresponding
8 Vicente Garzó
functional dependence of ji, P, and q on these fields are called constitutive
equations and define the transport coefficients of the mixture. The primary
feature of a hydrodynamic description is the reduction of the description from
many microscopic degrees of freedom to a set of equations involving only six
local fields. At a kinetic level, the constitutive equations are obtained when
one admits the existence of a normal solution to the Boltzmann equation
where the velocity distribution functions depend on r and t only through
their functional dependence on the fields, namely,
fi(r,v1, t) = fi[v1|ni(r, t), T (r, t),u(r, t)]. (21)
This normal solution is generated by the Chapman–Enskog method [4] con-
veniently adapted to dissipative dynamics. Since the method is based on an
expansion around the local version of the homogeneous state, let us charac-
terize it before considering inhomogeneous solutions.
3 Homogeneous States
In this Section we are interested in spatially homogeneous isotropic states. In
this case, we assume that the magnitude of the conservative external forces
is at least of first order in the spatial gradients (i.e., Fi = 0), so that Eq. (1)
becomes
(∂t + Fi) fi =
Jij [fi, fj ]. (22)
For elastic collisions (αij = 1) and in the absence of external forcing (Fi = 0),
it is well known that the long-time solution to (22) is a Maxwellian distribution
for each species at the same temperature T . However, if the particles collide
inelastically (αij < 1) and Fi = 0, a steady state is not possible in uniform
situations since the temperature decreases monotonically in time. This state is
usually referred to as the homogeneous cooling state (HCS). In this situation,
since ni is uniform and u = 0, the normal (hydrodynamic) solution to fi
requires that all its time dependence occurs only through the temperature
T (t). Consequently, fi(v, t) must be of the form [15]
fi(v, t) = niv
0 (t)Φi (v/v0(t)) , (23)
where v0(t) =
2T (t)(m1 +m2)/ (m1m2) is a thermal speed defined in terms
of the temperature T (t) of the mixture. The balance equations (12)–(14) to
this order become ∂txi = ∂tu = 0, and T
−1∂tT = −ζ, where the cooling rate
ζ is determined by Eq. (9). In addition, from Eqs. (10) and (22) (when Fi = 0)
one can derive the time evolution for the temperature ratio γ = T1(t)/T2(t):
∂t ln γ = ζ2 − ζ1. (24)
Kinetic Theory for Binary Granular Mixtures at Low-Density 9
The fact that the distributions fi depend on time only through T (t) necessarily
implies that the temperature ratio γ must be independent of time and so, Eq.
(24) gives the condition
ζ1 = ζ2 = ζ. (25)
In the elastic case, where fi is a Maxwellian distribution, the above condition
yields T1(t) = T2(t) = T (t) and the energy equipartition applies. However, in
the inelastic case, the equality of the cooling rates leads to different values for
the partial temperatures, even if one considers the Maxwellian approximation
to fi. Nevertheless, the constancy of γ assures that the time dependence of
the distributions is entirely through T (t), and in fact the partial temperatures
can be expressed in terms of the global temperature as
T1(t) =
1 + x1(γ − 1)
T (t), T2(t) =
1 + x1(γ − 1)
T (t). (26)
Just as for the single gas case [30, 32], the exact form of Φi has not yet been
found, although a good approximation for thermal velocities can be obtained
from an expansion in Sonine polynomials [15]. In the leading order, Φi is given
e−θiv
∗4 − (d+ 2)θiv
d(d+ 2)
Here, v∗ ≡ v/v0,
m−1j , (28)
and γi = Ti/T . The coefficients ci (which measure the deviation of Φi from
the reference Maxwellian) are determined consistently from the Boltzmann
equation. The explicit form of the approximation (27) provides detailed pre-
dictions for the temperature ratio T1/T2 (through calculation of the cooling
rates) and for the cumulants ci as functions of the mass ratio, size ratio, com-
position and coefficients of restitution [15]. The numerical accuracy of this
truncated Sonine expansion has been confirmed by comparison with Monte
Carlo [33] and MD [34] simulations.
As said in the Introduction, the existence of different temperatures for
each species has been observed in real experiments of driven steady states.
These states are achieved from external forces that do work at the same rate
as collisional cooling. In experiments this is accomplished by vibrating the
system so that it is locally driven at walls. Far from these walls a steady state
is reached whose properties are expected to be insensitive to the details of
the driving forces. Due to the technical difficulties involved in incorporating
oscillating boundary conditions, it is usual to introduce external forces (or
thermostats) acting locally on each particle. These forces are represented by
the operator Fi in Eq. (22) and depend on the state of the system. Two
10 Vicente Garzó
types of thermostats have been usually considered in the literature. One is a
deterministic thermostat widely used in nonequilibrium MD simulations for
ordinary fluids [35, 36]. The force is similar to a Stokes law drag force, linear
in the velocity, but with the opposite sign so that it heats rather than cools
the system. In this case, Fi is given by [21, 37]
Fifi(v) =
· [vfi(v)], (29)
where the friction constant in the force has been adjusted to get a constant
temperature in the long-time limit. It must be remarked that the correspond-
ing Boltzmann equation (22) for this Gaussian thermostat force is formally
identical with the Boltzmann equation in the HCS (i.e., with Fi = 0) when
both equations are written in terms of the reduced distributions Φi(v
∗) [37].
In particular, the dependence of γ on the parameters of the system is the same
with and without the Gaussian thermostat.
A second method of driving the system is by means of a stochastic
Langevin force representing Gaussian white noise [38]. The corresponding op-
erator Fi has a Fokker–Planck form [32]
Fifi(v) = −
fi(v), (30)
where for simplicity the covariance of the stochastic acceleration has been
taken to the same for each species [39, 40]. This requirement gives the steady
state condition
ζ2. (31)
The cooling rates ζi are no longer equal, in contrast to the HCS, and the de-
pendence of γ on the control parameters is different as well [34]. The procedure
for determining the temperature ratio and the cumulants ci is the same as in
the HCS state since the steady state distribution Φi can also be represented as
an expansion of the form (27) and the coefficients are now determined from
the solution to the Boltzmann equation (22). The condition (31) gives the
corresponding equation for the temperature ratio.
Figure 1 illustrates the differences between the HCS and the stochastic
steady state at the level of the temperature ratio T1/T2. We have considered
mixtures constituted by spheres (d = 3) of the same material [and so, αij = α,
and m1/m2 = (σ1/σ2)
3] and equal volumes of large and small particles [i.e.,
x2 = (σ1/σ2)
3x1]. Here, for the sake of simplicity, the cooling rates have been
analytically estimated by using Maxwellians (namely, by taking ci = 0) for
the distributions fi(v):
ζij →
4π(d−1)/2
njµjiσ
θi + θj
×(1 + αij)
(1 + αij)
θi + θj
. (32)
Kinetic Theory for Binary Granular Mixtures at Low-Density 11
1.0 1.5 2.0 2.5 3.0
Fig. 1. Temperature ratio T1/T2 versus the size ratio σ1/σ2 for αij ≡ α = 0.78 in
the case of mixtures constituted by particles of the same material and equal total
volumes of large and small particles. The lines are the kinetic theory results in (a)
the stochastic driving case and (b) the free cooling case, while the points refer to
MD simulations [41].
Simulation data recently obtained from MD simulations in agitated mixtures
have also been included [41]. The experimental value of the (common) coef-
ficient of restitution is α = 0.78. While a good agreement between kinetic
theory and MD simulations is found when the gas is assumed to be driven
by the stochastic thermostat, significant discrepancies appear in the undriven
(HCS) case, especially as the size ratio σ1/σ2 increases. These results con-
trast with the comparison made by Brey et al. [42] for agitated mixtures in
the tracer limit (x1 → 0), where the predictions of the temperature ratio from
kinetic theory based on the condition ζ1 = ζ2 compare quite well with MD
simulations. However, for the cases studied in Ref. [42] the conditions (25) and
(31) yield quite similar results for the dependence of T1/T2 on the parameters
of the system. The good agreement found in Fig. 1 between MD simulations
for agitated mixtures and kinetic theory suggests that the stochastic driving
condition can be considered as a plausible first approximation for qualitative
comparisons with experimental results [17].
4 Navier–Stokes Hydrodynamic Equations
We consider now a spatially inhomogeneous state created either by initial
preparation or by boundary conditions. We are interested in a hydrodynamic
description where the state of the system is completely specified through their
hydrodynamic fields. This implies that the latter dominate over other excita-
12 Vicente Garzó
tions for times large compared to the mean free time and for wavelengths large
compared to the mean free time. The hydrodynamic regime is characterized
by the existence of a normal solution to the Boltzmann equation which can be
explicitly obtained by means of the Chapman–Enskog method [4]. For small
spatial variations, the functional dependence (21) of the normal solution can
be made local in space and time through an expansion in gradients of the
fields:
fi = f
i + ǫ f
i + · · · , (33)
where each factor of ǫ means an implicit gradient of a hydrodynamic field.
The reference distribution function f
i (r,v, t) is the local version of the ho-
mogeneous distribution (23), namely, it is obtained from the homogeneous
distribution by replacing the temperature, densities and flow velocity by their
nonequilibrium values:
i (r,v, t) = ni(r, t)v
0 (T (r, t))Φi (V/v0(T (r, t))) , (34)
where V = v−u(r, t). The time derivatives of the fields are also expanded as
∂t = ∂
t + ǫ∂
t + · · · . The coefficients of the time derivative expansion are
identified from the balance equations (12)–(14) with a representation of the
fluxes and the cooling rate in the macroscopic balance equations as a similar
series through their definitions as functionals of fi.
This is the usual Chapman–Enkog method for solving kinetic equations
[4, 28, 43]. Nevertheless, the complexity introduced by the energy dissipation
in collisions has led to the introduction by some authors [14] of some additional
approximations, restricting the validity of most of the results to the small
inelasticity limit. Only very recently, explicit expressions for the fluxes to first
order in the gradients as explicit functions of the coefficients of restitution have
been obtained [19, 21, 44]. To Navier–Stokes order, the constitutive equations
for mass, momentum, and heat fluxes are given, respectively, by
1 = −
m1m2n
D∇x1 −
Dp∇p−
D′∇T +
χ1iFi, j
2 = −j
kℓ = p δkℓ − η
∇ℓuk +∇kuℓ −
δkℓ∇ · u
, (36)
(1) = −T 2D′′∇x1 − L∇p− λ∇T +
κiFi. (37)
The transport coefficients in these equations are the diffusion coefficient D,
the pressure diffusion coefficient Dp, the thermal diffusion coefficient D
′, the
mobility coefficient χij , the shear viscosity η, the Dufour coefficient D
′′, the
pressure energy coefficient L, the thermal conductivity λ, and the coefficient
κi. These coefficients are defined as
Kinetic Theory for Binary Granular Mixtures at Low-Density 13
D = −
dvV ·A1, (38)
Dp = −
dvV · B1, (39)
D′ = −
dvV · C1, (40)
χij =
dvV · E ij , (41)
η = −
(d− 1)(d+ 2)
dvVV : Di, (42)
D′′ = −
dv V 2V ·Ai, (43)
L = −
dv V 2V · Bi, (44)
λ = −
dv V 2V · Ci, (45)
dv V 2V · Eij . (46)
Here, Ai(V), Bi(V), Ci(V), Di(V), and Eij(V) are functions of the peculiar
velocity and the hydrodynamic fields. They obey the following set of coupled
linear integral equations:
ζ (T∂T + p∂p) + F
1 + L1
A1 +M1A2 = A1 +
× (pB1 + TC1) , (47a)
ζ (T∂T + p∂p) + F
2 + L2
A2 +M2A1 = A2 +
× (pB2 + TC2) , (47b)
ζ (T∂T + p∂p) + F
1 + L1 − 2c
B1+M1B2 = B1+
C1, (48a)
ζ (T∂T + p∂p) + F
2 + L2 − 2c
B2+M2B1 = B2+
C2, (48b)
14 Vicente Garzó
ζ (T∂T + p∂p) + F
1 + L1 −
C1+M1C2 = C1−
B1, (49a)
ζ (T∂T + p∂p) + F
2 + L2 −
C2+M2C1 = C2−
B2, (49b)
ζ (T∂T + p∂p) + F
1 + L1
D1 +M1D2 = D1, (50a)
ζ (T∂T + p∂p) + F
2 + L2
D2 +M2D1 = D2, (50b)
ζ (T∂T + p∂p) + F
1 + L1
E11 +M1E21 = E11, (51a)
ζ (T∂T + p∂p) + F
1 + L1
E12 +M1E22 = E12. (51b)
Here, we have introduced the linearized Boltzmann collision operators
L1X = −
J11[f
1 , X ] + J11[X, f
1 ] + J12[X, f
, (52)
M1X = −J12[f
1 , X ], (53)
with a similar definition for L2 andM2. In addition, c
ζ = ζ
(0) in the undriven
case while c
ζ = 0 in the driven case, where ζ
(0) is given by Eq. (9) to zeroth
order, i.e.,
ζ(0) = −
dvV 2Jij [v|f
i , f
j ] . (54)
In Eqs. (47a)–(51b) we have also introduced the operators
i X =
· (VX) , (Gaussian thermostat), (55a)
i X = −
X, (stochastic thermostat), (55b)
and the quantities
Ai(V) = −
V, (56)
Bi(V) = −
i V +
, (57)
Ci(V) =
V, (58)
Kinetic Theory for Binary Granular Mixtures at Low-Density 15
Di(V) = V
I, (59)
Eij(V) = −
δij −
. (60)
Here, I is the unit tensor in d dimensions and ρi = mini is the mass density
of species i. Upon writing Eqs. (47a)–(51b) use has been made of the fact
that there is no contribution to the cooling rate at this order, i.e., ζ(1) = 0.
As a consequence, F
i = 0. The property ζ
(1) = 0 is special of the low
density Boltzmann kinetic theory (since f
i does not contain any contribution
proportional to ∇ ·u), but such terms occur at higher densities [45, 46]. Note
that in the particular case of the gravitational force Fi = mig (where g is
the gravity acceleration), the combination m1E11+m2E12 = 0. This leads to
Eij = 0, and so there are no contributions to the mass and heat fluxes coming
from the external conservative forces.
In summary, the solutions to the Boltzmann equations to first order in the
spatial gradients are given by [19]
fi = f
i +Ai · ∇x1 +Bi · ∇p+ Ci · ∇T +Di : ∇u+
Eij · Fj . (61)
The solution to zeroth-order is obtained from Eq. (22) while the functions
{Ai,Bi,Ci,Di,Eij} are determined from the integral equations (47a)–(51b).
Once these equations are solved, the Navier-Stokes transport coefficients are
obtained from Eqs. (38)–(46) and the mass, momentum, and heat fluxes are
explicitly known. These fluxes, together with the macroscopic balance equa-
tions (12)–(14), provide the closed set of Navier-Stokes order hydrodynamic
equations for a granular binary mixture. All these results are still formally
exact and valid for arbitrary values of the coefficients of restitution.
However, explicit expressions for the Navier–Stokes transport coefficients
require to solve the integral equations (47a)–(51b). Accurate approximations
for {Ai,Bi,Ci,Di,Eij} may be obtained using low order truncation of ex-
pansions in a series of Sonine polynomials. The polynomials are defined with
respect to a Gaussian weight factor whose parameters are chosen such that
the leading term in the expansion yields the exact moments of the entire dis-
tribution with respect to 1, v, and v2. The procedure is similar to the one
followed for elastic collisions [4] and yields explicit expressions for the trans-
port coefficients in terms of the parameters of the mixture [19, 20, 44].
5 Comparison with Monte Carlo Simulations
As said before, the expressions derived for the Navier Stokes transport coef-
ficients are obtained by considering two different approximations. First, since
the deviation of f
i from its Maxwellian form is quite small in the region of
16 Vicente Garzó
thermal velocities, one uses the distribution (27) as a trial function for f
Second, one only considers the leading terms of an expansion of the distri-
bution f
i in Sonine polynomials. Both approximations allow one to offer
a simplified kinetic theory for a granular binary mixture. To assess the de-
gree of accuracy of these predictions, one resorts to numerical solutions of the
Boltzmann equation, such as those obtained from the Direct simulation Monte
Carlo (DSMC) method [47]. Although the method was originally devised for
normal fluids, its extension to granular gases is straightforward [37]. In this
Section we provide some comparisons between theory and numerical solutions
of the Boltzmann equation by means of the DSMC method in the cases of the
diffusion coefficient D (in the tracer limit) and the shear viscosity coefficient
η of a heated gas. Let us study each coefficient separately.
5.1 Tracer Diffusion Coefficient
We consider a free granular mixture (Fi = 0) in which one of the components
of the mixture (say, for instance, species 1) is present in tracer concentration
(x1 → 0). In this situation the diffusion coefficient of impurities in a granular
gas undergoing homogeneous cooling state can be measured in simulations
from the mean square displacement of the tracer particle after a time interval
t [48]:
〈|r(t) − r(0)|2〉 =
, (62)
where |r(t) − r(0)| is the distance travelled by the impurity from t = 0 until
time t. The relation (62) written in appropriate dimensionless variables to
eliminate the time dependence of D(t) can be used to measure by computer
simulations the diffusion coefficient [8, 22].
If the hydrodynamic description (or normal solution in the context of the
Chapman–Enskogmethod) applies, then the diffusion coefficientD(t) depends
on time only through its dependence on the temperature T (t). Dimensional
analysis shows that D(t) ∝
T (t). In this case, after a transient regime,
the reduced diffusion coefficient D∗ = (m1m2/ρ)D(t)ν0(t)/T (t) achieves a
time-independent value. Here, ν0(t) = nσ
12 v0(t) ∝
T (t) is an effective
collision frequency for hard spheres. The fact that D∗ reaches a constant
value for times large compared with the mean free path is closely related
with the validity of a hydrodynamic description for the system. In addition,
as has been recently shown [49], the dependence of D∗ on the mass ratio
m1/m2 and the coefficient of restitution α12 is only through the effective
mass m∗1 = m1 + (m1 +m2)(1 − α12)/(1 + α12).
The dependence of D∗ on the common coefficient of restitution αij ≡ α
is shown in Fig. 2 in the case of hard disks (d = 2) for three different sys-
tems. The symbols refer to DSMC simulations while the lines correspond to
the kinetic theory results obtained in the first Sonine approximation [19, 44].
MD results reported in Ref. [8] when impurities and particles of the gas are
Kinetic Theory for Binary Granular Mixtures at Low-Density 17
0.6 0.7 0.8 0.9 1.0
2.5 ω=µ
µ=1/4
Fig. 2. Plot of the reduced diffusion coefficient D∗ as a function of the (common)
coefficient of restitution α for binary mixtures with ω = µ in the case of a two-
dimensional system (d = 2). Here, ω ≡ σ1/σ2 and µ ≡ m1/m2. The symbols are
computer simulation results obtained from the mean square displacement and the
lines are the theoretical results obtained in the first Sonine approximation. The
DSMC results correspond to µ = 1/4 (•), µ = 4 (◦) and µ = 1 (⋄). Molecular
dynamics results reported in Ref. [8] for µ = 1 (△) have also been included.
mechanically equivalent have also been included. We observe that in the latter
case MD and DSMC results are consistent among themselves in the range of
values of α explored. This good agreement gives support to the applicability of
the inelastic Boltzmann equation beyond the quasielastic limit. It is apparent
that the agreement between the first Sonine approximation and simulation
results is excellent when impurities and particles of the gas are mechanically
equivalent and when impurities are much heavier and/or much larger than
the particles of the gas (Brownian limit). However, some discrepancies be-
tween simulation an theory are found with decreasing values of the mass ratio
µ ≡ m1/m2 and the size ratio ω ≡ σ1/σ2. These discrepancies are not easily
observed in Fig. 2 because of the small magnitude of D∗ for µ = 1/4. For
these systems, the second Sonine approximation [22] improves the qualita-
tive predictions over the first Sonine approximation for the cases in which
the gas particles are heavier and/or larger than impurities. This means that
the Sonine polynomial expansion exhibits a slow convergence for sufficiently
small values of the mass ratio µ and/or the size ratio ω. This tendency is also
present in the case of elastic systems [50].
18 Vicente Garzó
5.2 Shear Viscosity Coefficient of a Heated Gas
The shear viscosity η is perhaps the most widely studied transport coefficient
in granular fluids. This coefficient can be measured in computer simulations
in the special hydrodynamic state of uniform shear flow (USF). At a macro-
scopic level, this state is characterized by constant partial densities ni, uni-
form temperature T , and a linear flow velocity profile u1,k = u2,k = akℓrℓ,
akℓ = aδkxδℓy, a being the constant shear rate. In this state, the temperature
changes in time due to the competition between two mechanisms: on the one
hand, viscous heating and, on the other hand, energy dissipation in collisions.
In addition, the mass and heat fluxes vanish by symmetry reasons and the
(uniform) pressure tensor is the only nonzero flux of the problem. The relevant
balance equation is that for temperature, Eq. (14), which reduces to
∂tT +
aPxy = −(ζ − ξ)T, (63)
where
Pxy =
dVVxVyfi(V) (64)
is the xy-element of the pressure tensor.
For a granular fluid under USF and in the absence of a thermostatting
force (ξ = 0), the energy balance equation (63) leads to a steady state when
the viscous heating effect is exactly balanced by the collisional cooling. This
situation will be analyzed in Sec. 10. However, if for instance the mixture
is heated by the Gaussian thermostat (29) (with v → V), then the viscous
heating still prevails so that the temperature increases in time. In this case,
the collision frequency ν0(t) ∝
T (t) also grows with t and hence the reduced
shear rate a∗(t) = a/ν0(t) (which is the relevant nonequilibrium parameter
of the problem) monotonically decreases in time. Under these conditions, the
system asymptotically achieves a regime described by linear hydrodynamics
and the (reduced) shear viscosity η∗ = [ν0(t)/nT (t)]η(t) can be measured as
η∗ = − lim
P ∗xy
, (65)
where P ∗xy = Pxy/nT . This procedure allows one to identify the shear viscosity
of a granular mixture excited by the Gaussian external force (29) and compare
it with the predictions given by the Chapman–Enskog method.
In Fig. 3, we plot the ratio η∗(α)/η∗(1) versus the mass ratiom1/m2 in the
case of hard spheres (d = 3) for σ1/σ2 = 1, x1 =
, and three different values
of the (common) coefficient of restitution αij ≡ α. Here, η
∗(1) refers to the
elastic value of the shear viscosity coefficient. Again, the symbols represent
the simulation data obtained by numerically solving the Boltzmann equation
[21], while the lines refer to the theoretical results obtained from the Boltz-
mann equation in the first Sonine approximation. We see that in general the
Kinetic Theory for Binary Granular Mixtures at Low-Density 19
0 2 4 6 8 10
η*(α)/η
Fig. 3. Plot of the ratio η∗(α)/η∗(1) as a function of the mass ratio m1/m2 for
σ1/σ2 = n1/n2 = 1 and three different values of the (common) coefficient of resti-
tution α: α = 0.9 (circles), α = 0.8 (squares), and α = 0.7 (triangles). The lines are
the theoretical predictions and the symbols refer to the results obtained from the
DSMC method.
1 2 3 4 5
η*(α)/η
Fig. 4. Plot of the ratio η∗(α)/η∗(1) as a function of the size ratio σ1/σ2 for
m1/m2 = 4, n1/n2 = 1 and three different values of the (common) coefficient of
restitution α: α = 0.9 (circles), α = 0.8 (squares), and α = 0.7 (triangles). The lines
are the theoretical predictions and the symbols refer to the results obtained from
the DSMC method.
20 Vicente Garzó
deviation of η∗(α) from its functional form for elastic collisions is quite impor-
tant. This tendency becomes more significant as the mass disparity increases.
The agreement between the first Sonine approximation and simulation is seen
to be in general excellent. This agreement is similar to the one previously
found in the monocomponent case [7, 10, 51]. At a quantitative level, the dis-
crepancies between theory and simulation tend to increase as the coefficient
of restitution decreases, although these differences are quite small (say, for
instance, around 2% at α = 0.7 in the disparate mass case m1/m2 = 10).
The influence of the size ratio on the shear viscosity is shown in Fig. 4 for
m1/m2 = 4 and x1 =
[21]. We observe again a strong dependence of the
shear viscosity on dissipation. However, for a given value of α, the influence of
σ1/σ2 on η
∗ is weaker than the one found before in Fig. 3 for the mass ratio.
The agreement for both α = 0.9 and α = 0.8 is quite good, except for the
largest size ratio at α = 0.8. These discrepancies become more significant as
the dissipation increases, especially for mixtures of particles of very different
sizes. In summary, according to the comparison carried out in Figs. 3 and 4,
one can conclude that the agreement between theory and simulation extends
over a wide range values of the coefficient of restitution, indicating the relia-
bility of the first Sonine approximation for describing granular flows beyond
the quasielastic limit.
6 Einstein Relation in Granular Gases
The results presented in Section 5 give some support to the validity of the
hydrodynamic description to granular fluids. However, in spite of this sup-
port some care is warranted in extending properties of normal fluids to those
with inelastic collisions. Thus, for elastic collisions, in the case of an impurity
(tracer) particle immersed in a gas the response to an external force on the
impurity particle leads to a mobility coefficient proportional to the diffusion
coefficient. This is the usual Einstein relation [48], which is a consequence of
the fluctuation-dissipation theorem. A natural question is whether the Ein-
stein relation also applies for granular fluids.
To analyze it, let us consider the tracer limit (x1 → 0) and assume that
the current of impurities j
1 is only generated by the presence of a weak
concentration gradient ∇x1 and/or a weak external field F1 acting only on
the impurity particles. Under these conditions, Eq. (35) becomes
1 = −m1D∇x1 + χ11F1. (66)
The Einstein ratio ǫ′ between the diffusion coefficient D and the mobility
coefficient χ11 is defined as
ǫ′ = m1x1
, (67)
Kinetic Theory for Binary Granular Mixtures at Low-Density 21
where T ≃ T2 in the tracer limit. For elastic collisions, the Chapman–Enskog
results yield ǫ′ = 1. However, at finite inelasticity the relationship between D
and χ11 is no longer simple and, as expected, the Chapman–Enskog expres-
sions for D and χ11 in the case of an unforced granular gas [23] clearly show
that ǫ′ 6= 1. This means that the Einstein relation does not apply in granular
gases. The deviations of the (standard) Einstein ratio ǫ′ from unity has three
distinct origins: the absence of the Gibbs state (non-Gaussianity of the distri-
bution function of the HCS), time evolution of the granular temperature, and
the occurrence of different kinetic temperatures between the impurity and gas
particles. The second source of discrepancy can be avoided if the system is
driven by an external energy input to achieve a stationary state. With respect
to the third reason of violation, this could also be partially eliminated if the
temperature of the gas T is replaced by the temperature of the impurity T1
in the usual Einstein relation (67). This change yields the modified Einstein
ratio
ǫ = m1x1
T1χ11
. (68)
As a consequence, the only reason for which ǫ 6= 1 is due to the non-Maxwellian
behavior of the HCS distribution. Given that the deviations of the gas dis-
tribution f
2 from its Maxwellian form are small [15], the discrepancies of ǫ
from unity could be difficult to detect in computer simulations. This conclu-
sion agrees with recent MD simulations [52] of granular mixtures subjected to
the stochastic driving of the form (30), where no deviations from the (modi-
fied) Einstein relation ǫ = 1 have been observed for a wide range of values of
the coefficients of restitution and parameters of the system.
To illustrate the influence of dissipation on the Einstein ratio more gener-
ally, in Figs. 5 and 6 the Einstein ratio as given by (68) is plotted versus the
coefficient of restitution α12 for σ1/σ2 = 1 and different values of the mass ra-
tio m1/m2 and the coefficient of restitution α22. The results obtained by using
the Gaussian thermostat (29) are shown in Fig. 5, while Fig. 6 corresponds
to the results derived when the system is heated by the stochastic thermostat
(30) [24]. We observe that in general ǫ 6= 1, although its value is very close to
unity, especially in the case of the stochastic thermostat, where the deviations
from the Einstein relation are smaller than 1%. However, in the case of the
Gaussian thermostat the deviations from unity are about 8%, which could be
detected in computer simulations. Figures 5 and 6 also show the fact that the
transport properties are affected by the thermostat introduced so that the
latter does not play a neutral role in the problem [51].
7 Onsager’s Reciprocal Relations in Granular Gases
In the usual language of the linear irreversible thermodynamics for ordinary
fluids [53], the constitutive equations (35) and (37) for the mass flux and heat
flux in the absence of external forces can be written as
22 Vicente Garzó
0.0 0.2 0.4 0.6 0.8 1.0
Fig. 5. Plot of the modified Einstein ratio ǫ versus the coefficient of restitution α12
for the Gaussian thermostat in the cases: (a) α22 = α12, m1/m2 = 5 and σ1/σ2 = 1;
(b) α22 = α12, m1/m2 = 0.5 and σ1/σ2 = 1; and (c) α22 = 0.5, m1/m2 = 10 and
σ1/σ2 = 1.
0.0 0.2 0.4 0.6 0.8 1.0
(c) (b)
Fig. 6. Plot of the Einstein ratio ǫ versus the coefficient of restitution α12 for the
stochastic thermostat in the cases: (a) α22 = α12, m1/m2 = 5 and σ1/σ2 = 1; (b)
α22 = α12, m1/m2 = 0.5 and σ1/σ2 = 1; and (c) α22 = 0.5, m1/m2 = 10 and
σ1/σ2 = 1.
Kinetic Theory for Binary Granular Mixtures at Low-Density 23
i = −
− Liq
− Cp∇p, (69)
q = −Lqq∇T −
− C′p∇p, (70)
where
q ≡ q
(1) −
and (
∇ ln(xip), (72)
µi being the chemical potential per unit mass. In Eqs. (69) and (70), the
coefficients Lij are the so-called Onsager phenomenological coefficients and the
coefficients Cp and C
p can be expressed in terms of the transport coefficients
associated with the heat and mass fluxes. For elastic fluids, Onsager showed
[53] that time reversal invariance of the underlying microscopic equations of
motion implies important constraints on the above set of transport coefficients,
namely
Lij = Lji, Liq = Lqi, Cp = C
p = 0. (73)
The first two symmetries are called reciprocal relations as they relate trans-
port coefficients for different processes. The last two are statements that the
pressure gradient does not appear in any of the fluxes, even though it is ad-
mitted by symmetry. Even for a one component fluid, Onsager’s theorem is
significant as it leads to a new contribution to the heat flux proportional to
the density gradient [5]. Since there is no time reversal symmetry for granular
fluids, Eqs. (73) cannot be expected to apply. However, since explicit expres-
sions for all transport coefficients are at hand, the quantitative extent of the
violation can be explored.
To make connection with the expressions (35) and (37) for the mass and
heat fluxes, respectively, it is first necessary to transform Eqs. (69)–(71) to
the variables x1, p, T. Since ∇x1 = −∇x2, Eq. (72) implies
(∇µ1)T − (∇µ2)T
∇x1 +
(m2 −m1)∇ ln p
. (74)
The coefficients {Lij, Liq, Lqi, Lqq, Cp, C
p} then can be easily obtained in
terms of the Navier-Stokes transport coefficients introduced in Sec. 4. The
result is [20]
L11 = −L12 = −L21 =
m1m2ρ1ρ2
D, L1q = ρTD
′, (75)
Lq1 = −Lq2 =
T 2ρ1ρ2
D′′ −
Tρ1ρ2
(m2 −m1)D, (76)
24 Vicente Garzó
Lqq = λ−
m2 −m1
D′, (77)
(m2 −m1)D, (78)
C′p ≡ L−
m2 −m1
T 2(m2 −m1)D
′′. (79)
Onsager’s relation L12 = L21 holds since the diffusion coefficient D is
symmetric under the change 1 ↔ 2 [19]. However, in general L1q 6= Lq1, Cp 6=
0, and C′p 6= 0 [20]. The Chapman-Enskog results [19] show that there are only
two limit cases for which L1q−Lq1 = Cp = C
p = 0: (i) the elastic limit (αij =
1) with arbitrary values of masses, sizes and composition and (ii) the case
of mechanically equivalent particles with arbitrary values of the (common)
coefficient of restitution α ≡ αij . Beyond these limit cases, Onsager’s relations
do not apply. At macroscopic level the origin of this failure is due to the cooling
of the reference state as well as the occurrence of different kinetic temperatures
for both species.
8 Linearized Hydrodynamic Equations and Stability of
the Homogeneous Cooling State
As shown in Sec. 4, the Navier-Stokes constitutive equations (35)–(37) have
been expressed in terms of a set of experimentally accessible fields such as
the composition of species 1, x1, the pressure p, the mean flow field u, and
the granular temperature T . In terms of these variables and in the absence of
external forces, the macroscopic balance equations (12)–(14) become
Dtx1 +
n2m1m2
∇ · j1 = 0 , (80)
Dtp+ p∇ · u+
(∇ · q+ P : ∇u) = −ζp, (81)
Dtu+ ρ
−1∇ · P = 0 , (82)
DtT −
∇ · ji
(∇ · q+ P : ∇u) = −ζT . (83)
When the expressions (35)–(37) for the fluxes and the cooling rate ζ → ζ(0) are
substituted into the above exact balance equations (80)–(83) one gets a closed
set of hydrodynamic equations for x1,u, T, and p. These are the Navier–Stokes
hydrodynamic equations for a binary granular mixture:
Dtx1 =
n2m1m2
m1m2n
D∇x1 +
Dp∇p+
, (84)
Kinetic Theory for Binary Granular Mixtures at Low-Density 25
(Dt + ζ) p+
p∇ · u =
T 2D′′∇x1 + L∇p+ λ∇T
∇ℓuk +∇kuℓ −
δkℓ∇ · u
∇ℓuk,
(Dt + ζ) T +
p∇ · u = −
m2 −m1
m1m2n
D∇x1 +
T 2D′′∇x1 + L∇p+ λ∇T
∇ℓuk +∇kuℓ −
δkℓ∇ · u
∇ℓuk, (86)
Dtuℓ + ρ
−1∇ℓp = ρ
−1∇kη
∇ℓuk +∇kuℓ −
δkℓ∇ · u
. (87)
For the chosen set of fields, n = p/T and ρ = p [(m1 −m2)x1 +m2] /T . These
equations are exact to second order in the spatial gradients for a low density
Boltzmann gas. Note that in Eqs. (84)–(87) the second order contributions
to the cooling rate have been neglected. These second order terms have been
calculated for a monocomponent fluid [5] and found to be very small relative
to corresponding terms from the fluxes. Consequently, they have not been
considered in the hydrodynamic equations (84)–(87).
One of the main peculiarities of the granular gases (in contrast to ordinary
fluids) is the existence of non-trivial solutions to the Navier–Stokes equations
(84)–(87), even for spatially homogeneous states,
∂tx1H = 0, ∂tuH = 0, (88)
[∂t + ζ (x1H , TH , pH)]TH = 0, [∂t + ζ (x1H , TH , pH)] pH = 0, (89)
where the subscript H denotes the homogeneous state. Since the dependence
of the cooling rate ζ (x1H , TH , pH) on x1H , TH , pH is known [15, 44], these
first order nonlinear equations can be solved for the time dependence of the
homogeneous state. The result is the familiar Haff cooling law for T (t) at
constant density [1, 54]:
TH(t) =
TH(0)
[1 + ζ(0)t/2]
. (90)
As said before, each partial temperature Ti(t) has the same time dependence
but with a different value [15],
T1H(t) =
1 + x1(γ − 1)
TH(t), T2H(t) =
1 + x1(γ − 1)
TH(t), (91)
where γ = T1H(t)/T2H(t) is the time-independent temperature ratio.
26 Vicente Garzó
Nevertheless, the homogeneous cooling state (HCS) is unstable to suffi-
ciently long wavelength perturbations. For systems large enough to support
such spontaneous fluctuations, the HCS becomes inhomogeneous at long times.
This feature was first observed in MD simulations of free monocomponent
gases [55]. In MD simulations the inhomogeneities may grow by the forma-
tion of clusters, ultimately aggregating to a single large cluster [56]; if cluster
growth is suppressed, a vortex field may grow to the system size where peri-
odic boundary conditions can induce a transition to a state with macroscopic
shear. The mechanism responsible for the growth of inhomogeneities can be
understood at the level of the Navier–Stokes hydrodynamics, where linear
stability analysis shows two shear modes and a heat mode to be unstable
[5, 7, 54, 57].
The objective here is to extend this analysis to the case of a binary mixture.
To do that, we perform a linear stability analysis of the nonlinear hydrody-
namic equations (84)–(87) with respect to this HCS for small initial spatial
perturbations. For ordinary fluids such perturbations decay in time according
to the hydrodynamic modes of diffusion (shear, thermal, mass) and damped
sound propagation [3, 58, 59]. For inelastic collisions, the analysis is for fixed
coefficients of restitution in the long wavelength limit. As will be seen be-
low, the corresponding modes for a granular mixture are then quite different
from those for ordinary mixtures. In fact, an alternative study with fixed long
wavelength and coefficients of restitution approaching unity yields the usual
ordinary fluid modes. Consequently, the nature of the hydrodynamic modes
is non uniform with respect to the inelasticity and the wavelength of the
perturbation.
Let us assume that the deviations δyα(r, t) = yα(r, t) − yHα(t) are small.
Here, δyα(r, t) denotes the deviation of {x1,u, T, p} from their values in the
HCS. If the initial spatial perturbation is sufficiently small, then for some
initial time interval these deviations will remain small and the hydrodynamic
equations (84)–(87) can be linearized with respect to δyα(r, t). This leads to
a set of partial differential equations with coefficients that are independent
of space but which depend on time. As in the monocomponent case [5, 57],
this time dependence can be eliminated through a change in the time and
space variables, and a scaling of the hydrodynamic fields. We introduce the
following dimensionless space and time variables:
dt′ν0H(t
′), s =
ν0H(t)
v0H(t)
r, (92)
where ν0H(t) is an effective collision frequency for hard spheres and v0H =√
2TH(m1 +m2)/m1m2. Since {x1H ,uH , TH , pH} are evaluated in the HCS,
then Eqs. (88) and (89) hold. A set of Fourier transformed dimensionless
variables are then introduced as
ρk(τ) =
δx1k(τ)
, wk(τ) =
δuk(τ)
v0H(τ)
, θk(τ) =
δTk(τ)
TH(τ)
, Πk(τ) =
δpk(τ)
pH(τ)
Kinetic Theory for Binary Granular Mixtures at Low-Density 27
where δyαk ≡ {δx1k, δuk, δTk, δpk} is defined as
δyαk(τ) =
ds e−ik·sδyα(s, τ). (94)
Note that here the wave vector k is dimensionless.
In terms of the above variables, the transverse velocity components wk⊥ =
wk − (wk · k̂)k̂ (orthogonal to the wave vector k) decouple from the other
four modes and hence can be obtained more easily. They obey the equation
+ η∗k2
wk⊥ = 0, (95)
where ζ∗ = ζH/ν0H and
η, (96)
where ρH = m1n1H +m2n2H . The solution for wk⊥(τ) reads
wk⊥(τ) = wk⊥(0) exp[s⊥(k)τ ], (97)
where
s⊥(k) =
ζ∗ − η∗k2. (98)
This identifies d − 1 shear (transversal) modes. We see from Eq. (98) that
there exists a critical wave number kc⊥ given by
kc⊥ =
. (99)
This critical value separates two regimes: shear modes with k ≥ kc⊥ always
decay while those with k < kc⊥ grow exponentially.
The remaining modes are called longitudinal modes. They correspond to
the set {ρk, θk, Πk, wk||} where the longitudinal velocity component (parallel
to k) is wk|| = wk · k̂. These modes are the solutions of the linear equation
∂δzαk(τ)
αβ + ikM
αβ + k
δzβk(τ), (100)
where δzαk(τ) denotes now the four variables
ρk, θk, Πk, wk||
. The matrices
in Eq. (100) are given by
(0) =

0 0 0 0
ζ∗ −ζ∗ 0
ζ∗ −ζ∗ 0
0 0 0 1

, (101)
28 Vicente Garzó
(1) =
0 0 0 0
0 0 0 − 2
0 0 0 − d+2
0 0 − 1
x1µ+x2
 , (102)
(2) =
−D∗ −x−11 D
′∗ −x−11 D
′′∗ − 2
λ∗ − 2
0 0 0 − 2
(d− 1)η∗
 , (103)
where
21 = −x1
D′′∗ −
x1µ+ x2
, (104)
x1µ+ x2
D′∗ −
λ∗, (105)
23 = −
x1µ+ x2
D∗p. (106)
In these equations, µ = m1/m2, xi = niH/nH , and we have introduced the
reduced Navier–Stokes transport coefficients1
D, D∗p =
ρ2Hν0H
m1m2n
Dp, D
ρ2Hν0H
m1m2n
D′, (107)
D′′∗ =
ν0HTH
D′′, L∗ =
L, λ∗ =
λ. (108)
The longitudinal modes have the form exp[sn(k)τ ] with n = 1, 2, 3, 4,
where sn(k) are the eigenvalues of the matrix M(k) = M
(0)+ ikM(1)+k2M(2),
namely, they are the solutions of the quartic equation
det |M− sI| = 0. (109)
The solution to (109) for arbitrary values of k is quite intricate. It is instructive
to consider first the solutions to these equations in the extreme long wave-
length limit, k = 0. In this case, they are found to be the eigenvalues of the
matrix of M(0):
s(0)n =
0, 0,−
. (110)
Hence, at asymptotically long wavelengths (k = 0) the spectrum of the lin-
earized hydrodynamic equations (both transverse and longitudinal) is com-
prised of a decaying mode at −ζ∗/2, a two-fold degenerate mode at 0, and a
d-fold degenerate unstable mode at ζ∗/2. Consequently, some of the solutions
are unstable. The two zero eigenvalues represent marginal stability solutions,
1 Note that the definition for the reduced diffusion coefficient D∗ given here differs
from the one introduced in Sec. 5.1.
Kinetic Theory for Binary Granular Mixtures at Low-Density 29
0.0 0.1 0.2 0.3 0.4 0.5
-0.10
-0.05
Fig. 7. Dispersion relations for α = 0.9, x1 = 0.2, ω = 1 and µ = 4.
while the negative eigenvalue gives stable solutions. For general initial per-
turbations all modes are excited. These modes correspond to evolution of
the fluid due to uniform perturbations of the HCS, i.e., a global change in
the HCS parameters. The unstable modes are seen to arise from the initial
perturbations wk⊥(0) or wk||(0). The marginal modes correspond to changes
in the composition at fixed pressure, density, and velocity, and to changes in
Πk−θk at constant composition and velocity. The decaying mode corresponds
to changes in the temperature or pressure for Πk = θk. The unstable modes
may appear trivial since they are due entirely to the normalization of the
fluid velocity by the time dependent thermal velocity. However, this normal-
ization is required by the scaling of the entire set of equations to obtain time
independent coefficients.
The real parts of the modes s⊥(k) and sn(k) is illustrated in Fig. 7 in the
case of hard spheres (d = 3) for α ≡ αij = 0.9, σ1/σ2 = 1, x1 = 0.2, and
m1/m2 = 4. The k = 0 values correspond to five hydrodynamic modes with
two different degeneracies. The shear mode degeneracy remains at finite k but
the other is removed at any finite k. At sufficiently large k a pair of real modes
become equal and become a complex conjugate pair at all larger wave vectors,
like a sound mode. The smallest of the unstable modes is that associated with
the longitudinal velocity, which couples to the scalar hydrodynamic fields. It
becomes negative at a wave vector smaller than that of Eq. (99) and gives the
threshold for development of spatial instabilities.
The results obtained here for mixtures show no new surprises relative to
the case for a monocomponent gas [5, 54, 57], with only the addition of the
stable mass diffusion mode. Of course, the quantitative features can be quite
different since there are additional degrees of freedom with the parameter set
{x1H ,m1/m2, σ1/σ2, αij}. Also, the manner in which these linear instabilities
30 Vicente Garzó
are enhanced by the nonlinearities may be different from that for the one
component case [60].
9 Segregation in Granular Binary Mixtures: Thermal
Diffusion
The analysis of the linearized hydrodynamic equations for a granular binary
mixture has shown that the resulting equations exhibit a long wavelength
instability for d of the modes. These instabilities lead to the spontaneous for-
mation of velocity vortices and density clusters when the system evolves freely.
A phenomenon related with the density clustering is the separation or species
segregation. Segregation and mixing of dissimilar grains is perhaps one of the
most interesting problems in agitated granular mixtures. In some processes
it is a desired and useful effect to separate particles of different types, while
in other situations it is undesired and can be difficult to control. A variety
of mechanisms have been proposed to describe the separation of particles of
two sizes in a mixture of vertically shaken particles. Different mechanisms
include void filling, static compressive force, convection, condensation, ther-
mal diffusion, interstitial gas forcing, friction, and buoyancy [61]. However,
in spite of the extensive literature published in the past few years on this
subject, the problem is not completely understood yet. Among the different
competing mechanisms, thermal diffusion becomes one of the most relevant
at large shaking amplitude where the sample of macroscopic grains resembles
a granular gas. In this regime, binary collisions prevail and kinetic theory can
be quite useful to analyze the physical mechanisms involved in segregation
processes.
Thermal diffusion is caused by the relative motion of the components of a
mixture because of the presence of a temperature gradient. Due to this motion,
concentration gradients subsequently appear in the mixture producing diffu-
sion that tends to oppose those gradients. A steady state is finally achieved
in which the separation effect arising from thermal diffusion is compensated
by the diffusion effect. In these conditions, the so-called thermal diffusion fac-
tor Λij characterizes the amount of segregation parallel to the temperature
gradient. In this Section, the thermal diffusion factor is determined from the
Chapman–Enskog solution described before.
To make some contact with experiments, let us assume that the binary
granular mixture is in the presence of the gravitational field g = −gêz, where g
is a positive constant and êz is the unit vector in the positive direction of the z
axis. In experiments [41], the energy is usually supplied by vibrating horizontal
walls so that the system reaches a steady state. Here, instead of considering
oscillating boundary conditions, particles are assumed to be heated by the
action of the stochastic driving force (30), which mimics a thermal bath.
As said above, although the relation between this driven idealized method
with the use of locally driven wall forces is not completely understood, it
Kinetic Theory for Binary Granular Mixtures at Low-Density 31
must be remarked that in the case of boundary conditions corresponding to
a sawtooth vibration of one wall the condition to determine the temperature
ratio coincides with the one derived from the stochastic force [34]. The good
agreement between theory and simulation found in Fig. 1 for the temperature
ratio confirms this expectation.
The thermal diffusion factor Λij (i 6= j) is defined at the steady state in
which the mass fluxes ji vanish. Under these conditions, the factor Λij is given
through the relation [62]
− Λij∇ lnT =
∇xi, Λij + Λji = 0. (111)
The physical meaning of Λij can be described by considering a granular binary
mixture held between plates at different temperatures T (top plate) and T ′
(bottom plate) under gravity. For the sake of concreteness, we will assume
that gravity and thermal gradient point in parallel directions, i.e., the bottom
is hotter than the top (T ′ > T ). In addition, without loss of generality, we
also assume that σ1 > σ2. In the steady state, Eq. (111) describes how the
thermal field is related to the composition of the mixture. Assuming that Λ12 is
constant over the relevant ranges of temperature and composition, integration
of Eq. (111) yields
= Λ12 ln
, (112)
where xi refers to the mole fraction of species i at the top plate and x
i refers
to the mole fraction of species i at the bottom plate. Consequently, according
to Eq. (112), if Λ12 > 0, then x
1 < x1, while if Λ12 < 0, then x
1 > x1. In
summary, when Λ12 > 0, the larger particles accumulate at the top of the
sample (cold plate), while if Λ12 < 0, the larger particles accumulate at the
bottom of the sample (hot plate). The former situation is referred to as the
Brazil-nut effect (BNE) while the latter is called the reverse Brazil-nut effect
(RBNE).
The RBNE was first observed by Hong et al. [63] in MD simulations of
vertically vibrated systems. They proposed a very simple segregation criterion
that was later confirmed by Jenkins and Yoon [64] by using kinetic theory.
More recently, Breu et al. [65] have experimentally investigated conditions
under which the large particles sink to the bottom and claim that their ex-
periments confirm the theory of Hong et al. [63] provided a number of condi-
tions are chosen carefully. In addition to the vertically vibrated systems, some
works have also focused in the last few years on horizontally driven systems
showing some similarities to the BNE and its reverse form [66]. However, it
is important to note that the criterion given in Ref. [63] is based on some
drastic assumptions: elastic particles, homogeneous temperature, and energy
equipartition. These conditions preclude a comparison of the kinetic theory
derived here with the above simulations.
Some theoretical attempts to assess the influence of non-equipartition on
segregation have been recently published. Thus, Trujillo et al. [67] have derived
32 Vicente Garzó
0.5 0.6 0.7 0.8 0.9 1.0
Fig. 8. Phase diagram for BNE/RBNE for mixtures constituted by spheres (d = 3)
of the same mass density and equal total volumes of large and small particles. The
data points represent the MD simulation results [41] for α = 0.78 when convection
is suppressed. Points below (above) the curve correspond to RBNE (BNE).
an evolution equation for the relative velocity of the intruders starting from
the kinetic theory proposed by Jenkins and Yoon [64], which applies for weak
dissipation. They use constitutive relations for partial pressures that take into
account the breakdown of energy equipartition between the two species. How-
ever, the influence of temperature gradients, which exist in the vibro-fluidized
regime, is neglected in Ref. [67] because it is assumed that the pressure and
temperature are constant in the absence of the intruder. A more refined theory
has recently been provided by Brey et al. [42] in the case of a single intruder
in a vibrated granular mixture under gravity. The theory displayed in this
section covers some of the aspects not accounted for in the previous theories
[42, 64, 67] since it is based on a kinetic theory [19] that goes beyond the
quasi-elastic limit [64, 67] and applies for arbitrary composition x1 (and so,
it reduces to the results obtained in Ref. [42] when x1 → 0). This allows one
to assess the influence of composition and dissipation on thermal diffusion in
bi-disperse granular gases without any restriction on the parameter space of
the system.
To determine the dependence of the coefficient Λ12 on the parameters of
the mixture, we consider a non-convecting (u = 0) steady state with only
gradients along the vertical direction (z axis). In this case, the mass balance
equation (12) yields j1 = j2 = 0, while the momentum equation (13) gives
= −ρg. (113)
To first order in the spatial gradients, the constitutive equation for the mass
flux j1,z is given by Eq. (35), i.e.,
Kinetic Theory for Binary Granular Mixtures at Low-Density 33
j1,z = −
m1m2n
, (114)
where the susceptibility coefficient χij = 0 in the particular case of the grav-
itational force. The condition j1,z = 0 yields
m1m2np
m1m2p
, (115)
where use has been made of Eq. (113). Substitution of Eq. (115) into Eq.
(111) leads to
Λ12 =
D′ −Dpg
, (116)
where
n (∂T/∂z)
< 0 (117)
is the reduced gravity acceleration. Since the mutual diffusion coefficient D
is positive [19, 44], the sign of Λ12 is determined by the sign of the quantity
D′ −Dpg
∗. This result is general since it goes beyond the regime of density
considered.
1 2 3 4 5
Fig. 9. Phase diagram for BNE/RBNE in three dimensions for αij = 0.7 and three
values of composition: (a) x1 = 0, (b) x1 = 0.3, and (c) x1 = 0.7. Points below
(above) each curve correspond to RBNE (BNE).
To gain some insight into the explicit dependence of D′ and Dp on the
parameter space of the system, one has to resort to a kinetic theory descrip-
tion. For a low-density gas, the expressions of the coefficients D′ and Dp in
the first-Sonine approximation are given by
34 Vicente Garzó
1 2 3 4
Fig. 10. Phase diagram for BNE/RBNE in three dimensions for x1 =
and three
values of the (common) coefficient of restitution: (a) α = 0.9, (b) α = 0.8, and (c)
α = 0.5. Points below (above) each curve correspond to RBNE (BNE).
D′ = 0, Dp =
x2 + x1γ
, (118)
where ν is the (positive) collision frequency [44]
2π(d−1)/2
) nσd−112 v0(1 + α12)
θ1 + θ2
(x2µ21 + x1µ12) . (119)
Given that the driving stochastic term does not play a neutral role in the
transport, it must be remarked that the expressions for the transport coeffi-
cients obtained in the driven case slightly differ from the ones derived in the
free cooling case [19, 44].
Consequently, according to Eqs. (116) and (118), the sign of Λ12 is the
same as that of the pressure diffusion coefficient Dp. The condition Λ12 = 0
(or equivalently, Dp = 0) provides the criterion for the transition from BNE
to RBNE. Equation (118) shows that the sign of Dp is determined by the
value of the control parameter
. (120)
This parameter gives the mean square velocity of the large particles relative
to that of the small particles. Thus, when θ > 1 (θ < 1), the thermal diffusion
factor is positive (negative), which leads to BNE (RBNE). The criterion for
the transition condition from BNE to RBNE is θ = 1, i.e.,
. (121)
Kinetic Theory for Binary Granular Mixtures at Low-Density 35
In the case of equal granular temperatures (energy equipartition), θ → µ−1
and so segregation is predicted for particles that differ in mass, no matter
what their diameters may be [64]. It must be remarked that, due to the lack
of energy equipartition, the condition θ = 1 is rather complicated since it
involves all the parameter space of the system. In particular, even when the
species differ only by their respective coefficients of restitution they also seg-
regate when subject to a temperature gradient. This is a novel pure effect of
inelasticity on segregation [25, 68]. On the other hand, the criterion (121) for
the transition BNE⇐⇒RBNE is the same as the one found previously in Ref.
[67] when αij is close to 1 and in Ref. [42] in the intruder limit case (x1 → 0).
However, as said before, the results obtained here are more general since they
cover all the range of the parameter space of the system.
To illustrate size segregation driven by thermal diffusion, we consider mix-
tures constituted by spheres (d = 3) of the same material and equal total
volumes of large and small particles. In this case, m1/m2 = (σ1/σ2)
3 and
x2/x1 = (σ1/σ2)
3. Figure 8 shows the phase diagram BNE/RBNE for this
kind of systems. The data points represent the simulation results obtained by
Schröter et al. [41] for α = 0.78 in agitated mixtures constituted by parti-
cles of the same density. To the best of my knowledge, this is one of the few
experiments in which thermal diffusion has been isolated from the remain-
ing segregation mechanisms [61]. Our results show that, for a given value of
the coefficient of restitution, the RBNE is dominant at small diameter ra-
tios. However, since non-equipartition grows with increasing diameter ratio,
the system shows a crossover to BNE at sufficiently large diameter ratios.
This behavior agrees qualitatively well with the results reported in Ref. [41]
at large shaking amplitudes, where thermal diffusion becomes the relevant
segregation mechanism. At a quantitative level, we observe that the results
are also consistent with the simulation results reported in [41] when periodic
boundary conditions are used to suppress convection since they do not ob-
serve a change back to BNE for diameter ratios up to 3 (see red squares in
Fig. 11 of [41]). Although the parameter range explored in MD simulations is
smaller than the one analyzed here, one is tempted to extrapolate the simula-
tion data presented in Ref. [41] to roughly predict the transition value of the
diameter ratio at α = 0.78 (which is the value of the coefficient of restitution
considered in the simulations). Thus, if one extrapolates from the simulation
data at the diameter ratios of σ1/σ2 = 2 and σ1/σ2 = 3, one sees that the
transition from RBNE to BNE might be around σ1/σ2 = 10, which would
quantitatively agree with the results reported in Fig. 8. Figure 8 also shows
that the BNE is completely destroyed in the quasielastic limit (α ≃ 1).
Let us now investigate the influence of composition on segregation. Figure 9
shows a typical phase diagram in the three-dimensional case for αij ≡ α = 0.7
and three different values of the mole fraction x1. The lines separate the
regimes between BNE and RBNE. We observe that the composition of the
mixture has significant effects in reducing the BNE as the concentration of
larger particles increases. In addition, for a given value of composition, the
36 Vicente Garzó
transition from BNE to RBNE may occur following two paths: (i) along the
constant mass ratio m2/m1 with increasing size ratio σ1/σ2, and (ii) along
the constant size ratio with increasing mass ratio m2/m1. The influence of
dissipation on the phase diagrams BNE/RBNE is illustrated in Fig. 10 for
d = 3 in the case of an equimolar mixture (x1 =
) and three values of
the (common) coefficient of restitution α. We observe that the role played by
inelasticity is quite important since the regime of RBNE increases significantly
with dissipation. Similar results are found for other values of composition.
In summary, thermal diffusion (which is the relevant segregation mech-
anism in agitated granular mixtures at large shaking amplitudes) can been
analyzed by the Boltzmann kinetic theory. This theory is able to explain
some of the experimental and/or MD segregation results [41] observed within
the range of parameter space explored. A more quantitative comparison in
the dilute regime with MD simulations is needed to show the relevance of
the Boltzmann equation to analyze segregation driven by a thermal gradi-
ent. As said before, comparison with MD simulations in the tracer limit case
(x1 → 0) [42] for a dilute gas has shown the reliability of the inelastic Boltz-
mann equation to describe segregation. In this context, one expects that the
same agreement observed before in the intruder case [42] is maintained when
x1 is different from zero.
10 Steady States: Uniform Shear Flow
In the preceding sections, the Navier–Stokes equations (constitutive equa-
tions that are linear in the hydrodynamic gradients) have been shown to be
quite useful to describe appropriately several problems in granular mixtures.
However, under some circumstances large gradients occur and more complex
constitutive equations are required. The need for more complex constitutive
equations does not signal a breakdown of hydrodynamics [69], only a failure
of the Navier–Stokes approximation [70]. Although in this case the Chapman-
Enskog method can be carried out to second order in gradients (Burnett or-
der), it is likely that failure of the Navier-Stokes description signals the need
for other methods to construct the normal solution that are not based on a
small gradient expansion.
One of the most interesting problems in granular fluids is the simple or
uniform shear flow (USF) [13, 26]. As described in Sec. 5, this state is charac-
terized by uniform density and temperature and a simple shear with the local
velocity field given by u1,x = u2,x = ay, uy = uz = 0, where a is the constant
shear rate. The USF is a well-known nonequilibrium problem widely studied,
for both granular [26, 27, 71, 72, 73] and conventional [28, 74] gases. However,
the nature of this state is quite different in both systems since a steady state
is achieved for granular fluids when viscous (shear) heating is compensated
for by energy dissipation in collisions:
Kinetic Theory for Binary Granular Mixtures at Low-Density 37
aPxy = −
nTζ. (122)
This steady state is what we want to analyze in this section. The balance
equation (122) shows the intrinsic connection between the shear field and dis-
sipation in the system. This contrasts with the description of USF for elastic
fluids where a steady state is not possible unless an external thermostat is in-
troduced [28]. Note that the hydrodynamic steady shear flow state associated
with the condition (122) is inherently beyond the scope of the Navier–Stokes
or Newtonian hydrodynamic equations [29]. The reason for this is the exis-
tence of an internal mechanism, collisional cooling, that sets the strength of
the velocity gradient in the steady state. For normal fluids, this scale is set by
external sources (boundary conditions, driving forces) that can be controlled
to admit the conditions required for Navier–Stokes hydrodynamics. In con-
trast, collisional cooling is fixed by the mechanical properties of the particles
making up the fluid. This observation is significant because it prevents the
possibility of measuring the Newtonian shear viscosity for granular fluids in
the steady USF [72, 73]. More generally, it provides a caution regarding the
simulation of other steady states to study Navier–Stokes hydrodynamics when
the gradients are strongly correlated to the collisional cooling [29].
From a microscopic point of view, the simple shear flow problem becomes
spatially uniform in the local Lagrangian frame moving with the flow velocity
u. In this frame [28, 75, 76], the velocity distribution functions adopt the
form: fi(r,v) → fi(V), where Vk = vk − akℓrℓ is the peculiar velocity. Here,
akℓ = aδkxδℓy. Under these conditions, the set of Boltzmann kinetic equations
(with Fi = 0) for an isolated system reads
− aVy
fi(V) =
Jij [V|fi, fj] , (i = 1, 2). (123)
The most relevant transport properties in a shear flow problem are obtained
from the pressure tensor P = P1 + P2, where Pi is the partial pressure tensor
of the species i given by
Pi,kℓ = mi
dVVkVℓfi(V). (124)
The trace of Pi defines the partial temperatures Ti as Ti = TrPi/dni.
As said before, these temperatures measure the mean kinetic energy of each
species. The elements of the pressure tensor Pi can be obtained by multiplying
the Boltzmann equation (123) by miVV and integrating over V. The result
akmPi,mℓ + aℓmPi,mk =
Aij,kℓ, (125)
where we have introduced the collisional moments Aij as
38 Vicente Garzó
Aij,kℓ = mi
dVVkVℓJij [V|fi, fj ]. (126)
From Eq. (125), in particular, one gets the balance equation for the partial
temperature Ti
aPi,xy = −
piζi, (127)
where pi = niTi is the partial pressure of species i and ζi is defined by Eq.
(10). According to Eq. (127), the (steady) partial temperature in the simple
shear flow problem can be obtained by equating the viscous heating term
a|Pi,xy| to the collisional cooling term (d/2)piζi.
The determination of Aij requires the knowledge of the velocity distribu-
tion functions fi. This is quite a formidable task, even in the monocomponent
case [27]. However, as in the elastic case, one expects to get a good estimate
of Aij by using Grad’s approximation [43]:
f(V) → fi,M (V)
Ci,kℓVkVℓ
, (128)
where fi,M is a Maxwellian distribution at the temperature of the species i,
i.e.,
fi,M (V) = ni
. (129)
As happens in the case of homogeneous states, in general the three temper-
atures T , T1, and T2 are different in the inelastic case. For this reason we
choose the parameters in the Maxwellians so that it is normalized to ni and
provides the exact second moment of fi. The Maxwellians fi,M for the two
species can be quite different due to the temperature differences. This aspect
is essential in a two-temperature theory and has not been taken into account
in most of the previous studies [14, 71, 72]. The coefficient Ci can be identified
by requiring the moments with respect to VV of the trial function (128) to
be the same as those for the exact distribution fi. This leads to
− I (130)
With this approximation, the Boltzmann collisional moments Aij can be
explicitly evaluated. The result is [77, 78]
Aij = −
2π(d−1)/2
dΓ (d/2)
mininjµjiσ
(1 + αij)
Tj − Ti
(mj/mi)Ti + Tj
1− αij
1 + (miTj/mjTi)
Ci − Cj
1 + (mjTi/miTj)
2(d+ 2)
, (131)
Kinetic Theory for Binary Granular Mixtures at Low-Density 39
where
λij = 2µji
Tj − Ti
(mj/mi)Ti + Tj
(2d+ 3− 3αij). (132)
The partial cooling rates ζi can be easily obtained from Eqs. (10) and (131).
1.0 1.5 2.0 2.5 3.0 3.5 4.0
α=0.8
α=0.9
Fig. 11. Plot of the temperature ratio T1/T2 as a function of the size ratio σ1/σ2 =
(m1/m2)
1/2 for a two-dimensional system in the case x1 = 1/2 and two different
values of the (common) coefficient of restitution: α = 0.9 and α = 0.8. The solid
lines are the theoretical predictions based on Grad’s solution, while the symbols
refer to the DSMC results. The dashed lines correspond to the results obtained from
the stochastic thermostat condition (31).
Substitution of Eq. (131) into the set of equations (125) allows one to get
the partial pressure tensor Pi in terms of the temperature ratio γ = T1/T2
and the parameters of the mixture. The temperature ratio can be obtained
from Eq. (127) as
x2ζ2P1,xy
x1ζ1P2,xy
. (133)
When the expressions of Pi and ζi are used in Eq. (133), one gets a closed
equation for the temperature ratio γ, that can be solved numerically. In Fig.
11 we plot γ versus the diameter ratio σ1/σ2 for a two-dimensional (d =
2) granular gas with x1 = 1/2 and two different values of α. The symbols
refer to the simulation data obtained from the DSMC method [79]. Here,
we have assumed that the disks are made of the same material, and hence
αij = α and m1/m2 = (σ1/σ2)
2. The dependence of γ on σ1/σ2 obtained in
the homogeneous steady state driven by the stochastic thermostat (30) is also
included for comparison. It is clearly seen that the kinetic theory results based
40 Vicente Garzó
on Grad’s solution agree very well with simulation data, even for quite large
values of the size ratio. In addition, the thermostat results overestimate the
simulation ones (especially for large mass ratio), showing that the properties
of the system are not insensitive to the way at which the granular gas is driven.
� �� ���
Fig. 12. Plot of the reduced pressure p∗ versus the mass ratio µ = m1/m2 for a two-
dimensional system with σ1 = σ2, x1 = 1/2 and α = 0.9. The solid line corresponds
to the theoretical predictions derived from Grad’s solution, the dotted line refers
to the latter theory but using the expression of T1/T2 obtained from the stochastic
thermostat condition (31), and the dashed line is the result obtained from Grad’s
solution by assuming the equality of the partial temperatures (γ = 1). The symbols
are the DSMC results.
Let us now consider the transport coefficients. To analyze the rheologi-
cal properties in the steady state, it is convenient to introduce dimensionless
quantities. As usual [72], for a low-density gas we introduce the reduced pres-
sure p∗ and the reduced shear viscosity η∗ as
, (134)
, (135)
where η = −Pxy/a is the non-Newtonian shear viscosity, Pxy = P1,xy + P2,xy
and ν = [π(d−1)/2/Γ (d/2)]nσd−112 v0. In Figs. 12 and 13, we plot p
∗ and η∗,
respectively, as functions of the mass ratio µ = m1/m2 for an equal-size
(σ1 = σ2) binary mixture of disks (d = 2) with x1 = 1/2 and α = 0.9.
Kinetic Theory for Binary Granular Mixtures at Low-Density 41
� �� ���
Fig. 13. Plot of the reduced shear viscosity η∗ versus the mass ratio µ = m1/m2
for a two-dimensional system with σ1 = σ2, x1 = 1/2 and α = 0.9. The solid line
corresponds to the theoretical predictions derived from Grad’s solution, the dotted
line refers to the latter theory but using the expression of T1/T2 obtained from the
stochastic thermostat condition (31), and the dashed line is the result obtained from
Grad’s solution by assuming the equality of the partial temperatures (γ = 1). The
symbols are the DSMC results.
We have also included the predictions for p∗ and η∗ given by the kinetic
theory but taking the expression of γ derived when the system is driven by
the stochastic thermostat [40]. We observe again in both figures an excellent
agreement between the Boltzmann theory based on Grad’s solution and the
DSMC results, even for very disparate values of the mass ratio. With respect
to the influence of energy nonequipartition, Fig. 12 shows that p∗ presents a
non-monotonic behaviour with the mass ratio whereas the theoretical predic-
tions with the equipartition assumption monotonically increase with µ. In the
case of the shear viscosity, as seen in Fig. 13, both theories (with and with-
out energy nonequipartition) predict a non-monotonic dependence of η∗ on
µ. However, at a quantitative level, the influence of energy nonequipartition
is quite significant over the whole range of mass ratios considered. The non-
monotonic dependence of p∗ and η∗ on µ obtained here from the Boltzmann
kinetic theory also agrees qualitatively well with MD simulations carried out
for bidisperse dense systems [72]. Thus, for instance, the minimum values of p∗
and η∗ are located close to µ = 10 in both dilute and dense cases. Moreover,
the predictions for the transport properties given from the present theory
by taking the stochastic thermostat expression of γ are quite close to those
obtained from the actual value of γ, especially for large mass ratios.
42 Vicente Garzó
11 Summary and Concluding Remarks
The primary objective of this review has been to derive the Navier–Stokes hy-
drodynamic equations of a binary mixture of granular gases from the (inelas-
tic) Boltzmann kinetic theory. The Chapman–Enskog method [4, 43] is used
to solve the Boltzmann equation up to the first order in the spatial gradients
and the associated transport coefficients are given in terms of the solutions of
a set of linear integral equations. These equations have been approximately
solved by taking the leading terms in a Sonine polynomial expansion. Compar-
ison with controlled numerical simulations in some idealized conditions shows
quite a good agreement between theory and simulation even for strong dis-
sipation. This supports the idea that the hydrodynamic description (derived
from kinetic theory) appears to be a powerful tool for analysis and predictions
of rapid flow gas dynamics of polydisperse systems [12].
The reference state in the Chapman–Enskog expansion has been taken
to be an exact solution of the uniform Boltzmann equation. An interesting
and important result of this solution [15] is that the partial temperatures
(which measure the mean kinetic energy of each species) are different. This
does not mean that there are additional degrees of freedom since the partial
temperatures can be expressed in terms of the global temperature. This is
confirmed by noting that Haff’s cooling law [1] (in the free cooling case) is the
hydrodynamic mode at long wavelengths and MD simulations confirm that
the global temperature dominates after a transient period of a few collision
times [34]. In this case, only the global temperature should appear among
the hydrodynamic fields. Nevertheless, the species temperatures play a new
and interesting secondary role [20]. For an ordinary (molecular) gas, there is
a rapid velocity relaxation in each fluid cell to a local equilibrium state on
the time scale of a few collisions (e.g., as illustrated by the approach to Haff’s
law). Subsequently, the equilibration among cells occurs via the hydrodynamic
equations. In each cell the species velocity distributions are characterized by
the species temperatures. These are approximately the same due to equiparti-
tion, and the hydrodynamic relaxation occurs for the single common temper-
ature [43]. A similar rapid velocity relaxation occurs for granular gases in each
small cell, but to a universal state different from local equilibrium and one for
which equipartition no longer occurs. Hence, the species temperatures Ti are
different from each other and from the overall temperature T of the cell. Nev-
ertheless, the time dependence of all temperatures (in the free cooling case)
is the same in this and subsequent states, i.e., they are proportional to the
global temperature. This implies that the species temperatures do not provide
any new dynamical degree of freedom at the hydrodynamic stage. However,
they still characterize the shape of the partial velocity distributions and affect
the quantitative averages calculated with these distributions. The transport
coefficients for granular mixtures therefore have new quantitative effects aris-
ing from the time independent temperature ratios for each species [19]. This
view contrasts with some recent works [80], where additional equations for
Kinetic Theory for Binary Granular Mixtures at Low-Density 43
each species temperature have been included among the hydrodynamic set.
However, as mentioned before, this is an unnecessary complication, describing
additional kinetics beyond hydrodynamics that is relevant only on the time
scale of a few collisions.
Another important issue discussed here has been the applicability of the
Navier–Stokes transport coefficients since their expressions are not restricted
to weak inelasticity [12]. However, the Navier–Stokes hydrodynamic equations
themselves may or may not be limited with respect to inelasticity, depending
on the particular states studied. The Chapman–Enskog method assumes that
the relative changes of the hydrodynamic fields over distances of the order
of the mean free path are small. In the case of ordinary fluids this can be
controlled by the initial or boundary conditions. For granular gases the situ-
ation is more complicated since in some cases (e.g., steady states such as the
simple shear flow problem [29]) the boundary conditions imply a relationship
between dissipation and gradients so that both cannot be chosen indepen-
dently. In these cases, the Navier–Stokes approximation only holds for nearly
elastic particles. However, the transport coefficients characterizing the Navier–
Stokes hydrodynamic equations are nonlinear functions of the coefficients of
restitution, regardless the applicability of those equations.
In spite of the above cautions, the Navier–Stokes approximation is ap-
propriate and accurate for a wide class of flows. One group refers to spatial
perturbations of the homogeneous cooling state (HCS) for an isolated sys-
tem. Both MD and DSMC simulations [7] have confirmed the dependence
of the Navier–Stokes transport coefficients on the coefficient of restitution,
and application of the Navier–Stokes hydrodynamics with these coefficients
to describe cluster formation has also been confirmed quantitatively [60]. The
same kinetic theory results apply to driven systems as well. This is so since
the reference state is a local HCS whose parameters vary throughout the sys-
tem to match the physical values in each cell. Examples include application of
Navier–Stokes hydrodynamics from kinetic theory to symmetry breaking and
density/temperature profiles in vertical vibrated gases, for comparison with
simulation [81]. Similar comparison with Navier–Stokes hydrodynamics of the
latter and of supersonic flow past a wedge in real experiments has been given
[82, 83], showing both qualitative and quantitative agreement. In summary,
the Navier–Stokes equations with the constitutive equations presented here
remain an important and useful description for a wide class of granular flow,
although more limited than for normal gases.
The explicit knowledge of the transport coefficients and the cooling rate
allows one to make some applications of the Navier–Stokes hydrodynamic
equations. One of them has been to obtain the linear hydrodynamic equa-
tions for small perturbations of the homogenous cooling state. The resulting
equations exhibit a long wavelength instability for three of the modes. This
is quite similar to the case of a monocomponent granular gas [5, 54, 57], and
in fact the same modes are unstable here. The additional diffusion mode for
two species behaves as for a normal fluid.
44 Vicente Garzó
On the other hand, the constitutive equations for the mass and heat fluxes
of a granular binary mixture differ from those obtained for ordinary fluids [4].
This is because the usual restrictions of irreversible thermodynamics no longer
apply. These restrictions include Onsager’s reciprocal relations among various
transport coefficients and the extent to which these are violated has also been
shown here. Another application of the Navier–Stokes equations has been to
assess the violation of the Einstein relation between the diffusion and mobility
coefficients. In the undriven case, the analysis shows that this violation is
due to three independent reasons [23]: the absence of the Gibbs state, the
cooling of the reference state, and the occurrence of different temperatures for
the particle and surrounding fluid. However, when the mixture is subjected
to stochastic driving, a modified Einstein relation suggested by recent MD
simulations [52] has also been analyzed. In this case, the results show that
the deviations of the (modified) Einstein ratio from unity are in general very
small (less than 1%), in agreement with MD simulations [52].
Thermal diffusion becomes the relevant segregation mechanism in agitated
granular mixtures at large shaking amplitudes. In these conditions, the use of
the Boltzmann kinetic theory for low-density gases appears justified to under-
stand the influence of thermal gradient on segregation phenomena. The ther-
mal diffusion factor in a heated granular mixture has been explicitly evaluated
from the Chapman–Enskog solution to the Boltzmann equation. The results
show that the criterion for the transition Brazil-nut effect ⇐⇒reverse Brazil-
nut effect is provided by the control parameter θ = m2T1/m1T2 [25, 42, 67].
Given that the energy equipartition is broken, the condition θ = 1 is quite
complex since it involves all the parameters of the system: composition,
masses, sizes, and coefficients of restitution. The Boltzmann kinetic theory
results agree qualitatively well with recent MD simulations [41] within the
range of parameter space analyzed.
The hydrodynamic description also seems to be justified in the case of
steady states that are inherently beyond the scope of the Navier–Stokes hy-
drodynamic equations. The reason for this non-Newtonian behavior is the
existence of an internal mechanism, collisional cooling, that sets the scale of
the spatial gradients in the steady state. For ordinary fluids, this scale can be
externally controlled by external sources so that the conditions for Navier–
Stokes hydrodynamics apply. On the other hand, for granular gases, collisional
cooling is fixed by the mechanical properties of the particles of the system and
so the gas can depart from the Navier–Stokes description. One well-known ex-
ample of steady states is the simple or uniform shear flow (USF). However, in
spite of the extensive prior work on USF for granular fluids [27, 72, 73], the in-
herent non-Newtonian character of this state has not been conveniently taken
into account. In fact, MD simulations of steady USF have been used for gran-
ular fluids to measure the Newtonian or Navier–Stokes shear viscosity. The
results derived here from Grad’s solution and DSMC simulations show that
USF is an ideal testing ground for the study of rheology since any choice of the
shear rate and the coefficients of restitution αij will provide non-Newtonian
Kinetic Theory for Binary Granular Mixtures at Low-Density 45
effects. It is one of the fascinating features of granular fluids that phenom-
ena associated with complex fluids are more easily accessible than for simple
atomic fluids [12, 84].
Hydrodynamics derived from hard-sphere models have found widespread
use in the description of numerous industrial processes involving solid parti-
cles. Of particular relevance are high-speed, gas-solid flows as found in pneu-
matic conveyors (of ores, chemicals, grains, etc.) and fluidized beds (for fluid
catalytic cracking, power generation, granulation of pharmaceutical powders,
synthesis of fine chemicals like titania, etc.). Such descriptions are now stan-
dard features of commercial and research codes. Those codes rely upon accu-
rate transport properties and a first order objective is to assure this accuracy
from a careful theoretical treatment. As shown in this review, the price of this
approach, in contrast to more phenomenological approaches, is an increasing
complexity of the expressions as the systems become more complex.
The analysis carried out in this presentation has been focused on mix-
tures in the dilute regime, where the collisional transfer contributions to the
transport coefficients are neglected and only their kinetic contributions are
considered. A further step is to develop a theory for moderately dense gran-
ular mixtures. This will provide a fundamental basis for the application of
hydrodynamics under realistic conditions. Possible extension of the present
Boltzmann kinetic theory to higher densities can be done in the context of
the Enskog kinetic equation [43]. Preliminary results [85] have been restricted
to the uniform shear flow state to get directly the shear viscosity coefficient of
a heated granular mixture. The extension of this study [85] to states with gra-
dients of concentration, pressure, and temperature is somewhat intricate due
to subtleties associated with the spatial dependence of the pair correlations
functions considered in the revised Enskog theory. A future work is to ex-
tend the results derived for moderately dense mixtures of smooth elastic hard
spheres [86] to inelastic collisions. This would allow us to assess the influence
of density on the different problems addressed in this review. Of course, the
precise expressions for transport coefficients in this case will be even more
complex than for a dilute gas due to the expanded parameter space. However,
this complexity is not a problem for implementation in a code.
As shown along this overview, granular mixtures exhibit a wide range of
interesting phenomena for which the Navier–Stokes hydrodynamic equations
can be considered as an accurate and practical tool. However, due to their
complexity, many of their features are not fully understood. Kinetic theory
and hydrodynamics (in the broader sense) can be expected to provide some
insight into the understanding of such complex materials.
Acknowledgments
I want to acknowledge J. W. Dufty, J. M. Montanero, and A. Santos in
their roles as collaborators and critics for much of the material discussed here.
46 Vicente Garzó
Partial support of the Ministerio de Ciencia y Tecnoloǵıa (Spain) through
Grant No. FIS2004-01399 (partially financed by FEDER funds) and from the
European Community’s Human Potential Programme HPRN-CT-2002-00307
(DYGLAGEMEM) is also acknowledged.
References
1. P. K. Haff, J. Fluid Mech. 134, 187 (1983).
2. L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999).
3. See for instance, J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids
(Academic Press, London, 1986).
4. S. Chapman and T. G. Cowling, The Mathematical Theory of Nonuniform Gases
(Cambridge University Press, Cambridge, 1970).
5. J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998).
6. J. J. Brey and D. Cubero, in Granular Gases, edited by T. Pöschel and S. Luding
(Lectures Notes in Physics, Springer, Vol. 564, 2001), pp. 59–78.
7. J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Europhys. Lett. 48, 359 (1999).
8. J. J. Brey, M. J. Ruiz-Montero, D. Cubero, and R. Garćıa-Rojo, Phys. Fluids
12, 876 (2000).
9. J. J. Brey and M. J. Ruiz-Montero, Phys. Rev. E 70, 051301 (2004).
10. J. M. Montanero, A. Santos, and V. Garzó, in Rarefied Gas Dynamics 24, edited
by M. Capitelli (American Institute of Physics, Vol. 762, 2005), pp. 797–802;
preprint cond-mat/0411219.
11. J. W. Dufty, J. Phys.: Condens. Matt. 12, A47 (2000).
12. J. W. Dufty, in Recent Research Develop. in Stat. Phys. 2, 21 (2002) and preprint
cond-mat/0108444.
13. I. Goldhirsch, Annu. Rev. Fluid Mech. 35, 267 (2003).
14. J. T. Jenkins and F. Mancini, Phys. Fluids A 1, 2050 (1989); P. Zamankhan,
Phys. Rev. E 52, 4877 (1995); B. Arnarson and J. T. Willits, Phys. Fluids 10,
1324 (1998); J. T. Willits and B. Arnarson, ibid. 11, 3116 (1999); M. Alam, J.
T. Willits, B. Arnarson, and S. Luding, ibid. 14, 4085 (2002); B. Arnarson and
J. T. Jenkins, ibid. 16, 4543 (2004); D. Serero, I. Goldhirsch, S. H. Noskowicz,
and M.-L. Tan, J. Fluid Mech. 554, 237 (2006).
15. V. Garzó and J. W. Dufty, Phys. Rev. E 60, 5706 (1999).
16. See for instance, J. M. Montanero and V. Garzó, Gran. Matt. 4, 17 (2002); A.
Barrat and E. Trizac, ibid. 4, 57 (2002); S. R. Dahl, C. M. Hrenya, V. Garzó, and
J. W. Dufty, Phys. Rev. E 66, 041301 (2002); R. Pagnani, U. M. B. Marconi,
and A. Puglisi, Gran. Matt. 66, 051304 (2002); D. Paolotti, C. Cattuto, U. M.
B. Marconi, and A. Puglisi, ibid. 5, 75 (2003); P. Krouskop and J. Talbot, Phys.
Rev. E 68, 021304 (2003); H. Wang, G. Jin, and Y. Ma, ibid. 68, 031301 (2003);
J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, Phys. Rev. Lett. 95, 098001
(2005); Phys. Rev. E 73, 031301 (2006); M. Schröter, S. Ulrich, J. Kreft, J. B.
Swift, and H. L. Swinney, Phys. Rev. E 74, 011307 (2006).
17. R. D. Wildman and D. J. Parker, Phys. Rev. Lett. 88, 064301 (2002); K. Feitosa
and N. Menon, ibid. 88, 198301 (2002); M. Schröter, S. Ulrich, J. Kreft, J. B.
Swift, and H. L. Swinney, Phys. Rev. E 74, 011307 (2006).
18. J. Jenkins and F. Mancini, J. Appl. Mech. 54, 27 (1987).
http://arxiv.org/abs/cond-mat/0411219
http://arxiv.org/abs/cond-mat/0108444
Kinetic Theory for Binary Granular Mixtures at Low-Density 47
19. V. Garzó and J. W. Dufty, Phys. Fluids 14, 1476 (2002).
20. V. Garzó, J. M. Montanero, and J. W. Dufty, Phys. Fluids 18, 083305 (2006).
21. J. M. Montanero and V. Garzó, Phys. Rev. E 67, 021308 (2003).
22. V. Garzó and J. M. Montanero, Phys. Rev. E 69, 021301 (2004).
23. J. W. Dufty and V. Garzó, J. Stat. Phys. 105, 723 (2001).
24. V. Garzó, Physica A 343, 105 (2004).
25. V. Garzó, Europhys. Lett. 75, 521 (2006).
26. C. S. Campbell, Annu. Rev. Fluid Mech. 22, 57 (1990).
27. In the case of a monocomponent system, see for instance, C. K. K. Lun, S. B.
Savage, D. J. Jeffrey, and N. Chepurniy, J. Fluid Mech. 140, 223 (1984); J. T.
Jenkins and M. W. Richman, ibid. 192, 313 (1988); C. S. Campbell, ibid. 203,
449 (1989); M. A. Hopkins and H. H. Shen, ibid. 244, 477 (1992); P. J. Schmid
and H. K. Kytömaa, ibid. 264, 255 (1994); C. K. K. Lun and A. A. Bent, ibid.
258, 335 (1994); I. Goldhirsch and M. L. Tan, Phys. Fluids 8, 1752 (1996); N.
Sela, I. Goldhirsch, and S. H. Noskowicz, ibid. 8, 2337 (1997); J. J. Brey, M. J.
Ruiz-Montero, and F. Moreno, Phys. Rev. E 55, 2846 (1997); C.-S. Chou and
M. W. Richman, Physica A 259, 430 (1998); C.-S. Chou ibid. 287, 127 (2000);
290, 341 (2001); J. M. Montanero, V. Garzó, A. Santos, and J. J. Brey, J. Fluid
Mech. 389, 391 (1999); A. Astillero and A. Santos, Phys. Rev. E 72, 031309
(2005).
28. V. Garzó and A. Santos, Kinetic Theory of Gases in Shear Flows. Nonlinear
Transport (Kluwer Academic, Dordrecht, 2003).
29. A. Santos, V. Garzó, and J. W. Dufty, Phys. Rev. E 69, 061303 (2004).
30. A. Goldshtein and M. Shapiro, J. Fluid Mech. 282, 75 (1995).
31. J. J. Brey, J. W. Dufty, and A. Santos, J. Stat. Phys. 87, 1051 (1997).
32. T. P. C. van Noije and M. H. Ernst, Gran. Matt. 1, 57 (1998).
33. J. M. Montanero and V. Garzó, Gran. Matt. 4, 17 (2002).
34. S. R. Dahl, C. M. Hrenya, V. Garzó, and J. W. Dufty, Phys. Rev. E 66, 041301
(2002).
35. D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids
(Academic Press, London, 1990).
36. W. G. Hoover, Computational Statistical Mechanics (Elsevier, Amsterdam,
1991).
37. J. M. Montanero and A. Santos, Gran. Matt. 2, 53 (2000).
38. D. R. M. Williams and F. C. McKintosh, Phys. Rev. E 54, R9 (1996).
39. C. Henrique, G. Batrouni, and D. Bideau, Phys. Rev. E 63, 011304 (2000).
40. A. Barrat and E. Trizac, Gran. Matt. 4, 52 (2002).
41. M. Schröter, S. Ulrich, J. Kreft, J. B. Swift, and H. L. Swinney, Phys. Rev. E
74, 011307 (2006).
42. J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, Phys. Rev. Lett. 95, 098001
(2005); Phys. Rev. E 73, 031301 (2006).
43. J. Ferziger and H. Kaper, Mathematical Theory of Transport Processes in Gases
(North-Holland, Amsterdam, 1972).
44. V. Garzó and J. M. Montanero, preprint cond-mat/0604552.
45. V. Garzó and J.W. Dufty, Phys. Rev. E 59, 5895 (1999).
46. J. F. Lutsko, Phys. Rev. E 72, 021306 (2005).
47. G. A. Bird, Molecular Gas Dynamics and the Direct Simulation Monte Carlo of
Gas Flows (Clarendon, Oxford, 1994).
48. J. A. McLennan, Introduction to Nonequilibrium Statistical Mechanics (Prentice-
Hall, New Jersey, 1989).
http://arxiv.org/abs/cond-mat/0604552
48 Vicente Garzó
49. A. Santos and J. W. Dufty, Phys. Rev. Lett. 97, 058001 (2006).
50. M. López de Haro and E. G. D. Cohen, J. Chem. Phys. 80, 408 (1984).
51. V. Garzó and J. M. Montanero, Physica A 313, 336 (2002).
52. A. Barrat, V. Loreto, and A. Puglisi, Physica A 334, 513 (2004).
53. S. R. de Groot and P. Mazur, Nonequilibrium Thermodynamics (Dover, New
York, 1984).
54. N. Brilliantov and T. Pöschel, Kinetic Theory of Granular Gases (Oxford Uni-
versity Press, Oxford, 2004).
55. I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993); I. Goldhirsch,
M.L. Tan, and G. Zanetti, J. Sci. Comput. 8, 1 (1993).
56. S. Luding and H.J. Herrmann, Chaos 9, 673 (1999).
57. V. Garzó, Phys. Rev. E 72, 021106 (2005).
58. P. Résibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley,
NY, 1977).
59. J. P. Boon and S. Yip, Molecular Hydrodynamics (Dover, NY, 1980).
60. J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 60, 3150 (1999).
61. See, for instance, A. Kudrolli, Rep. Prog. Phys. 67, 209 (2004).
62. J. M. Kincaid, E. G. D. Cohen, and M. López de Haro, J. Chem. Phys. 86, 963
(1983).
63. D. C. Hong, P. V. Quinn, and S. Luding, Phys. Rev. Lett. 86, 3423 (2001).
64. J. T. Jenkins and D. K. Yoon, Phys. Rev. Lett. 88, 194301 (2002); D. K. Yoon
and J. T. Jenkins, Phys. Fluids 18, 073303 (2006).
65. A. P. J. Breu, H. M. Ensner, C. A. Kruelle, and I. Rehberg, Phys. Rev. Lett.
90, 014302 (2003).
66. T. Schautz, R. Brito, C. A. Kruelle, and I. Rehberg, Phys. Rev. Lett. 95, 028001
(2005) and references therein.
67. L. Trujillo, M. Alam, and H. J. Herrmann, Europhys. Lett. 64, 190 (2003).
68. D. Serero, I. Goldhirsch, S. H. Noskowicz, and M.-L. Tan, J. Fluid Mech. 554,
237 (2006).
69. M. L. Tan and I. Goldhirsch, Phys. Rev. Lett. 81, 3022 (1998).
70. J. W. Dufty and J. J. Brey, Phys. Rev. Lett. 82, 4566 (1999).
71. J. T. Willits and B. Arnarson, Phys. Fluids 11, 3116 (1999); M. Alam, J. Willits,
B. Arnarson, and S. Luding, ibid. 14, 4085 (2002).
72. M. Alam and S. Luding, Gran. Matt. 4, 139 (2002); J. Fluid Mech. 476, 69
(2003); Phys. Fluids 17, 063303 (2005).
73. M. Alam and C. M. Hrenya, Phys. Rev. E 63, 0613018 (2001); R. Clelland
and C. M. Hrenya, ibid. 65, 031301 (2002); S. R. Dahl, R. Clelland, and C.
M. Hrenya, Phys. Fluids 14, 1972 (2002); S. R. Dahl, R. Clelland, and C. M.
Hrenya, Powder Tech., 138, 7 (2003); S. R. Dahl and C. M. Hrenya, Phys. Fluids
16, 1 (2004); H. Iddir, H. Arastoopour, and C. M. Hrenya, Powder Tech. 151,
117 (2005).
74. See for instance, Nonlinear Fluid Behavior, edited by H. J. M. Hanley (North-
Holland, Amsterdam, 1983).
75. A. W. Lees and S. F. Edwards, J. Phys. C 5, 1921 (1972).
76. J. W. Dufty, A. Santos, J. J. Brey, and R. F. Rodŕıguez, Phys. Rev. A 33, 459
(1986).
77. J. M. Montanero and V. Garzó, Physica A 310, 17 (2002).
78. V. Garzó, Phys. Rev. E 66, 021308 (2002).
79. V. Garzó and J. M. Montanero, Gran. Matt. 5, 165 (2003).
Kinetic Theory for Binary Granular Mixtures at Low-Density 49
80. L. Huilin, D. Gidaspow, and E. Manger, Phys. Rev. E 64, 061301 (2001); M.F.
Ramahan, J. Naser, and P. J. Witt, Powder Tech. 138, 82 (2003); J. E. Galvin,
S. R. Dahl, and C. M. Hrenya, J. Fluid Mech. 528, 207 (2005).
81. J. J. Brey, M. J. Ruiz-Montero, F. Moreno, and R. Garćıa-Rojo, Phys. Rev. E
63, 061305 (2001); ibid. 65, 061302 (2002)
82. X. Yang, C. Huan, D. Candela, R. W. Mair, and R. L. Walsworth, Phys. Rev.
Lett. 88, 044301 (2002); C. Huan, X. Yang, D. Candela, R. W. Mair, and R. L.
Walsworth, Phys. Rev. E 69, 041302 (2004).
83. C. Bizon, M. D. Shattuck, J. B. Swift, and Harry L. Swinney, Phys. Rev. E 60,
4340 (1999); E. C. Rericha, C. Bizon, M. D. Shattuck, and H. L. Swinney, Phys.
Rev. Lett. 88, 014302 (2002).
84. I. Goldhirsch, in Granular Gases, edited by T. Pöschel and S. Luding (Lectures
Notes in Physics, Springer, Vol. 564, 2001), pp. 79–99.
85. V. Garzó and J. M. Montanero, Phys. Rev. E 68, 041302 (2003).
86. M. López de Haro, E. G. D. Cohen, and J. M. Kincaid, J. Chem. Phys. 78, 2746
(1983).
	Kinetic Theory for Binary Granular Mixtures at Low-Density
	Vicente Garzó

  Many features of granular media can be modelled as a fluid of hard spheres
with {\em inelastic} collisions. Under rapid flow conditions, the macroscopic
behavior of grains can be described through hydrodynamic equations. At
low-density, a fundamental basis for the derivation of the hydrodynamic
equations and explicit expressions for the transport coefficients appearing in
them is provided by the Boltzmann kinetic theory conveniently modified to
account for inelastic binary collisions. The goal of this chapter is to give an
overview of the recent advances made for binary granular gases by using kinetic
theory tools. Some of the results presented here cover aspects such as
transport properties, energy nonequipartition, instabilities, segregation or
mixing, non-Newtonian behavior, .... In addition, comparison of the analytical
results with those obtained from Monte Carlo and molecular dynamics simulations
is also carried out, showing the reliability of kinetic theory to describe
granular flows even for strong dissipation.

Introduction
Granular systems have attracted the attention of the physics community in
the past few years, in part because the behavior of these systems under many
conditions exhibit a great similarity to ordinary fluids [1]. These conditions
include rapid, dilute flows where the dominant transfer of momentum and
energy is through binary collisions of the grains. The main difference from
ordinary fluids is the absence of energy conservation, leading to both obvious
and subtle modifications of the usual macroscopic balance equations as well as
the constitutive equations for the irreversible fluxes. However, in spite of the
utility of hydrodynamics to describe rapid granular flows, there are still some
open questions about its domain of validity and the associated constitutive
equations appearing in the hydrodynamic equations [2].
To isolate the effects of such collisional dissipation from other important
properties of granular media, an idealized microscopic model system is usually
http://arxiv.org/abs/0704.1211v1
2 Vicente Garzó
considered: a system composed by smooth hard spheres with inelastic colli-
sions. As in the elastic case, the collisions are specified in terms of the change
in relative velocity at contact but with a decrease in the magnitude of the nor-
mal component measured by a positive coefficient of restitution α ≤ 1. This
parameter distinguishes the granular fluid (α < 1) from the ordinary fluid
(α = 1). Given that the hard sphere system with elastic collisions has been
widely studied for both equilibrium and non-equilibrium statistical mechanics
[3], it is tempting to apply the same methods for the case of inelastic collisions.
However, some care is warranted in translating properties of ordinary fluids to
granular fluids. In this presentation, a kinetic theory description based on the
Boltzmann kinetic equation (which applies at sufficiently low density) will be
considered as the appropriate tool to study granular flows from a microscopic
point of view.
Although many efforts have been devoted in the past few years to the
understanding of granular fluids, the derivation of the form of the constitu-
tive equations with explicit expressions for the transport coefficients is still a
subject of interest and controversy. The conditions to obtain a hydrodynamic
description are expected to be similar to those for normal fluids. For a given
initial state there are two stages of evolution. First, during the kinetic stage
there is rapid velocity relaxation to a “universal” velocity distribution that
depends on the average local density, temperature, and flow velocity. Subse-
quently, the hydrodynamic stage is described through a slower evolution of
these local hydrodynamic fields as they approach uniformity. The solution to
the Boltzmann equation in this second stage is said to be normal, where all
space and time dependence of the distribution function occurs through the
macroscopic hydrodynamic fields. The Chapman–Enskog method [4] provides
a constructive means to obtain an approximation to such a solution for states
whose spatial gradients are not too large. In this case, the explicit form of this
normal solution is given as a perturbation expansion in the spatial gradients of
the fields. This solution is then used to evaluate the fluxes in the macroscopic
balance equations in terms of these gradients. To lowest order the balance
equations become the granular Euler equations while to second order they are
the granular Navier–Stokes equations. In carrying out this analysis, explicit
forms for the transport coefficients are obtained as functions of the coefficient
of restitution and other parameters of the collision operator. In this general
context, the study of hydrodynamics for granular gases is the same as that
for ordinary fluids.
The derivation of hydrodynamics from the inelastic Boltzmann equation
has been widely covered in the case of a monodisperse gas where the particles
are of the same mass and size. As for elastic collisions, the transport coeffi-
cients are given in terms of the solutions of linear integral equations [5, 6],
which are approximately solved by using Sonine polynomial expansions. The
estimates for the transport coefficients provided by the Sonine solution com-
pare in general quite well with both direct Monte Carlo simulation (DSMC)
of the Boltzmann equation and molecular dynamics (MD) simulation of the
Kinetic Theory for Binary Granular Mixtures at Low-Density 3
gas, even for relatively strong degrees of dissipation [7, 8, 9, 10]. This good
agreement supports the formal theoretical analysis and the claim that hydro-
dynamics is not limited to nearly elastic particles [11, 12, 13].
Nevertheless, a real granular system is generally characterized by some
degrees of polydispersity in density and size, which leads to phenomena very
often observed in nature and experiments, such as separation or segregation.
Needless to say, the study of granular mixtures is much more complicated
than for a monodisperse gas since not only the number of transport coeffi-
cients in a multicomponent system is higher than that of a single gas, but
also they depend on parameters such as masses, sizes, composition as well
as several independent coefficients of restitution αij . Due to these difficulties,
studies for multicomponent gases are more scarce in the literature. Many of
the previous attempts [14] to derive hydrodynamics from kinetic theory were
carried out in the quasi-elastic limit where the equipartition of energy can be
considered as an acceptable assumption. In addition, according to this level of
approximation, the inelasticity is only accounted for by the presence of a sink
term in the energy balance equation, so that the expressions for the transport
coefficients are the same as those obtained for ordinary fluids. However, the
theoretical prediction of the failure of energy equipartition in multicomponent
granular gases [15] has been confirmed by computer simulations [16], and even
observed in real experiments [17].
Although the possibility of nonequipartition was already pointed out many
years ago [18], it has not been until recently that a systematic study of the
effect of nonequipartition on the Navier–Stokes hydrodynamic equations has
been carefuly carried out [19, 20]. These new equations and associated trans-
port coefficients provide a somewhat more stringent test of the analysis since
the parameter space is much larger. As in the monodisperse case, explicit ex-
pressions for the transport coefficients requires also to consider Sonine poly-
nomial expansions. The numerical accuracy of this Sonine expansion has been
confirmed by comparison with Monte Carlo simulations of the Boltzmann
equation in the cases of the shear viscosity [21] and the tracer diffusion [22]
coefficients. Exceptions to this agreement are extreme mass or size ratios and
strong dissipation, although these discrepancies between theory and simula-
tion diminish as one considers more terms in the Sonine polynomial approxi-
mation [22].
The explicit knowledge of the Navier-Stokes transport coefficients allows
quantitative application of the nonlinear hydrodynamic equations to a num-
ber of interesting problems for granular mixtures, such as to quantify the
violation of the Einstein relation [23, 24] or the Onsager reciprocal relations
[20], the stability analysis of the homogeneous cooling state [20], and segrega-
tion induced by a thermal gradient [25]. In all the cases, the analysis clearly
shows the important role played by the inelasticity in the different physical
situations.
The analogy between rapid granular flow and ordinary fluids can be also
extended to many other transport situations. A particularly simple case, al-
4 Vicente Garzó
lowing detailed analysis even in far from equilibrium conditions is the simple
or uniform shear flow (USF) problem. Macroscopically, it is characterized by
uniform density and temperature and a constant mean velocity profile. This
is a well-known non-equilibrium problem widely studied, for both granular
monodisperse [13, 26, 27] and ordinary gases [28]. Nevertheless, the nature of
this state is quite different in each system. While for elastic fluids the temper-
ature increases monotonically in time due to viscous heating, a steady state
is possible for granular media when the effect of the viscous heating is exactly
compensated by the dissipation in collisions. Thus, in the steady state, there
is an intrinsic connection between the shear field and dissipation so that the
collisional cooling sets the strength of the velocity gradient. As a consequence,
the USF state is inherently non-Newtonian and the rheological properties of
the system cannot be obtained from the Navier-Stokes description, at least
for finite dissipation [29].
The aim of this chapter is to offer a short review of recent results ob-
tained for binary granular mixtures from the Boltzmann kinetic theory. It is
structured as follows. The Boltzmann kinetic equation for a granular binary
mixture and its associated macroscopic balance equations are introduced in
Sec. 2. Section 3 deals with the solutions to the Boltzmann equation for homo-
geneous states in the free cooling case as well as when the mixture is heated
by an external thermostat. The Chapman–Enskog method around the local
version of the homogenous distributions obtained in Sec. 2 is applied in Sec. 3
to get the form of the Navier–Stokes hydrodynamic equations. Theoretical re-
sults for the diffusion and shear viscosity transport coefficients are compared
with simulation data in Sec. 5, while the Einstein and the Onsager relations for
granular mixtures are analyzed in Secs. 6 and 7, respectively. The dispersion
relations for the hydrodynamic equations linearized about the homogeneous
cooling state are obtained in Sec. 8, showing that the homogeneous reference
state is unstable to long wavelength perturbations. The conditions for sta-
bility are identified as functions of the wave vector, the dissipation, and the
parameters of the mixture. Segregation due to thermal diffusion is studied in
Sec. 9 by using the Navier–Stokes description. A new criterion for segregation
is found that is consistent with recent experimental results. Section 10 deals
with the USF problem for a granular mixture. Finally, the paper is closed in
Sec. 11 with a discussion of the results presented here.
Before ending this section, I want to remark that the present account is a
personal perspective based on the author’s work and that of his collaborators
so that no attempt is made to include the extensive related work of many
others in this field. The references given are selective and apologies are of-
fered at the outset to the many other important contributions not recognized
explicitly.
Kinetic Theory for Binary Granular Mixtures at Low-Density 5
2 Boltzmann Kinetic Equation for Binary Mixtures of
Inelastic Hard Spheres
Consider a binary granular mixture composed by smooth inelastic disks
(d = 2) or spheres (d = 3) of massesm1 andm2, and diameters σ1 and σ2. The
inelasticity of collisions among all pairs is characterized by three independent
constant coefficients of restitution α11, α22, and α12 = α21, where αij ≤ 1 is
the coefficient of restitution for collisions between particles of species i and j.
Since the spheres are assumed to be perfectly smooth, only the translational
degrees of freedom of grains are affected by dissipation. In the low density
regime, a simultaneous interaction of more than two particles is highly unlike
and so can be neglected. Consequently, in a dilute gas the interactions among
the particles reduce to a succession of binary collisions. At this level of de-
scription, all the relevant information on the state of the system is contained
in the one-body velocity distribution functions fi(r,v; t) (i = 1, 2) defined so
that fi(r,v; t)drdv is the most probable (or average) number of particles of
species i which at time t lie in the volume element dr centered at the point
r and moving with velocities in the range dv about v. For an inelastic gas,
the distributions fi(r,v; t) (i = 1, 2) for the two species satisfy the coupled
nonlinear Boltzmann equations [30, 31]
∂t + v · ∇+
fi(r,v, t) =
Jij [v|fi(t), fj(t)] , (1)
where the Boltzmann collision operator Jij [v|fi, fj ] is
Jij [v1|fi, fj] = σ
dσ̂Θ(σ̂ · g12)(σ̂ · v12)
α−2ij fi(r,v
1, t)fj(r,v
2, t)− fi(r,v1, t)fj(r,v2, t)
. (2)
In Eq. (2), d is the dimensionality of the system, σij = (σi + σj) /2, σ̂ is an
unit vector along the line of centers, Θ is the Heaviside step function, and
v12 = v1 −v2 is the relative velocity. The primes on the velocities denote the
initial values {v′1,v
2} that lead to {v1,v2} following a binary (restituting)
collision:
1 = v1 − µji
1 + α−1ij
(σ̂ · v12)σ̂,
2 = v2 + µij
1 + α−1ij
(σ̂ · v12)σ̂, (3)
where µij ≡ mi/ (mi +mj). In addition, Fi denotes an external conservative
force acting on species i (such as a gravity field) and Fi is an operator repre-
senting a possible effect of an external nonconservative forcing which injects
energy into the system to compensate for the energy dissipated by collisional
cooling. This type of force acts as a thermostat that tries to mimics a thermal
bath. Some explicit forms for the operator Fi will be chosen later.
6 Vicente Garzó
The relevant hydrodynamic fields for the mixture are the number densities
ni, the flow velocity u, and the temperature T . They are defined in terms of
moments of the velocity distribution functions fi as
dvfi(v) , (4)
dvvfi(v) , (5)
nT = p =
niTi =
dvV 2fi(v) , (6)
where V = v − u is the peculiar velocity, n = n1 + n2 is the total number
density, ρ = m1n1 + m2n2 is the total mass density, and p is the pressure.
Furthermore, the third equality of Eq. (6) defines the kinetic temperatures Ti
for each species, which measure their mean kinetic energies.
The collision operators conserve the particle number of each species and
the total momentum, but the total energy is not conserved:
dvJij [v|fi, fj ] = 0 , (7)
dvvJij [v|fi, fj] +mj
dvvJji[v|fj , fi] = 0 , (8)
dvV 2Jij [v|fi, fj] = −dnTζ . (9)
In Eq. (9), ζ is identified as the total “cooling rate” due to collisions among
all species. It measures the rate of energy loss due to dissipation. At a kinetic
level, it is also convenient to introduce the “cooling rates” ζi for the partial
temperatures Ti. They are defined as
ζij = −
dniTi
dvV 2Jij [v|fi, fj], (10)
where the second equality defines the quantities ζij . The total cooling rate ζ
can be written in terms of the partial cooling rates ζi as
ζ = T−1
xiTiζi, (11)
where xi = ni/n is the mole fraction of species i.
From Eqs. (7)–(9), the macroscopic balance equations for the number den-
sities ni, the total momentum density ρu and the energy density (d/2)nT can
be obtained. They are given, respectively, by [19]
Kinetic Theory for Binary Granular Mixtures at Low-Density 7
Dtni + ni∇ · u+
∇ · ji
= 0 , (12)
ρDtu+∇ · P =
niFi , (13)
DtT −
∇ · ji
∇ · q+ P : ∇u−
Fi · ji
= −(ζ − ξ)T . (14)
In the above equations, Dt = ∂t + u · ∇ is the material derivative,
ji = mi
dvV fi(v) (15)
is the mass flux for species i relative to the local flow,
dvVV fi(v) (16)
is the total pressure tensor, and
dv V 2V fi(v) (17)
is the total heat flux. On the right-hand side of the temperature equation
(14), the source term ξ (measuring the rate of heating due to the external
thermostat) is given by
ξ = −
dvV 2Fifi(v). (18)
In the balance equations (12)–(14) it is assumed that the external driving
thermostat does not change the number of particles of each species or the
total momentum, i.e., ∫
dvFifi(v) = 0, (19)
dv v Fifi(v) = 0. (20)
The utility of the balance equations (12)–(14) is limited without further
specification of the fluxes and the cooling rate, which in general have a com-
plex dependence on space and time. However, for sufficiently large space and
time scales, one expects that the system reaches a hydrodynamic regime in
which all the space and time dependence is given entirely through a functional
dependence on the six hydrodynamic fields ni, u, and T . The corresponding
8 Vicente Garzó
functional dependence of ji, P, and q on these fields are called constitutive
equations and define the transport coefficients of the mixture. The primary
feature of a hydrodynamic description is the reduction of the description from
many microscopic degrees of freedom to a set of equations involving only six
local fields. At a kinetic level, the constitutive equations are obtained when
one admits the existence of a normal solution to the Boltzmann equation
where the velocity distribution functions depend on r and t only through
their functional dependence on the fields, namely,
fi(r,v1, t) = fi[v1|ni(r, t), T (r, t),u(r, t)]. (21)
This normal solution is generated by the Chapman–Enskog method [4] con-
veniently adapted to dissipative dynamics. Since the method is based on an
expansion around the local version of the homogeneous state, let us charac-
terize it before considering inhomogeneous solutions.
3 Homogeneous States
In this Section we are interested in spatially homogeneous isotropic states. In
this case, we assume that the magnitude of the conservative external forces
is at least of first order in the spatial gradients (i.e., Fi = 0), so that Eq. (1)
becomes
(∂t + Fi) fi =
Jij [fi, fj ]. (22)
For elastic collisions (αij = 1) and in the absence of external forcing (Fi = 0),
it is well known that the long-time solution to (22) is a Maxwellian distribution
for each species at the same temperature T . However, if the particles collide
inelastically (αij < 1) and Fi = 0, a steady state is not possible in uniform
situations since the temperature decreases monotonically in time. This state is
usually referred to as the homogeneous cooling state (HCS). In this situation,
since ni is uniform and u = 0, the normal (hydrodynamic) solution to fi
requires that all its time dependence occurs only through the temperature
T (t). Consequently, fi(v, t) must be of the form [15]
fi(v, t) = niv
0 (t)Φi (v/v0(t)) , (23)
where v0(t) =
2T (t)(m1 +m2)/ (m1m2) is a thermal speed defined in terms
of the temperature T (t) of the mixture. The balance equations (12)–(14) to
this order become ∂txi = ∂tu = 0, and T
−1∂tT = −ζ, where the cooling rate
ζ is determined by Eq. (9). In addition, from Eqs. (10) and (22) (when Fi = 0)
one can derive the time evolution for the temperature ratio γ = T1(t)/T2(t):
∂t ln γ = ζ2 − ζ1. (24)
Kinetic Theory for Binary Granular Mixtures at Low-Density 9
The fact that the distributions fi depend on time only through T (t) necessarily
implies that the temperature ratio γ must be independent of time and so, Eq.
(24) gives the condition
ζ1 = ζ2 = ζ. (25)
In the elastic case, where fi is a Maxwellian distribution, the above condition
yields T1(t) = T2(t) = T (t) and the energy equipartition applies. However, in
the inelastic case, the equality of the cooling rates leads to different values for
the partial temperatures, even if one considers the Maxwellian approximation
to fi. Nevertheless, the constancy of γ assures that the time dependence of
the distributions is entirely through T (t), and in fact the partial temperatures
can be expressed in terms of the global temperature as
T1(t) =
1 + x1(γ − 1)
T (t), T2(t) =
1 + x1(γ − 1)
T (t). (26)
Just as for the single gas case [30, 32], the exact form of Φi has not yet been
found, although a good approximation for thermal velocities can be obtained
from an expansion in Sonine polynomials [15]. In the leading order, Φi is given
e−θiv
∗4 − (d+ 2)θiv
d(d+ 2)
Here, v∗ ≡ v/v0,
m−1j , (28)
and γi = Ti/T . The coefficients ci (which measure the deviation of Φi from
the reference Maxwellian) are determined consistently from the Boltzmann
equation. The explicit form of the approximation (27) provides detailed pre-
dictions for the temperature ratio T1/T2 (through calculation of the cooling
rates) and for the cumulants ci as functions of the mass ratio, size ratio, com-
position and coefficients of restitution [15]. The numerical accuracy of this
truncated Sonine expansion has been confirmed by comparison with Monte
Carlo [33] and MD [34] simulations.
As said in the Introduction, the existence of different temperatures for
each species has been observed in real experiments of driven steady states.
These states are achieved from external forces that do work at the same rate
as collisional cooling. In experiments this is accomplished by vibrating the
system so that it is locally driven at walls. Far from these walls a steady state
is reached whose properties are expected to be insensitive to the details of
the driving forces. Due to the technical difficulties involved in incorporating
oscillating boundary conditions, it is usual to introduce external forces (or
thermostats) acting locally on each particle. These forces are represented by
the operator Fi in Eq. (22) and depend on the state of the system. Two
10 Vicente Garzó
types of thermostats have been usually considered in the literature. One is a
deterministic thermostat widely used in nonequilibrium MD simulations for
ordinary fluids [35, 36]. The force is similar to a Stokes law drag force, linear
in the velocity, but with the opposite sign so that it heats rather than cools
the system. In this case, Fi is given by [21, 37]
Fifi(v) =
· [vfi(v)], (29)
where the friction constant in the force has been adjusted to get a constant
temperature in the long-time limit. It must be remarked that the correspond-
ing Boltzmann equation (22) for this Gaussian thermostat force is formally
identical with the Boltzmann equation in the HCS (i.e., with Fi = 0) when
both equations are written in terms of the reduced distributions Φi(v
∗) [37].
In particular, the dependence of γ on the parameters of the system is the same
with and without the Gaussian thermostat.
A second method of driving the system is by means of a stochastic
Langevin force representing Gaussian white noise [38]. The corresponding op-
erator Fi has a Fokker–Planck form [32]
Fifi(v) = −
fi(v), (30)
where for simplicity the covariance of the stochastic acceleration has been
taken to the same for each species [39, 40]. This requirement gives the steady
state condition
ζ2. (31)
The cooling rates ζi are no longer equal, in contrast to the HCS, and the de-
pendence of γ on the control parameters is different as well [34]. The procedure
for determining the temperature ratio and the cumulants ci is the same as in
the HCS state since the steady state distribution Φi can also be represented as
an expansion of the form (27) and the coefficients are now determined from
the solution to the Boltzmann equation (22). The condition (31) gives the
corresponding equation for the temperature ratio.
Figure 1 illustrates the differences between the HCS and the stochastic
steady state at the level of the temperature ratio T1/T2. We have considered
mixtures constituted by spheres (d = 3) of the same material [and so, αij = α,
and m1/m2 = (σ1/σ2)
3] and equal volumes of large and small particles [i.e.,
x2 = (σ1/σ2)
3x1]. Here, for the sake of simplicity, the cooling rates have been
analytically estimated by using Maxwellians (namely, by taking ci = 0) for
the distributions fi(v):
ζij →
4π(d−1)/2
njµjiσ
θi + θj
×(1 + αij)
(1 + αij)
θi + θj
. (32)
Kinetic Theory for Binary Granular Mixtures at Low-Density 11
1.0 1.5 2.0 2.5 3.0
Fig. 1. Temperature ratio T1/T2 versus the size ratio σ1/σ2 for αij ≡ α = 0.78 in
the case of mixtures constituted by particles of the same material and equal total
volumes of large and small particles. The lines are the kinetic theory results in (a)
the stochastic driving case and (b) the free cooling case, while the points refer to
MD simulations [41].
Simulation data recently obtained from MD simulations in agitated mixtures
have also been included [41]. The experimental value of the (common) coef-
ficient of restitution is α = 0.78. While a good agreement between kinetic
theory and MD simulations is found when the gas is assumed to be driven
by the stochastic thermostat, significant discrepancies appear in the undriven
(HCS) case, especially as the size ratio σ1/σ2 increases. These results con-
trast with the comparison made by Brey et al. [42] for agitated mixtures in
the tracer limit (x1 → 0), where the predictions of the temperature ratio from
kinetic theory based on the condition ζ1 = ζ2 compare quite well with MD
simulations. However, for the cases studied in Ref. [42] the conditions (25) and
(31) yield quite similar results for the dependence of T1/T2 on the parameters
of the system. The good agreement found in Fig. 1 between MD simulations
for agitated mixtures and kinetic theory suggests that the stochastic driving
condition can be considered as a plausible first approximation for qualitative
comparisons with experimental results [17].
4 Navier–Stokes Hydrodynamic Equations
We consider now a spatially inhomogeneous state created either by initial
preparation or by boundary conditions. We are interested in a hydrodynamic
description where the state of the system is completely specified through their
hydrodynamic fields. This implies that the latter dominate over other excita-
12 Vicente Garzó
tions for times large compared to the mean free time and for wavelengths large
compared to the mean free time. The hydrodynamic regime is characterized
by the existence of a normal solution to the Boltzmann equation which can be
explicitly obtained by means of the Chapman–Enskog method [4]. For small
spatial variations, the functional dependence (21) of the normal solution can
be made local in space and time through an expansion in gradients of the
fields:
fi = f
i + ǫ f
i + · · · , (33)
where each factor of ǫ means an implicit gradient of a hydrodynamic field.
The reference distribution function f
i (r,v, t) is the local version of the ho-
mogeneous distribution (23), namely, it is obtained from the homogeneous
distribution by replacing the temperature, densities and flow velocity by their
nonequilibrium values:
i (r,v, t) = ni(r, t)v
0 (T (r, t))Φi (V/v0(T (r, t))) , (34)
where V = v−u(r, t). The time derivatives of the fields are also expanded as
∂t = ∂
t + ǫ∂
t + · · · . The coefficients of the time derivative expansion are
identified from the balance equations (12)–(14) with a representation of the
fluxes and the cooling rate in the macroscopic balance equations as a similar
series through their definitions as functionals of fi.
This is the usual Chapman–Enkog method for solving kinetic equations
[4, 28, 43]. Nevertheless, the complexity introduced by the energy dissipation
in collisions has led to the introduction by some authors [14] of some additional
approximations, restricting the validity of most of the results to the small
inelasticity limit. Only very recently, explicit expressions for the fluxes to first
order in the gradients as explicit functions of the coefficients of restitution have
been obtained [19, 21, 44]. To Navier–Stokes order, the constitutive equations
for mass, momentum, and heat fluxes are given, respectively, by
1 = −
m1m2n
D∇x1 −
Dp∇p−
D′∇T +
χ1iFi, j
2 = −j
kℓ = p δkℓ − η
∇ℓuk +∇kuℓ −
δkℓ∇ · u
, (36)
(1) = −T 2D′′∇x1 − L∇p− λ∇T +
κiFi. (37)
The transport coefficients in these equations are the diffusion coefficient D,
the pressure diffusion coefficient Dp, the thermal diffusion coefficient D
′, the
mobility coefficient χij , the shear viscosity η, the Dufour coefficient D
′′, the
pressure energy coefficient L, the thermal conductivity λ, and the coefficient
κi. These coefficients are defined as
Kinetic Theory for Binary Granular Mixtures at Low-Density 13
D = −
dvV ·A1, (38)
Dp = −
dvV · B1, (39)
D′ = −
dvV · C1, (40)
χij =
dvV · E ij , (41)
η = −
(d− 1)(d+ 2)
dvVV : Di, (42)
D′′ = −
dv V 2V ·Ai, (43)
L = −
dv V 2V · Bi, (44)
λ = −
dv V 2V · Ci, (45)
dv V 2V · Eij . (46)
Here, Ai(V), Bi(V), Ci(V), Di(V), and Eij(V) are functions of the peculiar
velocity and the hydrodynamic fields. They obey the following set of coupled
linear integral equations:
ζ (T∂T + p∂p) + F
1 + L1
A1 +M1A2 = A1 +
× (pB1 + TC1) , (47a)
ζ (T∂T + p∂p) + F
2 + L2
A2 +M2A1 = A2 +
× (pB2 + TC2) , (47b)
ζ (T∂T + p∂p) + F
1 + L1 − 2c
B1+M1B2 = B1+
C1, (48a)
ζ (T∂T + p∂p) + F
2 + L2 − 2c
B2+M2B1 = B2+
C2, (48b)
14 Vicente Garzó
ζ (T∂T + p∂p) + F
1 + L1 −
C1+M1C2 = C1−
B1, (49a)
ζ (T∂T + p∂p) + F
2 + L2 −
C2+M2C1 = C2−
B2, (49b)
ζ (T∂T + p∂p) + F
1 + L1
D1 +M1D2 = D1, (50a)
ζ (T∂T + p∂p) + F
2 + L2
D2 +M2D1 = D2, (50b)
ζ (T∂T + p∂p) + F
1 + L1
E11 +M1E21 = E11, (51a)
ζ (T∂T + p∂p) + F
1 + L1
E12 +M1E22 = E12. (51b)
Here, we have introduced the linearized Boltzmann collision operators
L1X = −
J11[f
1 , X ] + J11[X, f
1 ] + J12[X, f
, (52)
M1X = −J12[f
1 , X ], (53)
with a similar definition for L2 andM2. In addition, c
ζ = ζ
(0) in the undriven
case while c
ζ = 0 in the driven case, where ζ
(0) is given by Eq. (9) to zeroth
order, i.e.,
ζ(0) = −
dvV 2Jij [v|f
i , f
j ] . (54)
In Eqs. (47a)–(51b) we have also introduced the operators
i X =
· (VX) , (Gaussian thermostat), (55a)
i X = −
X, (stochastic thermostat), (55b)
and the quantities
Ai(V) = −
V, (56)
Bi(V) = −
i V +
, (57)
Ci(V) =
V, (58)
Kinetic Theory for Binary Granular Mixtures at Low-Density 15
Di(V) = V
I, (59)
Eij(V) = −
δij −
. (60)
Here, I is the unit tensor in d dimensions and ρi = mini is the mass density
of species i. Upon writing Eqs. (47a)–(51b) use has been made of the fact
that there is no contribution to the cooling rate at this order, i.e., ζ(1) = 0.
As a consequence, F
i = 0. The property ζ
(1) = 0 is special of the low
density Boltzmann kinetic theory (since f
i does not contain any contribution
proportional to ∇ ·u), but such terms occur at higher densities [45, 46]. Note
that in the particular case of the gravitational force Fi = mig (where g is
the gravity acceleration), the combination m1E11+m2E12 = 0. This leads to
Eij = 0, and so there are no contributions to the mass and heat fluxes coming
from the external conservative forces.
In summary, the solutions to the Boltzmann equations to first order in the
spatial gradients are given by [19]
fi = f
i +Ai · ∇x1 +Bi · ∇p+ Ci · ∇T +Di : ∇u+
Eij · Fj . (61)
The solution to zeroth-order is obtained from Eq. (22) while the functions
{Ai,Bi,Ci,Di,Eij} are determined from the integral equations (47a)–(51b).
Once these equations are solved, the Navier-Stokes transport coefficients are
obtained from Eqs. (38)–(46) and the mass, momentum, and heat fluxes are
explicitly known. These fluxes, together with the macroscopic balance equa-
tions (12)–(14), provide the closed set of Navier-Stokes order hydrodynamic
equations for a granular binary mixture. All these results are still formally
exact and valid for arbitrary values of the coefficients of restitution.
However, explicit expressions for the Navier–Stokes transport coefficients
require to solve the integral equations (47a)–(51b). Accurate approximations
for {Ai,Bi,Ci,Di,Eij} may be obtained using low order truncation of ex-
pansions in a series of Sonine polynomials. The polynomials are defined with
respect to a Gaussian weight factor whose parameters are chosen such that
the leading term in the expansion yields the exact moments of the entire dis-
tribution with respect to 1, v, and v2. The procedure is similar to the one
followed for elastic collisions [4] and yields explicit expressions for the trans-
port coefficients in terms of the parameters of the mixture [19, 20, 44].
5 Comparison with Monte Carlo Simulations
As said before, the expressions derived for the Navier Stokes transport coef-
ficients are obtained by considering two different approximations. First, since
the deviation of f
i from its Maxwellian form is quite small in the region of
16 Vicente Garzó
thermal velocities, one uses the distribution (27) as a trial function for f
Second, one only considers the leading terms of an expansion of the distri-
bution f
i in Sonine polynomials. Both approximations allow one to offer
a simplified kinetic theory for a granular binary mixture. To assess the de-
gree of accuracy of these predictions, one resorts to numerical solutions of the
Boltzmann equation, such as those obtained from the Direct simulation Monte
Carlo (DSMC) method [47]. Although the method was originally devised for
normal fluids, its extension to granular gases is straightforward [37]. In this
Section we provide some comparisons between theory and numerical solutions
of the Boltzmann equation by means of the DSMC method in the cases of the
diffusion coefficient D (in the tracer limit) and the shear viscosity coefficient
η of a heated gas. Let us study each coefficient separately.
5.1 Tracer Diffusion Coefficient
We consider a free granular mixture (Fi = 0) in which one of the components
of the mixture (say, for instance, species 1) is present in tracer concentration
(x1 → 0). In this situation the diffusion coefficient of impurities in a granular
gas undergoing homogeneous cooling state can be measured in simulations
from the mean square displacement of the tracer particle after a time interval
t [48]:
〈|r(t) − r(0)|2〉 =
, (62)
where |r(t) − r(0)| is the distance travelled by the impurity from t = 0 until
time t. The relation (62) written in appropriate dimensionless variables to
eliminate the time dependence of D(t) can be used to measure by computer
simulations the diffusion coefficient [8, 22].
If the hydrodynamic description (or normal solution in the context of the
Chapman–Enskogmethod) applies, then the diffusion coefficientD(t) depends
on time only through its dependence on the temperature T (t). Dimensional
analysis shows that D(t) ∝
T (t). In this case, after a transient regime,
the reduced diffusion coefficient D∗ = (m1m2/ρ)D(t)ν0(t)/T (t) achieves a
time-independent value. Here, ν0(t) = nσ
12 v0(t) ∝
T (t) is an effective
collision frequency for hard spheres. The fact that D∗ reaches a constant
value for times large compared with the mean free path is closely related
with the validity of a hydrodynamic description for the system. In addition,
as has been recently shown [49], the dependence of D∗ on the mass ratio
m1/m2 and the coefficient of restitution α12 is only through the effective
mass m∗1 = m1 + (m1 +m2)(1 − α12)/(1 + α12).
The dependence of D∗ on the common coefficient of restitution αij ≡ α
is shown in Fig. 2 in the case of hard disks (d = 2) for three different sys-
tems. The symbols refer to DSMC simulations while the lines correspond to
the kinetic theory results obtained in the first Sonine approximation [19, 44].
MD results reported in Ref. [8] when impurities and particles of the gas are
Kinetic Theory for Binary Granular Mixtures at Low-Density 17
0.6 0.7 0.8 0.9 1.0
2.5 ω=µ
µ=1/4
Fig. 2. Plot of the reduced diffusion coefficient D∗ as a function of the (common)
coefficient of restitution α for binary mixtures with ω = µ in the case of a two-
dimensional system (d = 2). Here, ω ≡ σ1/σ2 and µ ≡ m1/m2. The symbols are
computer simulation results obtained from the mean square displacement and the
lines are the theoretical results obtained in the first Sonine approximation. The
DSMC results correspond to µ = 1/4 (•), µ = 4 (◦) and µ = 1 (⋄). Molecular
dynamics results reported in Ref. [8] for µ = 1 (△) have also been included.
mechanically equivalent have also been included. We observe that in the latter
case MD and DSMC results are consistent among themselves in the range of
values of α explored. This good agreement gives support to the applicability of
the inelastic Boltzmann equation beyond the quasielastic limit. It is apparent
that the agreement between the first Sonine approximation and simulation
results is excellent when impurities and particles of the gas are mechanically
equivalent and when impurities are much heavier and/or much larger than
the particles of the gas (Brownian limit). However, some discrepancies be-
tween simulation an theory are found with decreasing values of the mass ratio
µ ≡ m1/m2 and the size ratio ω ≡ σ1/σ2. These discrepancies are not easily
observed in Fig. 2 because of the small magnitude of D∗ for µ = 1/4. For
these systems, the second Sonine approximation [22] improves the qualita-
tive predictions over the first Sonine approximation for the cases in which
the gas particles are heavier and/or larger than impurities. This means that
the Sonine polynomial expansion exhibits a slow convergence for sufficiently
small values of the mass ratio µ and/or the size ratio ω. This tendency is also
present in the case of elastic systems [50].
18 Vicente Garzó
5.2 Shear Viscosity Coefficient of a Heated Gas
The shear viscosity η is perhaps the most widely studied transport coefficient
in granular fluids. This coefficient can be measured in computer simulations
in the special hydrodynamic state of uniform shear flow (USF). At a macro-
scopic level, this state is characterized by constant partial densities ni, uni-
form temperature T , and a linear flow velocity profile u1,k = u2,k = akℓrℓ,
akℓ = aδkxδℓy, a being the constant shear rate. In this state, the temperature
changes in time due to the competition between two mechanisms: on the one
hand, viscous heating and, on the other hand, energy dissipation in collisions.
In addition, the mass and heat fluxes vanish by symmetry reasons and the
(uniform) pressure tensor is the only nonzero flux of the problem. The relevant
balance equation is that for temperature, Eq. (14), which reduces to
∂tT +
aPxy = −(ζ − ξ)T, (63)
where
Pxy =
dVVxVyfi(V) (64)
is the xy-element of the pressure tensor.
For a granular fluid under USF and in the absence of a thermostatting
force (ξ = 0), the energy balance equation (63) leads to a steady state when
the viscous heating effect is exactly balanced by the collisional cooling. This
situation will be analyzed in Sec. 10. However, if for instance the mixture
is heated by the Gaussian thermostat (29) (with v → V), then the viscous
heating still prevails so that the temperature increases in time. In this case,
the collision frequency ν0(t) ∝
T (t) also grows with t and hence the reduced
shear rate a∗(t) = a/ν0(t) (which is the relevant nonequilibrium parameter
of the problem) monotonically decreases in time. Under these conditions, the
system asymptotically achieves a regime described by linear hydrodynamics
and the (reduced) shear viscosity η∗ = [ν0(t)/nT (t)]η(t) can be measured as
η∗ = − lim
P ∗xy
, (65)
where P ∗xy = Pxy/nT . This procedure allows one to identify the shear viscosity
of a granular mixture excited by the Gaussian external force (29) and compare
it with the predictions given by the Chapman–Enskog method.
In Fig. 3, we plot the ratio η∗(α)/η∗(1) versus the mass ratiom1/m2 in the
case of hard spheres (d = 3) for σ1/σ2 = 1, x1 =
, and three different values
of the (common) coefficient of restitution αij ≡ α. Here, η
∗(1) refers to the
elastic value of the shear viscosity coefficient. Again, the symbols represent
the simulation data obtained by numerically solving the Boltzmann equation
[21], while the lines refer to the theoretical results obtained from the Boltz-
mann equation in the first Sonine approximation. We see that in general the
Kinetic Theory for Binary Granular Mixtures at Low-Density 19
0 2 4 6 8 10
η*(α)/η
Fig. 3. Plot of the ratio η∗(α)/η∗(1) as a function of the mass ratio m1/m2 for
σ1/σ2 = n1/n2 = 1 and three different values of the (common) coefficient of resti-
tution α: α = 0.9 (circles), α = 0.8 (squares), and α = 0.7 (triangles). The lines are
the theoretical predictions and the symbols refer to the results obtained from the
DSMC method.
1 2 3 4 5
η*(α)/η
Fig. 4. Plot of the ratio η∗(α)/η∗(1) as a function of the size ratio σ1/σ2 for
m1/m2 = 4, n1/n2 = 1 and three different values of the (common) coefficient of
restitution α: α = 0.9 (circles), α = 0.8 (squares), and α = 0.7 (triangles). The lines
are the theoretical predictions and the symbols refer to the results obtained from
the DSMC method.
20 Vicente Garzó
deviation of η∗(α) from its functional form for elastic collisions is quite impor-
tant. This tendency becomes more significant as the mass disparity increases.
The agreement between the first Sonine approximation and simulation is seen
to be in general excellent. This agreement is similar to the one previously
found in the monocomponent case [7, 10, 51]. At a quantitative level, the dis-
crepancies between theory and simulation tend to increase as the coefficient
of restitution decreases, although these differences are quite small (say, for
instance, around 2% at α = 0.7 in the disparate mass case m1/m2 = 10).
The influence of the size ratio on the shear viscosity is shown in Fig. 4 for
m1/m2 = 4 and x1 =
[21]. We observe again a strong dependence of the
shear viscosity on dissipation. However, for a given value of α, the influence of
σ1/σ2 on η
∗ is weaker than the one found before in Fig. 3 for the mass ratio.
The agreement for both α = 0.9 and α = 0.8 is quite good, except for the
largest size ratio at α = 0.8. These discrepancies become more significant as
the dissipation increases, especially for mixtures of particles of very different
sizes. In summary, according to the comparison carried out in Figs. 3 and 4,
one can conclude that the agreement between theory and simulation extends
over a wide range values of the coefficient of restitution, indicating the relia-
bility of the first Sonine approximation for describing granular flows beyond
the quasielastic limit.
6 Einstein Relation in Granular Gases
The results presented in Section 5 give some support to the validity of the
hydrodynamic description to granular fluids. However, in spite of this sup-
port some care is warranted in extending properties of normal fluids to those
with inelastic collisions. Thus, for elastic collisions, in the case of an impurity
(tracer) particle immersed in a gas the response to an external force on the
impurity particle leads to a mobility coefficient proportional to the diffusion
coefficient. This is the usual Einstein relation [48], which is a consequence of
the fluctuation-dissipation theorem. A natural question is whether the Ein-
stein relation also applies for granular fluids.
To analyze it, let us consider the tracer limit (x1 → 0) and assume that
the current of impurities j
1 is only generated by the presence of a weak
concentration gradient ∇x1 and/or a weak external field F1 acting only on
the impurity particles. Under these conditions, Eq. (35) becomes
1 = −m1D∇x1 + χ11F1. (66)
The Einstein ratio ǫ′ between the diffusion coefficient D and the mobility
coefficient χ11 is defined as
ǫ′ = m1x1
, (67)
Kinetic Theory for Binary Granular Mixtures at Low-Density 21
where T ≃ T2 in the tracer limit. For elastic collisions, the Chapman–Enskog
results yield ǫ′ = 1. However, at finite inelasticity the relationship between D
and χ11 is no longer simple and, as expected, the Chapman–Enskog expres-
sions for D and χ11 in the case of an unforced granular gas [23] clearly show
that ǫ′ 6= 1. This means that the Einstein relation does not apply in granular
gases. The deviations of the (standard) Einstein ratio ǫ′ from unity has three
distinct origins: the absence of the Gibbs state (non-Gaussianity of the distri-
bution function of the HCS), time evolution of the granular temperature, and
the occurrence of different kinetic temperatures between the impurity and gas
particles. The second source of discrepancy can be avoided if the system is
driven by an external energy input to achieve a stationary state. With respect
to the third reason of violation, this could also be partially eliminated if the
temperature of the gas T is replaced by the temperature of the impurity T1
in the usual Einstein relation (67). This change yields the modified Einstein
ratio
ǫ = m1x1
T1χ11
. (68)
As a consequence, the only reason for which ǫ 6= 1 is due to the non-Maxwellian
behavior of the HCS distribution. Given that the deviations of the gas dis-
tribution f
2 from its Maxwellian form are small [15], the discrepancies of ǫ
from unity could be difficult to detect in computer simulations. This conclu-
sion agrees with recent MD simulations [52] of granular mixtures subjected to
the stochastic driving of the form (30), where no deviations from the (modi-
fied) Einstein relation ǫ = 1 have been observed for a wide range of values of
the coefficients of restitution and parameters of the system.
To illustrate the influence of dissipation on the Einstein ratio more gener-
ally, in Figs. 5 and 6 the Einstein ratio as given by (68) is plotted versus the
coefficient of restitution α12 for σ1/σ2 = 1 and different values of the mass ra-
tio m1/m2 and the coefficient of restitution α22. The results obtained by using
the Gaussian thermostat (29) are shown in Fig. 5, while Fig. 6 corresponds
to the results derived when the system is heated by the stochastic thermostat
(30) [24]. We observe that in general ǫ 6= 1, although its value is very close to
unity, especially in the case of the stochastic thermostat, where the deviations
from the Einstein relation are smaller than 1%. However, in the case of the
Gaussian thermostat the deviations from unity are about 8%, which could be
detected in computer simulations. Figures 5 and 6 also show the fact that the
transport properties are affected by the thermostat introduced so that the
latter does not play a neutral role in the problem [51].
7 Onsager’s Reciprocal Relations in Granular Gases
In the usual language of the linear irreversible thermodynamics for ordinary
fluids [53], the constitutive equations (35) and (37) for the mass flux and heat
flux in the absence of external forces can be written as
22 Vicente Garzó
0.0 0.2 0.4 0.6 0.8 1.0
Fig. 5. Plot of the modified Einstein ratio ǫ versus the coefficient of restitution α12
for the Gaussian thermostat in the cases: (a) α22 = α12, m1/m2 = 5 and σ1/σ2 = 1;
(b) α22 = α12, m1/m2 = 0.5 and σ1/σ2 = 1; and (c) α22 = 0.5, m1/m2 = 10 and
σ1/σ2 = 1.
0.0 0.2 0.4 0.6 0.8 1.0
(c) (b)
Fig. 6. Plot of the Einstein ratio ǫ versus the coefficient of restitution α12 for the
stochastic thermostat in the cases: (a) α22 = α12, m1/m2 = 5 and σ1/σ2 = 1; (b)
α22 = α12, m1/m2 = 0.5 and σ1/σ2 = 1; and (c) α22 = 0.5, m1/m2 = 10 and
σ1/σ2 = 1.
Kinetic Theory for Binary Granular Mixtures at Low-Density 23
i = −
− Liq
− Cp∇p, (69)
q = −Lqq∇T −
− C′p∇p, (70)
where
q ≡ q
(1) −
and (
∇ ln(xip), (72)
µi being the chemical potential per unit mass. In Eqs. (69) and (70), the
coefficients Lij are the so-called Onsager phenomenological coefficients and the
coefficients Cp and C
p can be expressed in terms of the transport coefficients
associated with the heat and mass fluxes. For elastic fluids, Onsager showed
[53] that time reversal invariance of the underlying microscopic equations of
motion implies important constraints on the above set of transport coefficients,
namely
Lij = Lji, Liq = Lqi, Cp = C
p = 0. (73)
The first two symmetries are called reciprocal relations as they relate trans-
port coefficients for different processes. The last two are statements that the
pressure gradient does not appear in any of the fluxes, even though it is ad-
mitted by symmetry. Even for a one component fluid, Onsager’s theorem is
significant as it leads to a new contribution to the heat flux proportional to
the density gradient [5]. Since there is no time reversal symmetry for granular
fluids, Eqs. (73) cannot be expected to apply. However, since explicit expres-
sions for all transport coefficients are at hand, the quantitative extent of the
violation can be explored.
To make connection with the expressions (35) and (37) for the mass and
heat fluxes, respectively, it is first necessary to transform Eqs. (69)–(71) to
the variables x1, p, T. Since ∇x1 = −∇x2, Eq. (72) implies
(∇µ1)T − (∇µ2)T
∇x1 +
(m2 −m1)∇ ln p
. (74)
The coefficients {Lij, Liq, Lqi, Lqq, Cp, C
p} then can be easily obtained in
terms of the Navier-Stokes transport coefficients introduced in Sec. 4. The
result is [20]
L11 = −L12 = −L21 =
m1m2ρ1ρ2
D, L1q = ρTD
′, (75)
Lq1 = −Lq2 =
T 2ρ1ρ2
D′′ −
Tρ1ρ2
(m2 −m1)D, (76)
24 Vicente Garzó
Lqq = λ−
m2 −m1
D′, (77)
(m2 −m1)D, (78)
C′p ≡ L−
m2 −m1
T 2(m2 −m1)D
′′. (79)
Onsager’s relation L12 = L21 holds since the diffusion coefficient D is
symmetric under the change 1 ↔ 2 [19]. However, in general L1q 6= Lq1, Cp 6=
0, and C′p 6= 0 [20]. The Chapman-Enskog results [19] show that there are only
two limit cases for which L1q−Lq1 = Cp = C
p = 0: (i) the elastic limit (αij =
1) with arbitrary values of masses, sizes and composition and (ii) the case
of mechanically equivalent particles with arbitrary values of the (common)
coefficient of restitution α ≡ αij . Beyond these limit cases, Onsager’s relations
do not apply. At macroscopic level the origin of this failure is due to the cooling
of the reference state as well as the occurrence of different kinetic temperatures
for both species.
8 Linearized Hydrodynamic Equations and Stability of
the Homogeneous Cooling State
As shown in Sec. 4, the Navier-Stokes constitutive equations (35)–(37) have
been expressed in terms of a set of experimentally accessible fields such as
the composition of species 1, x1, the pressure p, the mean flow field u, and
the granular temperature T . In terms of these variables and in the absence of
external forces, the macroscopic balance equations (12)–(14) become
Dtx1 +
n2m1m2
∇ · j1 = 0 , (80)
Dtp+ p∇ · u+
(∇ · q+ P : ∇u) = −ζp, (81)
Dtu+ ρ
−1∇ · P = 0 , (82)
DtT −
∇ · ji
(∇ · q+ P : ∇u) = −ζT . (83)
When the expressions (35)–(37) for the fluxes and the cooling rate ζ → ζ(0) are
substituted into the above exact balance equations (80)–(83) one gets a closed
set of hydrodynamic equations for x1,u, T, and p. These are the Navier–Stokes
hydrodynamic equations for a binary granular mixture:
Dtx1 =
n2m1m2
m1m2n
D∇x1 +
Dp∇p+
, (84)
Kinetic Theory for Binary Granular Mixtures at Low-Density 25
(Dt + ζ) p+
p∇ · u =
T 2D′′∇x1 + L∇p+ λ∇T
∇ℓuk +∇kuℓ −
δkℓ∇ · u
∇ℓuk,
(Dt + ζ) T +
p∇ · u = −
m2 −m1
m1m2n
D∇x1 +
T 2D′′∇x1 + L∇p+ λ∇T
∇ℓuk +∇kuℓ −
δkℓ∇ · u
∇ℓuk, (86)
Dtuℓ + ρ
−1∇ℓp = ρ
−1∇kη
∇ℓuk +∇kuℓ −
δkℓ∇ · u
. (87)
For the chosen set of fields, n = p/T and ρ = p [(m1 −m2)x1 +m2] /T . These
equations are exact to second order in the spatial gradients for a low density
Boltzmann gas. Note that in Eqs. (84)–(87) the second order contributions
to the cooling rate have been neglected. These second order terms have been
calculated for a monocomponent fluid [5] and found to be very small relative
to corresponding terms from the fluxes. Consequently, they have not been
considered in the hydrodynamic equations (84)–(87).
One of the main peculiarities of the granular gases (in contrast to ordinary
fluids) is the existence of non-trivial solutions to the Navier–Stokes equations
(84)–(87), even for spatially homogeneous states,
∂tx1H = 0, ∂tuH = 0, (88)
[∂t + ζ (x1H , TH , pH)]TH = 0, [∂t + ζ (x1H , TH , pH)] pH = 0, (89)
where the subscript H denotes the homogeneous state. Since the dependence
of the cooling rate ζ (x1H , TH , pH) on x1H , TH , pH is known [15, 44], these
first order nonlinear equations can be solved for the time dependence of the
homogeneous state. The result is the familiar Haff cooling law for T (t) at
constant density [1, 54]:
TH(t) =
TH(0)
[1 + ζ(0)t/2]
. (90)
As said before, each partial temperature Ti(t) has the same time dependence
but with a different value [15],
T1H(t) =
1 + x1(γ − 1)
TH(t), T2H(t) =
1 + x1(γ − 1)
TH(t), (91)
where γ = T1H(t)/T2H(t) is the time-independent temperature ratio.
26 Vicente Garzó
Nevertheless, the homogeneous cooling state (HCS) is unstable to suffi-
ciently long wavelength perturbations. For systems large enough to support
such spontaneous fluctuations, the HCS becomes inhomogeneous at long times.
This feature was first observed in MD simulations of free monocomponent
gases [55]. In MD simulations the inhomogeneities may grow by the forma-
tion of clusters, ultimately aggregating to a single large cluster [56]; if cluster
growth is suppressed, a vortex field may grow to the system size where peri-
odic boundary conditions can induce a transition to a state with macroscopic
shear. The mechanism responsible for the growth of inhomogeneities can be
understood at the level of the Navier–Stokes hydrodynamics, where linear
stability analysis shows two shear modes and a heat mode to be unstable
[5, 7, 54, 57].
The objective here is to extend this analysis to the case of a binary mixture.
To do that, we perform a linear stability analysis of the nonlinear hydrody-
namic equations (84)–(87) with respect to this HCS for small initial spatial
perturbations. For ordinary fluids such perturbations decay in time according
to the hydrodynamic modes of diffusion (shear, thermal, mass) and damped
sound propagation [3, 58, 59]. For inelastic collisions, the analysis is for fixed
coefficients of restitution in the long wavelength limit. As will be seen be-
low, the corresponding modes for a granular mixture are then quite different
from those for ordinary mixtures. In fact, an alternative study with fixed long
wavelength and coefficients of restitution approaching unity yields the usual
ordinary fluid modes. Consequently, the nature of the hydrodynamic modes
is non uniform with respect to the inelasticity and the wavelength of the
perturbation.
Let us assume that the deviations δyα(r, t) = yα(r, t) − yHα(t) are small.
Here, δyα(r, t) denotes the deviation of {x1,u, T, p} from their values in the
HCS. If the initial spatial perturbation is sufficiently small, then for some
initial time interval these deviations will remain small and the hydrodynamic
equations (84)–(87) can be linearized with respect to δyα(r, t). This leads to
a set of partial differential equations with coefficients that are independent
of space but which depend on time. As in the monocomponent case [5, 57],
this time dependence can be eliminated through a change in the time and
space variables, and a scaling of the hydrodynamic fields. We introduce the
following dimensionless space and time variables:
dt′ν0H(t
′), s =
ν0H(t)
v0H(t)
r, (92)
where ν0H(t) is an effective collision frequency for hard spheres and v0H =√
2TH(m1 +m2)/m1m2. Since {x1H ,uH , TH , pH} are evaluated in the HCS,
then Eqs. (88) and (89) hold. A set of Fourier transformed dimensionless
variables are then introduced as
ρk(τ) =
δx1k(τ)
, wk(τ) =
δuk(τ)
v0H(τ)
, θk(τ) =
δTk(τ)
TH(τ)
, Πk(τ) =
δpk(τ)
pH(τ)
Kinetic Theory for Binary Granular Mixtures at Low-Density 27
where δyαk ≡ {δx1k, δuk, δTk, δpk} is defined as
δyαk(τ) =
ds e−ik·sδyα(s, τ). (94)
Note that here the wave vector k is dimensionless.
In terms of the above variables, the transverse velocity components wk⊥ =
wk − (wk · k̂)k̂ (orthogonal to the wave vector k) decouple from the other
four modes and hence can be obtained more easily. They obey the equation
+ η∗k2
wk⊥ = 0, (95)
where ζ∗ = ζH/ν0H and
η, (96)
where ρH = m1n1H +m2n2H . The solution for wk⊥(τ) reads
wk⊥(τ) = wk⊥(0) exp[s⊥(k)τ ], (97)
where
s⊥(k) =
ζ∗ − η∗k2. (98)
This identifies d − 1 shear (transversal) modes. We see from Eq. (98) that
there exists a critical wave number kc⊥ given by
kc⊥ =
. (99)
This critical value separates two regimes: shear modes with k ≥ kc⊥ always
decay while those with k < kc⊥ grow exponentially.
The remaining modes are called longitudinal modes. They correspond to
the set {ρk, θk, Πk, wk||} where the longitudinal velocity component (parallel
to k) is wk|| = wk · k̂. These modes are the solutions of the linear equation
∂δzαk(τ)
αβ + ikM
αβ + k
δzβk(τ), (100)
where δzαk(τ) denotes now the four variables
ρk, θk, Πk, wk||
. The matrices
in Eq. (100) are given by
(0) =

0 0 0 0
ζ∗ −ζ∗ 0
ζ∗ −ζ∗ 0
0 0 0 1

, (101)
28 Vicente Garzó
(1) =
0 0 0 0
0 0 0 − 2
0 0 0 − d+2
0 0 − 1
x1µ+x2
 , (102)
(2) =
−D∗ −x−11 D
′∗ −x−11 D
′′∗ − 2
λ∗ − 2
0 0 0 − 2
(d− 1)η∗
 , (103)
where
21 = −x1
D′′∗ −
x1µ+ x2
, (104)
x1µ+ x2
D′∗ −
λ∗, (105)
23 = −
x1µ+ x2
D∗p. (106)
In these equations, µ = m1/m2, xi = niH/nH , and we have introduced the
reduced Navier–Stokes transport coefficients1
D, D∗p =
ρ2Hν0H
m1m2n
Dp, D
ρ2Hν0H
m1m2n
D′, (107)
D′′∗ =
ν0HTH
D′′, L∗ =
L, λ∗ =
λ. (108)
The longitudinal modes have the form exp[sn(k)τ ] with n = 1, 2, 3, 4,
where sn(k) are the eigenvalues of the matrix M(k) = M
(0)+ ikM(1)+k2M(2),
namely, they are the solutions of the quartic equation
det |M− sI| = 0. (109)
The solution to (109) for arbitrary values of k is quite intricate. It is instructive
to consider first the solutions to these equations in the extreme long wave-
length limit, k = 0. In this case, they are found to be the eigenvalues of the
matrix of M(0):
s(0)n =
0, 0,−
. (110)
Hence, at asymptotically long wavelengths (k = 0) the spectrum of the lin-
earized hydrodynamic equations (both transverse and longitudinal) is com-
prised of a decaying mode at −ζ∗/2, a two-fold degenerate mode at 0, and a
d-fold degenerate unstable mode at ζ∗/2. Consequently, some of the solutions
are unstable. The two zero eigenvalues represent marginal stability solutions,
1 Note that the definition for the reduced diffusion coefficient D∗ given here differs
from the one introduced in Sec. 5.1.
Kinetic Theory for Binary Granular Mixtures at Low-Density 29
0.0 0.1 0.2 0.3 0.4 0.5
-0.10
-0.05
Fig. 7. Dispersion relations for α = 0.9, x1 = 0.2, ω = 1 and µ = 4.
while the negative eigenvalue gives stable solutions. For general initial per-
turbations all modes are excited. These modes correspond to evolution of
the fluid due to uniform perturbations of the HCS, i.e., a global change in
the HCS parameters. The unstable modes are seen to arise from the initial
perturbations wk⊥(0) or wk||(0). The marginal modes correspond to changes
in the composition at fixed pressure, density, and velocity, and to changes in
Πk−θk at constant composition and velocity. The decaying mode corresponds
to changes in the temperature or pressure for Πk = θk. The unstable modes
may appear trivial since they are due entirely to the normalization of the
fluid velocity by the time dependent thermal velocity. However, this normal-
ization is required by the scaling of the entire set of equations to obtain time
independent coefficients.
The real parts of the modes s⊥(k) and sn(k) is illustrated in Fig. 7 in the
case of hard spheres (d = 3) for α ≡ αij = 0.9, σ1/σ2 = 1, x1 = 0.2, and
m1/m2 = 4. The k = 0 values correspond to five hydrodynamic modes with
two different degeneracies. The shear mode degeneracy remains at finite k but
the other is removed at any finite k. At sufficiently large k a pair of real modes
become equal and become a complex conjugate pair at all larger wave vectors,
like a sound mode. The smallest of the unstable modes is that associated with
the longitudinal velocity, which couples to the scalar hydrodynamic fields. It
becomes negative at a wave vector smaller than that of Eq. (99) and gives the
threshold for development of spatial instabilities.
The results obtained here for mixtures show no new surprises relative to
the case for a monocomponent gas [5, 54, 57], with only the addition of the
stable mass diffusion mode. Of course, the quantitative features can be quite
different since there are additional degrees of freedom with the parameter set
{x1H ,m1/m2, σ1/σ2, αij}. Also, the manner in which these linear instabilities
30 Vicente Garzó
are enhanced by the nonlinearities may be different from that for the one
component case [60].
9 Segregation in Granular Binary Mixtures: Thermal
Diffusion
The analysis of the linearized hydrodynamic equations for a granular binary
mixture has shown that the resulting equations exhibit a long wavelength
instability for d of the modes. These instabilities lead to the spontaneous for-
mation of velocity vortices and density clusters when the system evolves freely.
A phenomenon related with the density clustering is the separation or species
segregation. Segregation and mixing of dissimilar grains is perhaps one of the
most interesting problems in agitated granular mixtures. In some processes
it is a desired and useful effect to separate particles of different types, while
in other situations it is undesired and can be difficult to control. A variety
of mechanisms have been proposed to describe the separation of particles of
two sizes in a mixture of vertically shaken particles. Different mechanisms
include void filling, static compressive force, convection, condensation, ther-
mal diffusion, interstitial gas forcing, friction, and buoyancy [61]. However,
in spite of the extensive literature published in the past few years on this
subject, the problem is not completely understood yet. Among the different
competing mechanisms, thermal diffusion becomes one of the most relevant
at large shaking amplitude where the sample of macroscopic grains resembles
a granular gas. In this regime, binary collisions prevail and kinetic theory can
be quite useful to analyze the physical mechanisms involved in segregation
processes.
Thermal diffusion is caused by the relative motion of the components of a
mixture because of the presence of a temperature gradient. Due to this motion,
concentration gradients subsequently appear in the mixture producing diffu-
sion that tends to oppose those gradients. A steady state is finally achieved
in which the separation effect arising from thermal diffusion is compensated
by the diffusion effect. In these conditions, the so-called thermal diffusion fac-
tor Λij characterizes the amount of segregation parallel to the temperature
gradient. In this Section, the thermal diffusion factor is determined from the
Chapman–Enskog solution described before.
To make some contact with experiments, let us assume that the binary
granular mixture is in the presence of the gravitational field g = −gêz, where g
is a positive constant and êz is the unit vector in the positive direction of the z
axis. In experiments [41], the energy is usually supplied by vibrating horizontal
walls so that the system reaches a steady state. Here, instead of considering
oscillating boundary conditions, particles are assumed to be heated by the
action of the stochastic driving force (30), which mimics a thermal bath.
As said above, although the relation between this driven idealized method
with the use of locally driven wall forces is not completely understood, it
Kinetic Theory for Binary Granular Mixtures at Low-Density 31
must be remarked that in the case of boundary conditions corresponding to
a sawtooth vibration of one wall the condition to determine the temperature
ratio coincides with the one derived from the stochastic force [34]. The good
agreement between theory and simulation found in Fig. 1 for the temperature
ratio confirms this expectation.
The thermal diffusion factor Λij (i 6= j) is defined at the steady state in
which the mass fluxes ji vanish. Under these conditions, the factor Λij is given
through the relation [62]
− Λij∇ lnT =
∇xi, Λij + Λji = 0. (111)
The physical meaning of Λij can be described by considering a granular binary
mixture held between plates at different temperatures T (top plate) and T ′
(bottom plate) under gravity. For the sake of concreteness, we will assume
that gravity and thermal gradient point in parallel directions, i.e., the bottom
is hotter than the top (T ′ > T ). In addition, without loss of generality, we
also assume that σ1 > σ2. In the steady state, Eq. (111) describes how the
thermal field is related to the composition of the mixture. Assuming that Λ12 is
constant over the relevant ranges of temperature and composition, integration
of Eq. (111) yields
= Λ12 ln
, (112)
where xi refers to the mole fraction of species i at the top plate and x
i refers
to the mole fraction of species i at the bottom plate. Consequently, according
to Eq. (112), if Λ12 > 0, then x
1 < x1, while if Λ12 < 0, then x
1 > x1. In
summary, when Λ12 > 0, the larger particles accumulate at the top of the
sample (cold plate), while if Λ12 < 0, the larger particles accumulate at the
bottom of the sample (hot plate). The former situation is referred to as the
Brazil-nut effect (BNE) while the latter is called the reverse Brazil-nut effect
(RBNE).
The RBNE was first observed by Hong et al. [63] in MD simulations of
vertically vibrated systems. They proposed a very simple segregation criterion
that was later confirmed by Jenkins and Yoon [64] by using kinetic theory.
More recently, Breu et al. [65] have experimentally investigated conditions
under which the large particles sink to the bottom and claim that their ex-
periments confirm the theory of Hong et al. [63] provided a number of condi-
tions are chosen carefully. In addition to the vertically vibrated systems, some
works have also focused in the last few years on horizontally driven systems
showing some similarities to the BNE and its reverse form [66]. However, it
is important to note that the criterion given in Ref. [63] is based on some
drastic assumptions: elastic particles, homogeneous temperature, and energy
equipartition. These conditions preclude a comparison of the kinetic theory
derived here with the above simulations.
Some theoretical attempts to assess the influence of non-equipartition on
segregation have been recently published. Thus, Trujillo et al. [67] have derived
32 Vicente Garzó
0.5 0.6 0.7 0.8 0.9 1.0
Fig. 8. Phase diagram for BNE/RBNE for mixtures constituted by spheres (d = 3)
of the same mass density and equal total volumes of large and small particles. The
data points represent the MD simulation results [41] for α = 0.78 when convection
is suppressed. Points below (above) the curve correspond to RBNE (BNE).
an evolution equation for the relative velocity of the intruders starting from
the kinetic theory proposed by Jenkins and Yoon [64], which applies for weak
dissipation. They use constitutive relations for partial pressures that take into
account the breakdown of energy equipartition between the two species. How-
ever, the influence of temperature gradients, which exist in the vibro-fluidized
regime, is neglected in Ref. [67] because it is assumed that the pressure and
temperature are constant in the absence of the intruder. A more refined theory
has recently been provided by Brey et al. [42] in the case of a single intruder
in a vibrated granular mixture under gravity. The theory displayed in this
section covers some of the aspects not accounted for in the previous theories
[42, 64, 67] since it is based on a kinetic theory [19] that goes beyond the
quasi-elastic limit [64, 67] and applies for arbitrary composition x1 (and so,
it reduces to the results obtained in Ref. [42] when x1 → 0). This allows one
to assess the influence of composition and dissipation on thermal diffusion in
bi-disperse granular gases without any restriction on the parameter space of
the system.
To determine the dependence of the coefficient Λ12 on the parameters of
the mixture, we consider a non-convecting (u = 0) steady state with only
gradients along the vertical direction (z axis). In this case, the mass balance
equation (12) yields j1 = j2 = 0, while the momentum equation (13) gives
= −ρg. (113)
To first order in the spatial gradients, the constitutive equation for the mass
flux j1,z is given by Eq. (35), i.e.,
Kinetic Theory for Binary Granular Mixtures at Low-Density 33
j1,z = −
m1m2n
, (114)
where the susceptibility coefficient χij = 0 in the particular case of the grav-
itational force. The condition j1,z = 0 yields
m1m2np
m1m2p
, (115)
where use has been made of Eq. (113). Substitution of Eq. (115) into Eq.
(111) leads to
Λ12 =
D′ −Dpg
, (116)
where
n (∂T/∂z)
< 0 (117)
is the reduced gravity acceleration. Since the mutual diffusion coefficient D
is positive [19, 44], the sign of Λ12 is determined by the sign of the quantity
D′ −Dpg
∗. This result is general since it goes beyond the regime of density
considered.
1 2 3 4 5
Fig. 9. Phase diagram for BNE/RBNE in three dimensions for αij = 0.7 and three
values of composition: (a) x1 = 0, (b) x1 = 0.3, and (c) x1 = 0.7. Points below
(above) each curve correspond to RBNE (BNE).
To gain some insight into the explicit dependence of D′ and Dp on the
parameter space of the system, one has to resort to a kinetic theory descrip-
tion. For a low-density gas, the expressions of the coefficients D′ and Dp in
the first-Sonine approximation are given by
34 Vicente Garzó
1 2 3 4
Fig. 10. Phase diagram for BNE/RBNE in three dimensions for x1 =
and three
values of the (common) coefficient of restitution: (a) α = 0.9, (b) α = 0.8, and (c)
α = 0.5. Points below (above) each curve correspond to RBNE (BNE).
D′ = 0, Dp =
x2 + x1γ
, (118)
where ν is the (positive) collision frequency [44]
2π(d−1)/2
) nσd−112 v0(1 + α12)
θ1 + θ2
(x2µ21 + x1µ12) . (119)
Given that the driving stochastic term does not play a neutral role in the
transport, it must be remarked that the expressions for the transport coeffi-
cients obtained in the driven case slightly differ from the ones derived in the
free cooling case [19, 44].
Consequently, according to Eqs. (116) and (118), the sign of Λ12 is the
same as that of the pressure diffusion coefficient Dp. The condition Λ12 = 0
(or equivalently, Dp = 0) provides the criterion for the transition from BNE
to RBNE. Equation (118) shows that the sign of Dp is determined by the
value of the control parameter
. (120)
This parameter gives the mean square velocity of the large particles relative
to that of the small particles. Thus, when θ > 1 (θ < 1), the thermal diffusion
factor is positive (negative), which leads to BNE (RBNE). The criterion for
the transition condition from BNE to RBNE is θ = 1, i.e.,
. (121)
Kinetic Theory for Binary Granular Mixtures at Low-Density 35
In the case of equal granular temperatures (energy equipartition), θ → µ−1
and so segregation is predicted for particles that differ in mass, no matter
what their diameters may be [64]. It must be remarked that, due to the lack
of energy equipartition, the condition θ = 1 is rather complicated since it
involves all the parameter space of the system. In particular, even when the
species differ only by their respective coefficients of restitution they also seg-
regate when subject to a temperature gradient. This is a novel pure effect of
inelasticity on segregation [25, 68]. On the other hand, the criterion (121) for
the transition BNE⇐⇒RBNE is the same as the one found previously in Ref.
[67] when αij is close to 1 and in Ref. [42] in the intruder limit case (x1 → 0).
However, as said before, the results obtained here are more general since they
cover all the range of the parameter space of the system.
To illustrate size segregation driven by thermal diffusion, we consider mix-
tures constituted by spheres (d = 3) of the same material and equal total
volumes of large and small particles. In this case, m1/m2 = (σ1/σ2)
3 and
x2/x1 = (σ1/σ2)
3. Figure 8 shows the phase diagram BNE/RBNE for this
kind of systems. The data points represent the simulation results obtained by
Schröter et al. [41] for α = 0.78 in agitated mixtures constituted by parti-
cles of the same density. To the best of my knowledge, this is one of the few
experiments in which thermal diffusion has been isolated from the remain-
ing segregation mechanisms [61]. Our results show that, for a given value of
the coefficient of restitution, the RBNE is dominant at small diameter ra-
tios. However, since non-equipartition grows with increasing diameter ratio,
the system shows a crossover to BNE at sufficiently large diameter ratios.
This behavior agrees qualitatively well with the results reported in Ref. [41]
at large shaking amplitudes, where thermal diffusion becomes the relevant
segregation mechanism. At a quantitative level, we observe that the results
are also consistent with the simulation results reported in [41] when periodic
boundary conditions are used to suppress convection since they do not ob-
serve a change back to BNE for diameter ratios up to 3 (see red squares in
Fig. 11 of [41]). Although the parameter range explored in MD simulations is
smaller than the one analyzed here, one is tempted to extrapolate the simula-
tion data presented in Ref. [41] to roughly predict the transition value of the
diameter ratio at α = 0.78 (which is the value of the coefficient of restitution
considered in the simulations). Thus, if one extrapolates from the simulation
data at the diameter ratios of σ1/σ2 = 2 and σ1/σ2 = 3, one sees that the
transition from RBNE to BNE might be around σ1/σ2 = 10, which would
quantitatively agree with the results reported in Fig. 8. Figure 8 also shows
that the BNE is completely destroyed in the quasielastic limit (α ≃ 1).
Let us now investigate the influence of composition on segregation. Figure 9
shows a typical phase diagram in the three-dimensional case for αij ≡ α = 0.7
and three different values of the mole fraction x1. The lines separate the
regimes between BNE and RBNE. We observe that the composition of the
mixture has significant effects in reducing the BNE as the concentration of
larger particles increases. In addition, for a given value of composition, the
36 Vicente Garzó
transition from BNE to RBNE may occur following two paths: (i) along the
constant mass ratio m2/m1 with increasing size ratio σ1/σ2, and (ii) along
the constant size ratio with increasing mass ratio m2/m1. The influence of
dissipation on the phase diagrams BNE/RBNE is illustrated in Fig. 10 for
d = 3 in the case of an equimolar mixture (x1 =
) and three values of
the (common) coefficient of restitution α. We observe that the role played by
inelasticity is quite important since the regime of RBNE increases significantly
with dissipation. Similar results are found for other values of composition.
In summary, thermal diffusion (which is the relevant segregation mech-
anism in agitated granular mixtures at large shaking amplitudes) can been
analyzed by the Boltzmann kinetic theory. This theory is able to explain
some of the experimental and/or MD segregation results [41] observed within
the range of parameter space explored. A more quantitative comparison in
the dilute regime with MD simulations is needed to show the relevance of
the Boltzmann equation to analyze segregation driven by a thermal gradi-
ent. As said before, comparison with MD simulations in the tracer limit case
(x1 → 0) [42] for a dilute gas has shown the reliability of the inelastic Boltz-
mann equation to describe segregation. In this context, one expects that the
same agreement observed before in the intruder case [42] is maintained when
x1 is different from zero.
10 Steady States: Uniform Shear Flow
In the preceding sections, the Navier–Stokes equations (constitutive equa-
tions that are linear in the hydrodynamic gradients) have been shown to be
quite useful to describe appropriately several problems in granular mixtures.
However, under some circumstances large gradients occur and more complex
constitutive equations are required. The need for more complex constitutive
equations does not signal a breakdown of hydrodynamics [69], only a failure
of the Navier–Stokes approximation [70]. Although in this case the Chapman-
Enskog method can be carried out to second order in gradients (Burnett or-
der), it is likely that failure of the Navier-Stokes description signals the need
for other methods to construct the normal solution that are not based on a
small gradient expansion.
One of the most interesting problems in granular fluids is the simple or
uniform shear flow (USF) [13, 26]. As described in Sec. 5, this state is charac-
terized by uniform density and temperature and a simple shear with the local
velocity field given by u1,x = u2,x = ay, uy = uz = 0, where a is the constant
shear rate. The USF is a well-known nonequilibrium problem widely studied,
for both granular [26, 27, 71, 72, 73] and conventional [28, 74] gases. However,
the nature of this state is quite different in both systems since a steady state
is achieved for granular fluids when viscous (shear) heating is compensated
for by energy dissipation in collisions:
Kinetic Theory for Binary Granular Mixtures at Low-Density 37
aPxy = −
nTζ. (122)
This steady state is what we want to analyze in this section. The balance
equation (122) shows the intrinsic connection between the shear field and dis-
sipation in the system. This contrasts with the description of USF for elastic
fluids where a steady state is not possible unless an external thermostat is in-
troduced [28]. Note that the hydrodynamic steady shear flow state associated
with the condition (122) is inherently beyond the scope of the Navier–Stokes
or Newtonian hydrodynamic equations [29]. The reason for this is the exis-
tence of an internal mechanism, collisional cooling, that sets the strength of
the velocity gradient in the steady state. For normal fluids, this scale is set by
external sources (boundary conditions, driving forces) that can be controlled
to admit the conditions required for Navier–Stokes hydrodynamics. In con-
trast, collisional cooling is fixed by the mechanical properties of the particles
making up the fluid. This observation is significant because it prevents the
possibility of measuring the Newtonian shear viscosity for granular fluids in
the steady USF [72, 73]. More generally, it provides a caution regarding the
simulation of other steady states to study Navier–Stokes hydrodynamics when
the gradients are strongly correlated to the collisional cooling [29].
From a microscopic point of view, the simple shear flow problem becomes
spatially uniform in the local Lagrangian frame moving with the flow velocity
u. In this frame [28, 75, 76], the velocity distribution functions adopt the
form: fi(r,v) → fi(V), where Vk = vk − akℓrℓ is the peculiar velocity. Here,
akℓ = aδkxδℓy. Under these conditions, the set of Boltzmann kinetic equations
(with Fi = 0) for an isolated system reads
− aVy
fi(V) =
Jij [V|fi, fj] , (i = 1, 2). (123)
The most relevant transport properties in a shear flow problem are obtained
from the pressure tensor P = P1 + P2, where Pi is the partial pressure tensor
of the species i given by
Pi,kℓ = mi
dVVkVℓfi(V). (124)
The trace of Pi defines the partial temperatures Ti as Ti = TrPi/dni.
As said before, these temperatures measure the mean kinetic energy of each
species. The elements of the pressure tensor Pi can be obtained by multiplying
the Boltzmann equation (123) by miVV and integrating over V. The result
akmPi,mℓ + aℓmPi,mk =
Aij,kℓ, (125)
where we have introduced the collisional moments Aij as
38 Vicente Garzó
Aij,kℓ = mi
dVVkVℓJij [V|fi, fj ]. (126)
From Eq. (125), in particular, one gets the balance equation for the partial
temperature Ti
aPi,xy = −
piζi, (127)
where pi = niTi is the partial pressure of species i and ζi is defined by Eq.
(10). According to Eq. (127), the (steady) partial temperature in the simple
shear flow problem can be obtained by equating the viscous heating term
a|Pi,xy| to the collisional cooling term (d/2)piζi.
The determination of Aij requires the knowledge of the velocity distribu-
tion functions fi. This is quite a formidable task, even in the monocomponent
case [27]. However, as in the elastic case, one expects to get a good estimate
of Aij by using Grad’s approximation [43]:
f(V) → fi,M (V)
Ci,kℓVkVℓ
, (128)
where fi,M is a Maxwellian distribution at the temperature of the species i,
i.e.,
fi,M (V) = ni
. (129)
As happens in the case of homogeneous states, in general the three temper-
atures T , T1, and T2 are different in the inelastic case. For this reason we
choose the parameters in the Maxwellians so that it is normalized to ni and
provides the exact second moment of fi. The Maxwellians fi,M for the two
species can be quite different due to the temperature differences. This aspect
is essential in a two-temperature theory and has not been taken into account
in most of the previous studies [14, 71, 72]. The coefficient Ci can be identified
by requiring the moments with respect to VV of the trial function (128) to
be the same as those for the exact distribution fi. This leads to
− I (130)
With this approximation, the Boltzmann collisional moments Aij can be
explicitly evaluated. The result is [77, 78]
Aij = −
2π(d−1)/2
dΓ (d/2)
mininjµjiσ
(1 + αij)
Tj − Ti
(mj/mi)Ti + Tj
1− αij
1 + (miTj/mjTi)
Ci − Cj
1 + (mjTi/miTj)
2(d+ 2)
, (131)
Kinetic Theory for Binary Granular Mixtures at Low-Density 39
where
λij = 2µji
Tj − Ti
(mj/mi)Ti + Tj
(2d+ 3− 3αij). (132)
The partial cooling rates ζi can be easily obtained from Eqs. (10) and (131).
1.0 1.5 2.0 2.5 3.0 3.5 4.0
α=0.8
α=0.9
Fig. 11. Plot of the temperature ratio T1/T2 as a function of the size ratio σ1/σ2 =
(m1/m2)
1/2 for a two-dimensional system in the case x1 = 1/2 and two different
values of the (common) coefficient of restitution: α = 0.9 and α = 0.8. The solid
lines are the theoretical predictions based on Grad’s solution, while the symbols
refer to the DSMC results. The dashed lines correspond to the results obtained from
the stochastic thermostat condition (31).
Substitution of Eq. (131) into the set of equations (125) allows one to get
the partial pressure tensor Pi in terms of the temperature ratio γ = T1/T2
and the parameters of the mixture. The temperature ratio can be obtained
from Eq. (127) as
x2ζ2P1,xy
x1ζ1P2,xy
. (133)
When the expressions of Pi and ζi are used in Eq. (133), one gets a closed
equation for the temperature ratio γ, that can be solved numerically. In Fig.
11 we plot γ versus the diameter ratio σ1/σ2 for a two-dimensional (d =
2) granular gas with x1 = 1/2 and two different values of α. The symbols
refer to the simulation data obtained from the DSMC method [79]. Here,
we have assumed that the disks are made of the same material, and hence
αij = α and m1/m2 = (σ1/σ2)
2. The dependence of γ on σ1/σ2 obtained in
the homogeneous steady state driven by the stochastic thermostat (30) is also
included for comparison. It is clearly seen that the kinetic theory results based
40 Vicente Garzó
on Grad’s solution agree very well with simulation data, even for quite large
values of the size ratio. In addition, the thermostat results overestimate the
simulation ones (especially for large mass ratio), showing that the properties
of the system are not insensitive to the way at which the granular gas is driven.
� �� ���
Fig. 12. Plot of the reduced pressure p∗ versus the mass ratio µ = m1/m2 for a two-
dimensional system with σ1 = σ2, x1 = 1/2 and α = 0.9. The solid line corresponds
to the theoretical predictions derived from Grad’s solution, the dotted line refers
to the latter theory but using the expression of T1/T2 obtained from the stochastic
thermostat condition (31), and the dashed line is the result obtained from Grad’s
solution by assuming the equality of the partial temperatures (γ = 1). The symbols
are the DSMC results.
Let us now consider the transport coefficients. To analyze the rheologi-
cal properties in the steady state, it is convenient to introduce dimensionless
quantities. As usual [72], for a low-density gas we introduce the reduced pres-
sure p∗ and the reduced shear viscosity η∗ as
, (134)
, (135)
where η = −Pxy/a is the non-Newtonian shear viscosity, Pxy = P1,xy + P2,xy
and ν = [π(d−1)/2/Γ (d/2)]nσd−112 v0. In Figs. 12 and 13, we plot p
∗ and η∗,
respectively, as functions of the mass ratio µ = m1/m2 for an equal-size
(σ1 = σ2) binary mixture of disks (d = 2) with x1 = 1/2 and α = 0.9.
Kinetic Theory for Binary Granular Mixtures at Low-Density 41
� �� ���
Fig. 13. Plot of the reduced shear viscosity η∗ versus the mass ratio µ = m1/m2
for a two-dimensional system with σ1 = σ2, x1 = 1/2 and α = 0.9. The solid line
corresponds to the theoretical predictions derived from Grad’s solution, the dotted
line refers to the latter theory but using the expression of T1/T2 obtained from the
stochastic thermostat condition (31), and the dashed line is the result obtained from
Grad’s solution by assuming the equality of the partial temperatures (γ = 1). The
symbols are the DSMC results.
We have also included the predictions for p∗ and η∗ given by the kinetic
theory but taking the expression of γ derived when the system is driven by
the stochastic thermostat [40]. We observe again in both figures an excellent
agreement between the Boltzmann theory based on Grad’s solution and the
DSMC results, even for very disparate values of the mass ratio. With respect
to the influence of energy nonequipartition, Fig. 12 shows that p∗ presents a
non-monotonic behaviour with the mass ratio whereas the theoretical predic-
tions with the equipartition assumption monotonically increase with µ. In the
case of the shear viscosity, as seen in Fig. 13, both theories (with and with-
out energy nonequipartition) predict a non-monotonic dependence of η∗ on
µ. However, at a quantitative level, the influence of energy nonequipartition
is quite significant over the whole range of mass ratios considered. The non-
monotonic dependence of p∗ and η∗ on µ obtained here from the Boltzmann
kinetic theory also agrees qualitatively well with MD simulations carried out
for bidisperse dense systems [72]. Thus, for instance, the minimum values of p∗
and η∗ are located close to µ = 10 in both dilute and dense cases. Moreover,
the predictions for the transport properties given from the present theory
by taking the stochastic thermostat expression of γ are quite close to those
obtained from the actual value of γ, especially for large mass ratios.
42 Vicente Garzó
11 Summary and Concluding Remarks
The primary objective of this review has been to derive the Navier–Stokes hy-
drodynamic equations of a binary mixture of granular gases from the (inelas-
tic) Boltzmann kinetic theory. The Chapman–Enskog method [4, 43] is used
to solve the Boltzmann equation up to the first order in the spatial gradients
and the associated transport coefficients are given in terms of the solutions of
a set of linear integral equations. These equations have been approximately
solved by taking the leading terms in a Sonine polynomial expansion. Compar-
ison with controlled numerical simulations in some idealized conditions shows
quite a good agreement between theory and simulation even for strong dis-
sipation. This supports the idea that the hydrodynamic description (derived
from kinetic theory) appears to be a powerful tool for analysis and predictions
of rapid flow gas dynamics of polydisperse systems [12].
The reference state in the Chapman–Enskog expansion has been taken
to be an exact solution of the uniform Boltzmann equation. An interesting
and important result of this solution [15] is that the partial temperatures
(which measure the mean kinetic energy of each species) are different. This
does not mean that there are additional degrees of freedom since the partial
temperatures can be expressed in terms of the global temperature. This is
confirmed by noting that Haff’s cooling law [1] (in the free cooling case) is the
hydrodynamic mode at long wavelengths and MD simulations confirm that
the global temperature dominates after a transient period of a few collision
times [34]. In this case, only the global temperature should appear among
the hydrodynamic fields. Nevertheless, the species temperatures play a new
and interesting secondary role [20]. For an ordinary (molecular) gas, there is
a rapid velocity relaxation in each fluid cell to a local equilibrium state on
the time scale of a few collisions (e.g., as illustrated by the approach to Haff’s
law). Subsequently, the equilibration among cells occurs via the hydrodynamic
equations. In each cell the species velocity distributions are characterized by
the species temperatures. These are approximately the same due to equiparti-
tion, and the hydrodynamic relaxation occurs for the single common temper-
ature [43]. A similar rapid velocity relaxation occurs for granular gases in each
small cell, but to a universal state different from local equilibrium and one for
which equipartition no longer occurs. Hence, the species temperatures Ti are
different from each other and from the overall temperature T of the cell. Nev-
ertheless, the time dependence of all temperatures (in the free cooling case)
is the same in this and subsequent states, i.e., they are proportional to the
global temperature. This implies that the species temperatures do not provide
any new dynamical degree of freedom at the hydrodynamic stage. However,
they still characterize the shape of the partial velocity distributions and affect
the quantitative averages calculated with these distributions. The transport
coefficients for granular mixtures therefore have new quantitative effects aris-
ing from the time independent temperature ratios for each species [19]. This
view contrasts with some recent works [80], where additional equations for
Kinetic Theory for Binary Granular Mixtures at Low-Density 43
each species temperature have been included among the hydrodynamic set.
However, as mentioned before, this is an unnecessary complication, describing
additional kinetics beyond hydrodynamics that is relevant only on the time
scale of a few collisions.
Another important issue discussed here has been the applicability of the
Navier–Stokes transport coefficients since their expressions are not restricted
to weak inelasticity [12]. However, the Navier–Stokes hydrodynamic equations
themselves may or may not be limited with respect to inelasticity, depending
on the particular states studied. The Chapman–Enskog method assumes that
the relative changes of the hydrodynamic fields over distances of the order
of the mean free path are small. In the case of ordinary fluids this can be
controlled by the initial or boundary conditions. For granular gases the situ-
ation is more complicated since in some cases (e.g., steady states such as the
simple shear flow problem [29]) the boundary conditions imply a relationship
between dissipation and gradients so that both cannot be chosen indepen-
dently. In these cases, the Navier–Stokes approximation only holds for nearly
elastic particles. However, the transport coefficients characterizing the Navier–
Stokes hydrodynamic equations are nonlinear functions of the coefficients of
restitution, regardless the applicability of those equations.
In spite of the above cautions, the Navier–Stokes approximation is ap-
propriate and accurate for a wide class of flows. One group refers to spatial
perturbations of the homogeneous cooling state (HCS) for an isolated sys-
tem. Both MD and DSMC simulations [7] have confirmed the dependence
of the Navier–Stokes transport coefficients on the coefficient of restitution,
and application of the Navier–Stokes hydrodynamics with these coefficients
to describe cluster formation has also been confirmed quantitatively [60]. The
same kinetic theory results apply to driven systems as well. This is so since
the reference state is a local HCS whose parameters vary throughout the sys-
tem to match the physical values in each cell. Examples include application of
Navier–Stokes hydrodynamics from kinetic theory to symmetry breaking and
density/temperature profiles in vertical vibrated gases, for comparison with
simulation [81]. Similar comparison with Navier–Stokes hydrodynamics of the
latter and of supersonic flow past a wedge in real experiments has been given
[82, 83], showing both qualitative and quantitative agreement. In summary,
the Navier–Stokes equations with the constitutive equations presented here
remain an important and useful description for a wide class of granular flow,
although more limited than for normal gases.
The explicit knowledge of the transport coefficients and the cooling rate
allows one to make some applications of the Navier–Stokes hydrodynamic
equations. One of them has been to obtain the linear hydrodynamic equa-
tions for small perturbations of the homogenous cooling state. The resulting
equations exhibit a long wavelength instability for three of the modes. This
is quite similar to the case of a monocomponent granular gas [5, 54, 57], and
in fact the same modes are unstable here. The additional diffusion mode for
two species behaves as for a normal fluid.
44 Vicente Garzó
On the other hand, the constitutive equations for the mass and heat fluxes
of a granular binary mixture differ from those obtained for ordinary fluids [4].
This is because the usual restrictions of irreversible thermodynamics no longer
apply. These restrictions include Onsager’s reciprocal relations among various
transport coefficients and the extent to which these are violated has also been
shown here. Another application of the Navier–Stokes equations has been to
assess the violation of the Einstein relation between the diffusion and mobility
coefficients. In the undriven case, the analysis shows that this violation is
due to three independent reasons [23]: the absence of the Gibbs state, the
cooling of the reference state, and the occurrence of different temperatures for
the particle and surrounding fluid. However, when the mixture is subjected
to stochastic driving, a modified Einstein relation suggested by recent MD
simulations [52] has also been analyzed. In this case, the results show that
the deviations of the (modified) Einstein ratio from unity are in general very
small (less than 1%), in agreement with MD simulations [52].
Thermal diffusion becomes the relevant segregation mechanism in agitated
granular mixtures at large shaking amplitudes. In these conditions, the use of
the Boltzmann kinetic theory for low-density gases appears justified to under-
stand the influence of thermal gradient on segregation phenomena. The ther-
mal diffusion factor in a heated granular mixture has been explicitly evaluated
from the Chapman–Enskog solution to the Boltzmann equation. The results
show that the criterion for the transition Brazil-nut effect ⇐⇒reverse Brazil-
nut effect is provided by the control parameter θ = m2T1/m1T2 [25, 42, 67].
Given that the energy equipartition is broken, the condition θ = 1 is quite
complex since it involves all the parameters of the system: composition,
masses, sizes, and coefficients of restitution. The Boltzmann kinetic theory
results agree qualitatively well with recent MD simulations [41] within the
range of parameter space analyzed.
The hydrodynamic description also seems to be justified in the case of
steady states that are inherently beyond the scope of the Navier–Stokes hy-
drodynamic equations. The reason for this non-Newtonian behavior is the
existence of an internal mechanism, collisional cooling, that sets the scale of
the spatial gradients in the steady state. For ordinary fluids, this scale can be
externally controlled by external sources so that the conditions for Navier–
Stokes hydrodynamics apply. On the other hand, for granular gases, collisional
cooling is fixed by the mechanical properties of the particles of the system and
so the gas can depart from the Navier–Stokes description. One well-known ex-
ample of steady states is the simple or uniform shear flow (USF). However, in
spite of the extensive prior work on USF for granular fluids [27, 72, 73], the in-
herent non-Newtonian character of this state has not been conveniently taken
into account. In fact, MD simulations of steady USF have been used for gran-
ular fluids to measure the Newtonian or Navier–Stokes shear viscosity. The
results derived here from Grad’s solution and DSMC simulations show that
USF is an ideal testing ground for the study of rheology since any choice of the
shear rate and the coefficients of restitution αij will provide non-Newtonian
Kinetic Theory for Binary Granular Mixtures at Low-Density 45
effects. It is one of the fascinating features of granular fluids that phenom-
ena associated with complex fluids are more easily accessible than for simple
atomic fluids [12, 84].
Hydrodynamics derived from hard-sphere models have found widespread
use in the description of numerous industrial processes involving solid parti-
cles. Of particular relevance are high-speed, gas-solid flows as found in pneu-
matic conveyors (of ores, chemicals, grains, etc.) and fluidized beds (for fluid
catalytic cracking, power generation, granulation of pharmaceutical powders,
synthesis of fine chemicals like titania, etc.). Such descriptions are now stan-
dard features of commercial and research codes. Those codes rely upon accu-
rate transport properties and a first order objective is to assure this accuracy
from a careful theoretical treatment. As shown in this review, the price of this
approach, in contrast to more phenomenological approaches, is an increasing
complexity of the expressions as the systems become more complex.
The analysis carried out in this presentation has been focused on mix-
tures in the dilute regime, where the collisional transfer contributions to the
transport coefficients are neglected and only their kinetic contributions are
considered. A further step is to develop a theory for moderately dense gran-
ular mixtures. This will provide a fundamental basis for the application of
hydrodynamics under realistic conditions. Possible extension of the present
Boltzmann kinetic theory to higher densities can be done in the context of
the Enskog kinetic equation [43]. Preliminary results [85] have been restricted
to the uniform shear flow state to get directly the shear viscosity coefficient of
a heated granular mixture. The extension of this study [85] to states with gra-
dients of concentration, pressure, and temperature is somewhat intricate due
to subtleties associated with the spatial dependence of the pair correlations
functions considered in the revised Enskog theory. A future work is to ex-
tend the results derived for moderately dense mixtures of smooth elastic hard
spheres [86] to inelastic collisions. This would allow us to assess the influence
of density on the different problems addressed in this review. Of course, the
precise expressions for transport coefficients in this case will be even more
complex than for a dilute gas due to the expanded parameter space. However,
this complexity is not a problem for implementation in a code.
As shown along this overview, granular mixtures exhibit a wide range of
interesting phenomena for which the Navier–Stokes hydrodynamic equations
can be considered as an accurate and practical tool. However, due to their
complexity, many of their features are not fully understood. Kinetic theory
and hydrodynamics (in the broader sense) can be expected to provide some
insight into the understanding of such complex materials.
Acknowledgments
I want to acknowledge J. W. Dufty, J. M. Montanero, and A. Santos in
their roles as collaborators and critics for much of the material discussed here.
46 Vicente Garzó
Partial support of the Ministerio de Ciencia y Tecnoloǵıa (Spain) through
Grant No. FIS2004-01399 (partially financed by FEDER funds) and from the
European Community’s Human Potential Programme HPRN-CT-2002-00307
(DYGLAGEMEM) is also acknowledged.
References
1. P. K. Haff, J. Fluid Mech. 134, 187 (1983).
2. L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999).
3. See for instance, J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids
(Academic Press, London, 1986).
4. S. Chapman and T. G. Cowling, The Mathematical Theory of Nonuniform Gases
(Cambridge University Press, Cambridge, 1970).
5. J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998).
6. J. J. Brey and D. Cubero, in Granular Gases, edited by T. Pöschel and S. Luding
(Lectures Notes in Physics, Springer, Vol. 564, 2001), pp. 59–78.
7. J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Europhys. Lett. 48, 359 (1999).
8. J. J. Brey, M. J. Ruiz-Montero, D. Cubero, and R. Garćıa-Rojo, Phys. Fluids
12, 876 (2000).
9. J. J. Brey and M. J. Ruiz-Montero, Phys. Rev. E 70, 051301 (2004).
10. J. M. Montanero, A. Santos, and V. Garzó, in Rarefied Gas Dynamics 24, edited
by M. Capitelli (American Institute of Physics, Vol. 762, 2005), pp. 797–802;
preprint cond-mat/0411219.
11. J. W. Dufty, J. Phys.: Condens. Matt. 12, A47 (2000).
12. J. W. Dufty, in Recent Research Develop. in Stat. Phys. 2, 21 (2002) and preprint
cond-mat/0108444.
13. I. Goldhirsch, Annu. Rev. Fluid Mech. 35, 267 (2003).
14. J. T. Jenkins and F. Mancini, Phys. Fluids A 1, 2050 (1989); P. Zamankhan,
Phys. Rev. E 52, 4877 (1995); B. Arnarson and J. T. Willits, Phys. Fluids 10,
1324 (1998); J. T. Willits and B. Arnarson, ibid. 11, 3116 (1999); M. Alam, J.
T. Willits, B. Arnarson, and S. Luding, ibid. 14, 4085 (2002); B. Arnarson and
J. T. Jenkins, ibid. 16, 4543 (2004); D. Serero, I. Goldhirsch, S. H. Noskowicz,
and M.-L. Tan, J. Fluid Mech. 554, 237 (2006).
15. V. Garzó and J. W. Dufty, Phys. Rev. E 60, 5706 (1999).
16. See for instance, J. M. Montanero and V. Garzó, Gran. Matt. 4, 17 (2002); A.
Barrat and E. Trizac, ibid. 4, 57 (2002); S. R. Dahl, C. M. Hrenya, V. Garzó, and
J. W. Dufty, Phys. Rev. E 66, 041301 (2002); R. Pagnani, U. M. B. Marconi,
and A. Puglisi, Gran. Matt. 66, 051304 (2002); D. Paolotti, C. Cattuto, U. M.
B. Marconi, and A. Puglisi, ibid. 5, 75 (2003); P. Krouskop and J. Talbot, Phys.
Rev. E 68, 021304 (2003); H. Wang, G. Jin, and Y. Ma, ibid. 68, 031301 (2003);
J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, Phys. Rev. Lett. 95, 098001
(2005); Phys. Rev. E 73, 031301 (2006); M. Schröter, S. Ulrich, J. Kreft, J. B.
Swift, and H. L. Swinney, Phys. Rev. E 74, 011307 (2006).
17. R. D. Wildman and D. J. Parker, Phys. Rev. Lett. 88, 064301 (2002); K. Feitosa
and N. Menon, ibid. 88, 198301 (2002); M. Schröter, S. Ulrich, J. Kreft, J. B.
Swift, and H. L. Swinney, Phys. Rev. E 74, 011307 (2006).
18. J. Jenkins and F. Mancini, J. Appl. Mech. 54, 27 (1987).
http://arxiv.org/abs/cond-mat/0411219
http://arxiv.org/abs/cond-mat/0108444
Kinetic Theory for Binary Granular Mixtures at Low-Density 47
19. V. Garzó and J. W. Dufty, Phys. Fluids 14, 1476 (2002).
20. V. Garzó, J. M. Montanero, and J. W. Dufty, Phys. Fluids 18, 083305 (2006).
21. J. M. Montanero and V. Garzó, Phys. Rev. E 67, 021308 (2003).
22. V. Garzó and J. M. Montanero, Phys. Rev. E 69, 021301 (2004).
23. J. W. Dufty and V. Garzó, J. Stat. Phys. 105, 723 (2001).
24. V. Garzó, Physica A 343, 105 (2004).
25. V. Garzó, Europhys. Lett. 75, 521 (2006).
26. C. S. Campbell, Annu. Rev. Fluid Mech. 22, 57 (1990).
27. In the case of a monocomponent system, see for instance, C. K. K. Lun, S. B.
Savage, D. J. Jeffrey, and N. Chepurniy, J. Fluid Mech. 140, 223 (1984); J. T.
Jenkins and M. W. Richman, ibid. 192, 313 (1988); C. S. Campbell, ibid. 203,
449 (1989); M. A. Hopkins and H. H. Shen, ibid. 244, 477 (1992); P. J. Schmid
and H. K. Kytömaa, ibid. 264, 255 (1994); C. K. K. Lun and A. A. Bent, ibid.
258, 335 (1994); I. Goldhirsch and M. L. Tan, Phys. Fluids 8, 1752 (1996); N.
Sela, I. Goldhirsch, and S. H. Noskowicz, ibid. 8, 2337 (1997); J. J. Brey, M. J.
Ruiz-Montero, and F. Moreno, Phys. Rev. E 55, 2846 (1997); C.-S. Chou and
M. W. Richman, Physica A 259, 430 (1998); C.-S. Chou ibid. 287, 127 (2000);
290, 341 (2001); J. M. Montanero, V. Garzó, A. Santos, and J. J. Brey, J. Fluid
Mech. 389, 391 (1999); A. Astillero and A. Santos, Phys. Rev. E 72, 031309
(2005).
28. V. Garzó and A. Santos, Kinetic Theory of Gases in Shear Flows. Nonlinear
Transport (Kluwer Academic, Dordrecht, 2003).
29. A. Santos, V. Garzó, and J. W. Dufty, Phys. Rev. E 69, 061303 (2004).
30. A. Goldshtein and M. Shapiro, J. Fluid Mech. 282, 75 (1995).
31. J. J. Brey, J. W. Dufty, and A. Santos, J. Stat. Phys. 87, 1051 (1997).
32. T. P. C. van Noije and M. H. Ernst, Gran. Matt. 1, 57 (1998).
33. J. M. Montanero and V. Garzó, Gran. Matt. 4, 17 (2002).
34. S. R. Dahl, C. M. Hrenya, V. Garzó, and J. W. Dufty, Phys. Rev. E 66, 041301
(2002).
35. D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids
(Academic Press, London, 1990).
36. W. G. Hoover, Computational Statistical Mechanics (Elsevier, Amsterdam,
1991).
37. J. M. Montanero and A. Santos, Gran. Matt. 2, 53 (2000).
38. D. R. M. Williams and F. C. McKintosh, Phys. Rev. E 54, R9 (1996).
39. C. Henrique, G. Batrouni, and D. Bideau, Phys. Rev. E 63, 011304 (2000).
40. A. Barrat and E. Trizac, Gran. Matt. 4, 52 (2002).
41. M. Schröter, S. Ulrich, J. Kreft, J. B. Swift, and H. L. Swinney, Phys. Rev. E
74, 011307 (2006).
42. J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, Phys. Rev. Lett. 95, 098001
(2005); Phys. Rev. E 73, 031301 (2006).
43. J. Ferziger and H. Kaper, Mathematical Theory of Transport Processes in Gases
(North-Holland, Amsterdam, 1972).
44. V. Garzó and J. M. Montanero, preprint cond-mat/0604552.
45. V. Garzó and J.W. Dufty, Phys. Rev. E 59, 5895 (1999).
46. J. F. Lutsko, Phys. Rev. E 72, 021306 (2005).
47. G. A. Bird, Molecular Gas Dynamics and the Direct Simulation Monte Carlo of
Gas Flows (Clarendon, Oxford, 1994).
48. J. A. McLennan, Introduction to Nonequilibrium Statistical Mechanics (Prentice-
Hall, New Jersey, 1989).
http://arxiv.org/abs/cond-mat/0604552
48 Vicente Garzó
49. A. Santos and J. W. Dufty, Phys. Rev. Lett. 97, 058001 (2006).
50. M. López de Haro and E. G. D. Cohen, J. Chem. Phys. 80, 408 (1984).
51. V. Garzó and J. M. Montanero, Physica A 313, 336 (2002).
52. A. Barrat, V. Loreto, and A. Puglisi, Physica A 334, 513 (2004).
53. S. R. de Groot and P. Mazur, Nonequilibrium Thermodynamics (Dover, New
York, 1984).
54. N. Brilliantov and T. Pöschel, Kinetic Theory of Granular Gases (Oxford Uni-
versity Press, Oxford, 2004).
55. I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993); I. Goldhirsch,
M.L. Tan, and G. Zanetti, J. Sci. Comput. 8, 1 (1993).
56. S. Luding and H.J. Herrmann, Chaos 9, 673 (1999).
57. V. Garzó, Phys. Rev. E 72, 021106 (2005).
58. P. Résibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley,
NY, 1977).
59. J. P. Boon and S. Yip, Molecular Hydrodynamics (Dover, NY, 1980).
60. J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 60, 3150 (1999).
61. See, for instance, A. Kudrolli, Rep. Prog. Phys. 67, 209 (2004).
62. J. M. Kincaid, E. G. D. Cohen, and M. López de Haro, J. Chem. Phys. 86, 963
(1983).
63. D. C. Hong, P. V. Quinn, and S. Luding, Phys. Rev. Lett. 86, 3423 (2001).
64. J. T. Jenkins and D. K. Yoon, Phys. Rev. Lett. 88, 194301 (2002); D. K. Yoon
and J. T. Jenkins, Phys. Fluids 18, 073303 (2006).
65. A. P. J. Breu, H. M. Ensner, C. A. Kruelle, and I. Rehberg, Phys. Rev. Lett.
90, 014302 (2003).
66. T. Schautz, R. Brito, C. A. Kruelle, and I. Rehberg, Phys. Rev. Lett. 95, 028001
(2005) and references therein.
67. L. Trujillo, M. Alam, and H. J. Herrmann, Europhys. Lett. 64, 190 (2003).
68. D. Serero, I. Goldhirsch, S. H. Noskowicz, and M.-L. Tan, J. Fluid Mech. 554,
237 (2006).
69. M. L. Tan and I. Goldhirsch, Phys. Rev. Lett. 81, 3022 (1998).
70. J. W. Dufty and J. J. Brey, Phys. Rev. Lett. 82, 4566 (1999).
71. J. T. Willits and B. Arnarson, Phys. Fluids 11, 3116 (1999); M. Alam, J. Willits,
B. Arnarson, and S. Luding, ibid. 14, 4085 (2002).
72. M. Alam and S. Luding, Gran. Matt. 4, 139 (2002); J. Fluid Mech. 476, 69
(2003); Phys. Fluids 17, 063303 (2005).
73. M. Alam and C. M. Hrenya, Phys. Rev. E 63, 0613018 (2001); R. Clelland
and C. M. Hrenya, ibid. 65, 031301 (2002); S. R. Dahl, R. Clelland, and C.
M. Hrenya, Phys. Fluids 14, 1972 (2002); S. R. Dahl, R. Clelland, and C. M.
Hrenya, Powder Tech., 138, 7 (2003); S. R. Dahl and C. M. Hrenya, Phys. Fluids
16, 1 (2004); H. Iddir, H. Arastoopour, and C. M. Hrenya, Powder Tech. 151,
117 (2005).
74. See for instance, Nonlinear Fluid Behavior, edited by H. J. M. Hanley (North-
Holland, Amsterdam, 1983).
75. A. W. Lees and S. F. Edwards, J. Phys. C 5, 1921 (1972).
76. J. W. Dufty, A. Santos, J. J. Brey, and R. F. Rodŕıguez, Phys. Rev. A 33, 459
(1986).
77. J. M. Montanero and V. Garzó, Physica A 310, 17 (2002).
78. V. Garzó, Phys. Rev. E 66, 021308 (2002).
79. V. Garzó and J. M. Montanero, Gran. Matt. 5, 165 (2003).
Kinetic Theory for Binary Granular Mixtures at Low-Density 49
80. L. Huilin, D. Gidaspow, and E. Manger, Phys. Rev. E 64, 061301 (2001); M.F.
Ramahan, J. Naser, and P. J. Witt, Powder Tech. 138, 82 (2003); J. E. Galvin,
S. R. Dahl, and C. M. Hrenya, J. Fluid Mech. 528, 207 (2005).
81. J. J. Brey, M. J. Ruiz-Montero, F. Moreno, and R. Garćıa-Rojo, Phys. Rev. E
63, 061305 (2001); ibid. 65, 061302 (2002)
82. X. Yang, C. Huan, D. Candela, R. W. Mair, and R. L. Walsworth, Phys. Rev.
Lett. 88, 044301 (2002); C. Huan, X. Yang, D. Candela, R. W. Mair, and R. L.
Walsworth, Phys. Rev. E 69, 041302 (2004).
83. C. Bizon, M. D. Shattuck, J. B. Swift, and Harry L. Swinney, Phys. Rev. E 60,
4340 (1999); E. C. Rericha, C. Bizon, M. D. Shattuck, and H. L. Swinney, Phys.
Rev. Lett. 88, 014302 (2002).
84. I. Goldhirsch, in Granular Gases, edited by T. Pöschel and S. Luding (Lectures
Notes in Physics, Springer, Vol. 564, 2001), pp. 79–99.
85. V. Garzó and J. M. Montanero, Phys. Rev. E 68, 041302 (2003).
86. M. López de Haro, E. G. D. Cohen, and J. M. Kincaid, J. Chem. Phys. 78, 2746
(1983).
	Kinetic Theory for Binary Granular Mixtures at Low-Density
	Vicente Garzó
