A Explicit expressions in the Sonine approximation

# Transport coefficients of a granular gas of inelastic rough hard spheres

## Abstract

The Boltzmann equation for inelastic and rough hard spheres is considered as a model of a dilute granular gas. In this model, the collisions are characterized by constant coefficients of normal and tangential restitution and hence the translational and rotational degrees of freedom are coupled. A normal solution to the Boltzmann equation is obtained by means of the Chapman–Enskog method for states near the homogeneous cooling state. The analysis is carried out to first order in the spatial gradients of the number density, the flow velocity, and the granular temperature. The constitutive equations for the momentum and heat fluxes and for the cooling rate are derived, and the associated transport coefficients are expressed in terms of the solutions of linear integral equations. For practical purposes, a first Sonine approximation is used to obtain explicit expressions of the transport coefficients as nonlinear functions of both coefficients of restitution and the moment of inertia. Known results for purely smooth inelastic spheres and perfectly elastic and rough spheres are recovered in the appropriate limits.

###### pacs:
45.70.Mg,05.20.Dd,51.10.+y,05.60.-k

## I Introduction

As is well known, the prototypical model of a granular gas is a system composed of smooth, frictionless hard spheres which collide inelastically with a constant coefficient of normal restitution Brilliantov and Pöschel (2004); Campbell (1990); Goldhirsch (2003). In the dilute and moderately dense regimes, the microscopic description of the gas is given by the one-particle velocity distribution function obeying the (inelastic) Boltzmann and Enskog kinetic equations Brilliantov and Pöschel (2004); Goldshtein and Shapiro (1995); Brey et al. (1997). At a more phenomenological level, the gas can also be described by the Navier–Stokes–Fourier (NSF) hydrodynamic equations for the densities of mass, momentum, and energy with appropriate constitutive equations for the stress tensor, heat flux, and cooling rate. The Chapman–Enskog method Chapman and Cowling (1970); Ferziger and Kaper (1972) bridges the gap between the kinetic and hydrodynamic descriptions, thus providing explicit expressions for the NSF transport coefficients in terms of the coefficient of normal restitution. This task was first accomplished in the quasismooth limit Lun et al. (1984); Jenkins and Richman (1985a); Sela and Goldhirsch (1998), the results being subsequently extended to finite degree of inelasticity for monocomponent Brey et al. (1998); Garzó and Dufty (1999) and multicomponent Garzó and Dufty (2002); Garzó et al. (2007a, b); Murray et al. (2012) granular gases.

In spite of the interest and success of the smooth hard-sphere model of granular gases, grains in nature are typically frictional, and hence energy transfer between the translational and rotational degrees of freedom occurs upon particle collisions. The simplest model accounting for particle roughness (and thus including the particle angular velocity as an additional mechanical variable) neglects the effect of sliding collisions and is characterized by a constant coefficient of tangential restitution Jenkins and Richman (1985b). This parameter ranges from (perfectly smooth spheres) to (perfectly rough spheres). A more sophisticated model incorporates the Coulomb friction coefficient as a third collision constant, so that collisions become sliding beyond a certain impact parameter Walton (1993); Foerster et al. (1994); Goldhirsch et al. (2005a, b). On the other hand, this three-parameter model significantly complicates the kinetic description, while the simpler two-parameter model captures the essential features of granular flows when particle rotations are relevant. This explains the wide use of the latter model in the literature Jenkins and Richman (1985b); Lun and Savage (1987); Campbell (1989); Lun (1991); Lun and Bent (1994); Luding (1995); Goldshtein and Shapiro (1995); Lun (1996); Huthmann and Zippelius (1997); Zamankhan et al. (1998); McNamara and Luding (1998); Luding et al. (1998); Herbst et al. (2000); Aspelmeier et al. (2001); Mitarai et al. (2002); Cafiero et al. (2002); Jenkins and Zhang (2002); Polashenski et al. (2002); Goldhirsch et al. (2005a, b); Zippelius (2006); Brilliantov et al. (2007); Gayen and Alam (2008); Kranz et al. (2009); Santos et al. (2010); Nott (2011); Santos et al. (2011); Bodrova and Brilliantov (2012); Mitrano et al. (2013); Vega Reyes et al. (2014).

Needless to say, an important challenge is the derivation of the NSF hydrodynamic equations of a granular gas of inelastic rough hard spheres, with explicit expressions for the transport coefficients as functions of and . Previous attempts have been restricted to nearly elastic collisions () and either nearly smooth particles () Jenkins and Richman (1985b); Goldhirsch et al. (2005a, b) or nearly perfectly rough particles () Jenkins and Richman (1985b); Lun (1991). The goal of this paper is to uncover the whole range of values of the two coefficients of restitution and and derive explicit expressions for the NSF transport coefficients of a dilute granular gas beyond the above limiting situations.

In the case of conventional gases (i.e., when energy is conserved upon collisions), the set of hydrodynamic variables is related to densities of conserved quantities, namely, the particle density (conservation of mass), the flow velocity (conservation of momentum), and the temperature (conservation of energy). If the particles are perfectly elastic and rough (), what is conserved is the sum of the translational and the rotational kinetic energies and thus the granular temperature has translational () and rotational () contributions Pidduck (1922); McCoy et al. (1966); Chapman and Cowling (1970). Moreover, since the angular velocity of the particles is not a collisional invariant, the mean spin is not included either in the set of hydrodynamic variables or in the definition of the rotational contribution () to the temperature Pidduck (1922); McCoy et al. (1966); Chapman and Cowling (1970).

For granular gases, although the total kinetic energy is dissipated by collisions, the granular temperature is typically included as a hydrodynamic field in most studies, and, consequently, a sink term appears in the corresponding balance equation. For nearly smooth spheres (), some authors Jenkins and Richman (1985b); Goldhirsch et al. (2005a, b) have chosen (in addition to and ) the two partial contributions and to the temperature, as well as the mean spin , as hydrodynamic variables. On the other hand, in this paper we will choose as hydrodynamic fields for dissipative gases the same as in conservative systems, i.e., , , and . In this way, the hydrodynamic description encompasses the conservative gases as special limits. The advantage of the choice of the set instead of the set is analogous to the advantage of the set instead of in a binary mixture of inelastic smooth hard spheres, as discussed in Refs. Garzó et al. (2006); Dufty and Brey (2011).

The plan of the paper is as follows. Section II is devoted to the definition of the model of inelastic rough hard spheres and their description in the low-density regime by means of the Boltzmann equation. The exact balance equations for the densities of mass, momentum, and energy are obtained from the Boltzmann equation, and the associated fluxes of momentum and energy, as well as the cooling rate, are identified. The so-called homogeneous cooling state (HCS) is studied in Sec. III. Special attention is paid to the time evolution of the mean spin . While the temperature ratio reaches a well-defined value for long times, the ratio ( being the moment of inertia) decays to zero with a characteristic time typically smaller (except near ) than the relaxation time of the temperature ratio (see Fig. 2). This clearly justifies the exclusion of as a hydrodynamic field. Next, the Chapman–Enskog method is applied in Sec. IV to derive the linear integral equations for the velocity-dependent functions characterizing the distribution function to first-order in the hydrodynamic gradients. In Sec. V, the NSF transport coefficients (shear and bulk viscosities, thermal conductivity, Dufour-like coefficient, and cooling rate transport coefficient) are expressed in terms of integrals involving the solutions of the linear integral equations. Next, a first Sonine polynomial approximation is used to obtain practical results from this formulation, thus providing explicit forms for the transport coefficients as nonlinear functions of both and (see Table 1). Some technical details of the calculations are relegated to the Appendix. The results are discussed in Sec. VI, where known expressions in the limiting cases of inelastic smooth spheres (, ) and perfectly elastic and rough spheres () are recovered. The intricate dependence of the transport coefficients on both coefficients of restitution is illustrated by some representative cases. Finally, the paper closes with some concluding remarks in Sec. VII.

## Ii Granular gas of inelastic rough hard spheres: Boltzmann description

### ii.1 Collision rules

Let and denote the linear and the angular precollisional velocities, respectively, of two rough spherical particles with the same mass , diameter , and moment of inertia , while and correspond to their postcollisional velocities. The pre- and postcollisional velocities are related by

 mv′=mv−Q, Iω′=Iω−σ2ˆσ×Q, (1a) mv′1=mv1+Q, Iω′1=Iω1−σ2ˆσ×Q, (1b)

where denotes the impulse exerted by the unlabeled particle on the labeled one and is the unit collision vector joining the centers of the two colliding spheres and pointing from the center of the unlabeled particle to the center of the labeled one. Furthermore, the relationship between the center-of-mass relative velocities and the relative velocities of the points of the spheres which are in contact during a binary encounter are

 ¯¯¯g=g−σ2ˆσ×(ω+ω1), (2a) ¯¯¯g′=g′−σ2ˆσ×(ω′+ω′1). (2b)

Combining Eqs. (1) and (2), one obtains

 ¯¯¯g′= Extra open brace or missing close brace = ¯¯¯g−2mκ+1κQ+2mκ(ˆσ⋅Q)ˆσ, (3)

where in the second step use has been made of the mathematical property and

 κ=4Imσ2 (4)

is a dimensionless moment of inertia, which may vary from zero to a maximum value of , the former corresponding to a concentration of the mass at the center of the sphere, while the latter corresponds to a concentration of the mass on the surface of the sphere. The value refers to a uniform distribution of the mass in the sphere.

The inelastic collisions of rough spherical particles are characterized by the relationships

 ˆσ⋅¯¯¯g′=−α(ˆσ⋅¯¯¯g),ˆσ×¯¯¯g′=−β(ˆσ×¯¯¯g), (5)

where and are the normal and tangential restitution coefficients, respectively. For an elastic collision of perfectly smooth spheres one has and , while and for an elastic encounter of perfectly rough spherical particles.

Insertion of Eq. (3) into Eq. (5) gives and , where the following abbreviations have been introduced:

 ˜α≡1+α2,˜β≡1+β2κκ+1. (6)

Therefore, the impulse can be expressed as

 Q= m˜α(ˆσ⋅¯g)ˆσ−m˜βˆσ×(ˆσ×¯¯¯g) = m˜α(ˆσ⋅g)ˆσ−m˜βˆσ×(ˆσ×g+σω+ω12). (7)

Equations (1) and (7) express the postcollisional velocities in terms of the precollisional velocities and of the collision vector not (a). From these results it is easy to obtain that the change of the total (translational plus rotational) kinetic energy reads

 ΔK= m2(v′2+v′21−v2−v21)+I2(ω′2+ω′21−ω2−ω21) = −m1−β24κκ+1[ˆσ×(ˆσ×g+σω+ω12)]2 Extra open brace or missing close brace (8)

The right-hand side vanishes for elastic collisions of perfectly smooth spheres and for elastic collisions of perfectly rough spherical particles . In those cases the total energy is conserved in a collision.

Apart from the linear momentum, the angular momentum is conserved, namely,

 m(r×v+r1×v1)+I(ω+ω1)= m(r×v′+r1×v′1) +I(ω′+ω′1), (9)

where and are the position vectors of the two colliding particles.

### ii.2 Boltzmann equation

A direct encounter is characterized by the precollisional velocities , by the postcollisional velocities , and by the collision vector . For a restitution encounter the pre- and postcollisional velocities are denoted by and , respectively, and the collision vector by . It is easy to verify the relationship . The modulus of the Jacobian of the transformation is given by

 ∣∣∣∂(v∗,v∗1;ω∗,ω∗1)∂(v,v1;ω,ω1)∣∣∣=1αβ2 (10)

Thus,

 (ˆσ∗⋅g∗)dv∗dω∗dv∗1dω∗1=1α2β2(ˆσ⋅g)dvdωdv1dω1. (11)

From Eq. (11) we may infer that the Boltzmann equation for granular gases of rough spherical particles without external forces and torques is given by

 ∂f∂t+v⋅∇f=J[f,f], (12a) J[f,f]=σ2∫dv1∫dω1∫+dˆσ(ˆσ⋅g)(f∗1f∗α2β2−f1f), (12b)

where is the one-particle distribution function in the phase space spanned by the positions and the linear and angular velocities of the particles. As usual, in Eq. (12b) the notation , , has been employed. Also, the subscript () in the integration over denotes the constraint .

The so-called transfer equation is obtained from the multiplication of the Boltzmann equation by an arbitrary function and integration of the resulting equation over all values of the velocities and , yielding

 ∂t(n⟨ψ⟩)+∇⋅(n⟨vψ⟩)−n⟨(∂t+v⋅∇)ψ⟩=J[ψ|f,f], (13)

where

 n(r,t)=∫dv∫dωf(r,v,ω,t) (14)

is the local number density,

 Extra open brace or missing close brace (15)

is the local average value of , and

 J[ψ|f,f]≡ ∫dv∫dωψ(r,v,ω,t)J[f,f] = σ22∫dv∫dω∫dv1∫dω1∫+dˆσ(ˆσ⋅g) ×(ψ′1+ψ′−ψ1−ψ)f1f (16)

is the collisional production term of . In the second step of Eq. (II.2) we have used the relationship (11) and the standard symmetry properties of the collision term.

### ii.3 Balance equations

A macroscopic description of a rarefied granular gas of rough spheres can be characterized by the following basic fields: the particle number density [see Eq. (14)], the hydrodynamic flow velocity , and the granular temperature . The two latter quantities are defined as

 u=⟨v⟩,T=12(Tt+Tr), (17)

where

 Tt=m3⟨V2⟩,Tr=I3⟨ω2⟩ (18)

are the (partial) translational and rotational temperatures, respectively, and the averages are defined by Eq. (15). In Eq. (18), is the (translational) peculiar velocity. On the other hand, in the definition of we have chosen not to refer the angular velocities to the mean value

 Ω=⟨ω⟩ (19)

because the latter is not a conserved quantity. Had we defined the granular temperature as with , then would not be a conserved quantity in the case of completely rough and elastic collisions (), even though in that case [see Eq. (8)]

The balance equations for the basic fields are obtained from the transfer equation (13) with the following choices for the arbitrary function :

1. Balance of particle number density (),

 Dtn+n∇⋅u=0. (20)
2. Balance of momentum density (),

 ρDtu+∇⋅P=0. (21)
3. Balance of temperature (),

 DtT+13n(∇⋅q+P:∇u)+Tζ=0. (22)

In the above balance equations, is the mass density, denotes the material time derivative, and the following quantities have been introduced: the pressure tensor

 Pij=ρ⟨ViVj⟩, (23)

the heat flux vector

 q=qt+qr (24)

with

 qt=ρ2⟨V2V⟩,qr=In2⟨ω2V⟩, (25)

and the cooling rate

 ζ=Tt2Tζt+Tr2Tζr (26)

with

 ζt=−m3nTtJ[v2|f,f], (27a) ζr=−I3nTrJ[ω2|f,f]. (27b)

For further use, we note that the hydrostatic pressure is defined as one third of the trace of the pressure tensor, so that

 p=nTt. (28)

Also, the balance equations for the partial temperatures are

 DtTt+23n(∇⋅qt+P:∇u)+Ttζt=0, (29)
 DtTr+23n∇⋅qr+Trζr=0. (30)

Combination of Eqs. (29) and (30) yields Eq. (22).

It is worth noting that the conservation of angular momentum, Eq. (9), does not generate a balance equation with for the quantity because of the difference between the centers of the two colliding spheres.

## Iii Homogeneous cooling state

Before considering the transport properties in inhomogeneous states, it is convenient to characterize the main properties of homogeneous states. In those cases, Eqs. (20) and (21) imply and , while Eqs. (29) and (30) become

 ∂tTt+Ttζt=0,∂tTr+Trζr=0. (31)

The exact forms for the cooling rates and cannot be determined, unless the distribution function is known. Good estimates of those quantities can be obtained by assuming the approximation

 f(v,ω)→(m2πTt)3/2e−mV2/2Ttfr(ω), (32)

where

 fr(ω)=∫dvf(v,ω) (33)

is the marginal distribution of angular velocities. Equation (32) can be justified by maximum-entropy arguments, except that the explicit expression of does not need to be specified. Substitution of (32) into Eqs. (27) and evaluation of the collision integrals given by Eq. (II.2) yields Santos et al. (2010)

 ζt= 512[1−α2+κ1+κ(1−β2)−κ(1+κ)2(1+β)2 ×θ(1+X−1θ)]ν, (34)
 ζr= 5121+β1+κ[(1−β)(1+X)+κ1+κ(1+β) ×(1+X−1θ)]ν, (35)

where

 θ≡TrTt,X≡IΩ23Tr, (36)

and

 ν≡165σ2n√πTt/m (37)

is an effective collision frequency. Note that . The total cooling rate is, according to Eq. (26),

 ζ=51211+θ[1−α2+1−β21+κθ(κθ+1+X)]ν. (38)

Since Eqs. (34) and (35) involve the norm of the mean angular velocity, we need to complement Eq. (31) with the evolution equation for . By using again the approximation (32) one obtains Santos et al. (2010)

 ∂tΩ2+2ζΩΩ2=0,ζΩ=561+β1+κν. (39)

By introducing the time variable , which measures the (nominal) number of collisions per particle from the initial time to time , the solution to Eq. (39) is

 Ω2(s)=Ω2(0)e−2ζ∗Ωs,ζ∗Ω≡ζΩ/ν. (40)

This allows us to solve the set of coupled equations (31) to obtain the time dependence of the temperatures and . Both quantities asymptotically decrease in time due to energy dissipation. On the other hand, the relevant quantity is the temperature ratio . Analogously, rather than the time decay of , and since also tends to decay in time, the relevant quantity is the ratio . The evolution equations for both quantities are

 ∂sθ+(ζ∗r−ζ∗t)θ=0, (41a) ∂sX+(2ζ∗Ω−ζ∗r)X=0, (41b)

where and . Since

 2ζ∗Ω−ζ∗r= 5121+β(1+κ)2[3−X+β(1+X)+2κ(1−X) +κ1+βθ] (42)

is positive definite, it follows that monotonically, no matter the initial condition. On the other hand, the evolution equation for admits a nonzero stationary solution given by the condition . Such a solution is

 θ∞=√1+h2+h (43)

with

 h≡1+κ2κ(1+β)[(1+κ)1−α21+β−(1−κ)(1−β)]. (44)

A standard linear stability analysis of the stationary solution and shows that the two associated eigenvalues are

 ℓ1= θ∞∂(ζ∗r−ζ∗t)∂θ∣∣∣θ=θ∞,X=0 = 512κ(1+β1+κ)21+θ2∞θ∞, (45)
 ℓ2= 2ζ∗Ω−ζ∗r∣∣θ=θ∞,X=0 = 5121+β(1+κ)2(3+β+2κ+κ1+βθ∞). (46)

As expected, both eigenvalues are positive definite, what proves the stability of the stationary solution against homogeneous perturbations. The time evolution of is governed by the eigenvalue only, so that the relaxation time (in units of number of collisions per particle) of is . As for , its relaxation time is if and if . Figure 1 shows the dependence of both relaxation times on roughness for the representative case . Except in a narrow roughness region adjacent to , one has , so that relaxes more rapidly than .

As an illustration of the evolution of and , Fig. 2 shows both quantities for and , , , and , starting from the initial condition , . It is quite apparent that reaches its stationary value (43) in about collisions per particle, while goes to zero in a significantly shorter period.

The special quasismooth limit deserves further comments Santos (2011). In that case, the asymptotic rotational-translational temperature ratio diverges as , the two eigenvalues becoming , which is just the cooling rate of perfectly smooth spheres van Noije and Ernst (1998), and . Therefore, if , relaxes to after a finite number of collisions per particle on the order of . On the other hand, if , evolves in two well-defined stages. The first stage lasts a characteristic time and is very similar to that of the case . This first stage is followed by a much slower relaxation (with collisions per particle) toward the asymptotic value. This singular scenario in the quasismooth limit is illustrated by Fig. 3 for and .

The stationary solution (43) represents the value of the temperature ratio in the HCS Huthmann and Zippelius (1997); Luding et al. (1998); Herbst et al. (2000); Aspelmeier et al. (2001); Zippelius (2006). In such a state the whole time dependence of the distribution function only occurs through the granular temperature , so that the Boltzmann equation (12a) becomes

 −ζT∂Tf=J[f,f], (47)

where, according to Eq. (38), with

 ζ∗=51211+θ∞[1−α2+(1−β2)θ∞+κ1+κ]. (48)

The dependence of the HCS (reduced) cooling rate on both and is displayed in Fig. 4. As can be observed, at a given value of a maximum of occurs around . In the region close to , where , Eq. (48) becomes . Therefore, .

The HCS distribution function has the scaling form

 f(v,ω,t)=n(mIτtτr)3/2[T(t)]−3ϕ(c(t),w(t)), (49)

where

 c(t)=V√2τtT(t)/m,w(t)=ω√2τrT(t)/m, (50)

with the time-independent temperature ratios

 τt≡Tt(t)T(t)=21+θ∞,τr≡Tr(t)T(t)=2θ∞1+θ∞. (51)

Hence, according to Eq. (49), one has the relation

 T∂Tf=−12(∂∂V⋅V+∂∂ω⋅ω)f. (52)

In the particular case of perfectly elastic and smooth particles (, ), the rotational velocities must be ignored (i.e., , ) and the solution of Eq. (47) is just the equilibrium distribution of translational velocities. In the other conservative case of perfectly elastic and rough spheres (), the solution of Eq. (47) is the equilibrium distribution with a common temperature (i.e., ). In the general dissipative case, however, the solution of Eq. (47) is not exactly known, although good approximations are given in the form of a two-temperature Maxwellian multiplied by truncated Sonine polynomial expansions Huthmann and Zippelius (1997); McNamara and Luding (1998); Luding et al. (1998); Herbst et al. (2000); Aspelmeier et al. (2001); Cafiero et al. (2002); Zippelius (2006); Brilliantov et al. (2007); Kranz et al. (2009); Bodrova and Brilliantov (2012); Santos et al. (2011); Vega Reyes et al. (2014).

## Iv Chapman–Enskog method

### iv.1 Outline of the method

To determine the distribution function from the Chapman–Enskog method Chapman and Cowling (1970); Ferziger and Kaper (1972), we write the Boltzmann equation (12a) as

 Dtf+ϵV⋅∇f=J[f,f], (53)

where is a uniformity parameter (set equal to unity at the end of the calculations) measuring the strength of the spatial gradients. According to the method, the distribution function and the material time derivative are expanded in terms of the parameter as follows,

 f=f(0)+ϵf(1)+ϵ2f(2)+⋯, (54)
 Dt=D(0)t+ϵD(1)t+ϵ2D(2)t+⋯. (55)

Substitution of Eq. (54) into the definitions of the fluxes and the cooling rate gives

 Pij=p(0)δij+ϵP(1)ij+ϵ2P(2)ij+⋯, (56)
 q=ϵq(1)+ϵ2q(2)+⋯, (57)
 ζ=