Energy production rates in fluid mixtures of inelastic rough hard spheres

# Energy production rates in fluid mixtures of inelastic rough hard spheres

Andrés Santos,*)*)*)E-mail: andres@unex.es Gilberto M. Kremer†)†)†)E-mail: kremer@fisica.ufpr.br and Vicente Garzó‡)‡)‡)E-mail: vicenteg@unex.es

## 1 Introduction

A granular fluid is usually modeled as a system of identical, inelastic smooth hard spheres with a constant coefficient of normal restitution . Despite its simplicity, this model has been useful to capture the basic properties of granular flows.[1] On the other hand, the model can be made closer to reality by introducing more ingredients, such as coefficients of normal restitution depending on the impact velocity,[2] presence of an interstitial fluid,[3] non-spherical shapes,[4] polydispersity,[5] or roughness.[6] Of course, the few citations in the preceding sentence are just representative of many works reporting features not accounted for by the simple monodisperse smooth-sphere model.

In this paper we will focus on the two latter ingredients, namely polydispersity and roughness. These properties are especially relevant, not only because beads and grains are unavoidably polydisperse and rough, but also because any of these ingredients unveils an inherent breakdown of energy equipartition in granular fluids, even in homogeneous and isotropic states. In the case of multi-component granular fluids, most of the studies have considered the inelastic smooth-sphere description. Some of the problems addressed include non-equipartition in homogeneous states,[7, 8, 9, 10] Navier–Stokes transport coefficients,[11, 12, 13, 14, 15] and segregation phenomena.[16, 17, 18, 19, 20]

Concerning the case of inelastic rough spheres, most of the works we are aware of are restricted to monodisperse systems.[21, 22, 23, 24, 25, 30, 26, 27, 31, 32, 33, 34, 35, 39, 36, 28, 29, 38, 40, 37, 41] Analogously to what happens with the coefficient of normal restitution , the simplest model accounting for friction during collisions assumes a constant coefficient of tangential restitution . While is a positive quantity smaller than or equal to 1 (the value corresponding to elastic spheres), the parameter lies in the range between (perfectly smooth spheres) to (perfectly rough spheres). Except for and , the total kinetic energy is not conserved in a collision. Some of the early attempts to develop a kinetic theory for rough spheres were carried out by Jenkins and Richman[21] and Lun and Savage,[22] who applied their approaches to the simple shear flow problem. The influence of roughness in shear flows has also been studied by several authors,[23, 24, 25, 26, 27, 28, 29] usually assuming that the spheres are nearly smooth and nearly elastic. In an extensive paper,[30] Goldshtein and Shapiro obtained the collisional energy production rates associated with the translational and rotational degrees of freedom by using Maxwellian forms for the distribution functions. The result was applied to the evaluation of the ratio between the translational () and rotational () temperatures in the homogeneous cooling state (HCS). The time evolution of the ratio towards its HCS asymptotic value has been widely analyzed, both theoretically and by means of molecular dynamics, by Luding, Zippelius, and co-workers.[31, 32, 33, 34, 35, 36, 37] Other studies involving roughness include vibration with rough walls,[38] a micropolar fluid model for granular flows on a slope,[39] derivation of hydrodynamic constitutive equations from the Boltzmann equation for nearly smooth, nearly elastic granular gases,[40] and correlations between the rotation axis and the translational direction.[41]

The studies about multi-component rough-sphere systems are much scarcer. To the best of our knowledge, only the case of a fixed particle immersed in a bath of thermalized point particles has been addressed.[42, 43, 44] On the other hand, the general case of a polydisperse system made of mobile particles of different coefficients of normal and tangential restitution ( and ) has not been studied yet. In this paper, we address one of the basic aspects of the problem, namely those related to the partition of the total energy. In order to characterize the effect of collisions on energy partition, we focus on the partial productions rates and measuring the rate of change of the translational and rotational kinetic energies, respectively, of particles of component due to collisions with particles of component . A combination of and gives the total cooling rate of the mixture.

Starting from the collisional rules worked out in Section 2, the collisional rates of change of momentum (linear and angular) and energy (translational and rotational) are expressed as linear combinations of two-body average values in Section 3. These averages are evaluated in terms of the partial temperatures , , , and by assuming a maximum-entropy two-body distribution in Section 4. The expressions for , , and obtained in Section 5 extend previous results derived for monodisperse rough spheres[30, 33] and for polydisperse smooth spheres.[7] An application to the HCS of a binary mixture is carried out in Section 6 with some illustrative examples. Finally, some concluding remarks are presented in Section 7.

## 2 Collisional rules

Let us consider the collision between two hard spheres of masses and , diameters and , and moments of inertia and . The latter two quantities can be equivalently characterized by the dimensionless parameters

 κi≡4Iimiσ2i,κj≡4Ijmjσ2j. (2.1)

The value of depends on the mass distribution within the sphere and runs from the extreme values (mass concentrated on the center) to (mass concentrated on the surface). If the mass is uniformly distributed, then . Let us denote by the pre-collisional relative velocity of the center of mass of both spheres and by and the respective pre-collisional angular velocities. This is sketched in Fig. 1, where is the unit vector pointing from the center of to the center of . The velocities of the points of the spheres which are in contact during the collision are

 wi=vi−σi2ˆσ×ωi,wj=vj+σj2ˆσ×ωj, (2.2)

the corresponding relative velocity being

 wij=vij−ˆσ×Sij,Sij≡σi2ωi+σj2ωj. (2.3)

Post-collisional velocities will be denoted by primes. Conservation of linear and angular momenta yields[37]

 miv′i+mjv′j=mivi+mjvj, (2.4)
 (2.5)

Equations (2) and (2) imply that

 v′i=vi−1miQij,v′j=vj+1mjQij, (2.6)
 ω′i=ωi−σi2Iiˆσ×Qij,ω′j=ωj−σj2Ijˆσ×Qij, (2.7)

where is the impulse exerted by particle on particle . Therefore,

 v′ij=vij−1mijQij,w′ij=wij−1mijQij+1mijκijˆσ×(ˆσ×Qij), (2.8)

where

 mij≡mimjmi+mj,κij≡κiκjmi+mjκimi+κjmj (2.9)

are the reduced mass and a sort of reduced inertia-moment parameter, respectively.

To close the collisional rules, we need to express in terms of the pre-collisional velocities and the unit vector . To that end, we relate the normal (i.e., parallel to ) and tangential (i.e., orthogonal to ) components of the relative velocities and by

 ˆσ⋅w′ij=−αijˆσ⋅wij,ˆσ×w′ij=−βijˆσ×wij. (2.10)

Here, and are the coefficients of normal and tangential restitution, respectively. The former coefficient ranges from (perfectly inelastic particles) to (perfectly elastic particles), while the latter runs from (perfectly smooth particles) to (perfectly rough particles). Inserting the second equality of Eq. (2) into Eq. (2) one simply gets and , where

 ˜αij≡mij(1+αij),˜βij≡mijκij1+κij(1+βij). (2.11)

Therefore,

 Qij=˜αij(vij⋅ˆσ)ˆσ+˜βij[vij−ˆσ×Sij−(vij⋅ˆσ)ˆσ], (2.12)

where use has been made of the mathematical property . Note that in the special case of perfectly smooth spheres () one has , so that . In that case, according to Eq. (2), the angular velocities of the two colliding spheres are unaltered by the collision.

The total kinetic energy before collision is

 Eij=mi2v2i+mj2v2j+Ii2ω2i+Ij2ω2j. (2.13)

It can be checked (see the Appendix) that

 E′ij−Eij = −mij2κij1+κij(1−β2ij)[(ˆσ×vij)2+(ˆσ×Sij)2+2(ˆσ×vij)⋅Sij] (2.14) −mij2(1−α2ij)(ˆσ⋅vij)2.

The right-hand side is a negative definite quantity. Thus, we observe that energy is conserved only if the particles are elastic () and either perfectly smooth () or perfectly rough (). Otherwise, and kinetic energy is dissipated upon collisions.

Equations (2), (2), and (2) give the direct collisional rules. The restituting collisional rules are

 v′′i=vi−1miQ−ij,v′′j=vj+1mjQ−ij, (2.15)
 ω′′i=ωi−σi2Iiˆσ×Q−ij,ω′′j=ωj−σj2Ijˆσ×Q−ij, (2.16)

where

 Q−ij=˜αijαij(vij⋅ˆσ)ˆσ+˜βijβij[vij−ˆσ×Sij−(vij⋅ˆσ)ˆσ]. (2.17)

Here the double primes denote pre-collisional quantities giving rise to unprimed quantities as post-collisional values. The modulus of the Jacobian of the transformation between pre- and post-collisional velocities is

 ∣∣ ∣∣∂(v′i,ω′i,v′j,ω′j)∂(vi,ωi,vj,ωj)∣∣ ∣∣=∣∣ ∣∣∂(vi,ωi,vj,ωj)∂(v′′i,ω′′i,v′′j,ω′′j)∣∣ ∣∣=αijβ2ij. (2.18)

## 3 Collisional rates of change

Let be the two-body distribution function with the normalization condition

 ∫dri∫dvi∫dωi∫drj∫dvj∫dωjf(2)ij(ri,vi,ωi;rj,vj,ωj;t)=NiNj, (3.1)

being the number of spheres of component . The one-body distribution function is

 (3.2)

The marginal distribution functions associated with the translational and rotational degrees of freedom are

 ftri(ri,vi;t)=∫dωifi(ri,vi,ωi;t),froti(ri,ωi;t)=∫dvifi(ri,vi,ωi;t). (3.3)

Given a one-body function , we define its average as

 ⟨ψ(vi,ωi)⟩≡1ni∫dvi∫dωiψ(vi,ωi)fi(vi,ωi),ni=∫dvi∫dωifi(vi,ωi), (3.4)

where is the number density of component and, for the sake of brevity, we have omitted the spatial and temporal arguments. In particular, one can define partial temperatures associated with the translational and rotational degrees of freedom as

 Ttri=mi3⟨(vi−u)2⟩,Troti=Ii3⟨ω2i⟩, (3.5)

where

 u=∑imini⟨vi⟩∑imini (3.6)

is the flow velocity. Note that in the definition of we have not referred the angular velocities to any average value because of the lack of invariance under the addition of a common vector to every angular velocity. The global temperature is

 T=∑ini2n(Ttri+Troti), (3.7)

where is the total number density.

By starting from the Liouville equation and following standard steps, one can derive the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy.[45]  The first equation of the hierarchy reads

 ∂tfi(ri,vi,ωi;t)+vi⋅∇fi(ri,vi,ωi;t)=∑jJij[ri,vi,ωi;t|f(2)ij], (3.8)

where

 Jij[ri,vi,ωi;t|f(2)ij] = σ2ij∫dvj∫dωj∫dˆσ\mathchar28930\relax(vij⋅ˆσ)(vij⋅ˆσ) (3.9) ×[1α2ijβ2ijf(2)ij(ri,v′′i,ω′′i;ri−σijˆσ,v′′j,ω′′j;t) −f(2)ij(ri,vi,ωi;ri+σijˆσ,vj,ωj;t)]

is the collision operator. Here, and use has been made of Eqs. (2) and (2.18).

Multiplying both sides of Eq. (3.8) by and integrating over and one gets

 ∂tni⟨ψ(vi,ωi)⟩+∇⋅ni⟨viψ(vi,ωi)⟩=∑jJij[ψ(vi,ωi)|f(2)ij], (3.10)

where

 Jij[ψ(vi,ωi)|f(2)ij] ≡ ∫dvi∫dωiψ(vi,ωi)Jij[vi,ωi|f(2)ij] = σ2ij∫dvi∫dωi∫dvj∫dωj∫dˆσ\mathchar28930\relax(vij⋅ˆσ)(vij⋅ˆσ) ×f(2)ij(ri,vi,ωi;ri+σijˆσ,vj,ωj)[ψ(v′i,ω′i)−ψ(vi,ωi)].

Thus, is the rate of change of the quantity due to collisions with particles of component . This rate of change is a functional of the two-body distribution function , as indicated by the notation, and it is in general a rather intricate quantity. The physically important cases are . The corresponding rates of change are obtained by inserting Eqs. (.5)–(.8) into Eq. (LABEL:3). Note that so far all the results are formally exact.

To proceed, let us make the approximation

 Jij[ψ(vi,ωi)|f(2)ij]≈Jij[ψ(vi,ωi)|¯f(2)ij], (3.12)

where

 ¯f(2)ij(ri,vi,ωi;vj,ωj)≡∫dˆσ\mathchar28930\relax(vij⋅ˆσ)(vij⋅ˆσ)f(2)ij(ri,vi,ωi;ri+σijˆσ,vj,ωj)∫dˆσ\mathchar28930\relax(vij⋅ˆσ)(vij⋅ˆσ) (3.13)

is the orientational average of the pre-collisional distribution . Thus, Eq. (3.12) replaces a detailed functional of by a simpler one where the solid angle integral

 ∫dˆσ\mathchar28930\relax(vij⋅ˆσ)(vij⋅ˆσ)[ψ(v′i,ω′i)−ψ(vi,ωi)] (3.14)

can be evaluated independently of . It is important to borne in mind that the approximation (3.12) is much weaker than the approximation . Notwithstanding this, the equality holds if (a) the system is homogeneous and isotropic (regardless of the reduced densities and ), in which case only depends on , or (b) the system is in the Boltzmann limit (, ), in which case one can formally take in the contact value of . Therefore, the approximation (3.12) is justified if the density of the granular gas and/or its inhomogeneities are small enough so the value of at contact is hardly dependent on the relative orientation of the two colliding spheres.

In the remainder of this Section we particularize to and express the rates of change in terms of two-body averages of the form

 ⟨A(vi,ωi;vj,ωj)⟩ ≡ 1ninj∫dvi∫dωi∫dvj∫dωjA(vi,ωi;vj,ωj) (3.15) ×¯f(2)ij(vi,ωi;vj,ωj).

The results are (see the Appendix)

 n−1iJij[mivi]=−njσ2ijπ(˜αij+˜βij2⟨vijvij⟩−2˜βij3⟨vij×Sij⟩), (3.16)
 n−1iJij[Iiωi]=−njσ2ijσiπ8˜βij[3⟨vijSij⟩−⟨v−1ij(vij⋅Sij)vij⟩], (3.17)
 n−1iJij[miv2i] = −njσ2ijπ[(˜αij+˜βij)⟨vijvi⋅vij⟩+4˜βij3⟨Sij⋅(vi×vj)⟩ (3.18) −˜α2ij+˜β2ij2mi⟨v3ij⟩−3˜β2ij4mi⟨vijS2ij⟩+˜β2ij4mi⟨v−1ij(vij⋅Sij)2⟩⎤⎥⎦,
 n−1iJij[Iiω2i] = −njσ2ijπ4˜βij{3σi⟨vijωi⋅Sij⟩−σi⟨v−1ij(vij⋅Sij)(vij⋅ωi)⟩ (3.19) −˜βijmiκi[2⟨v3ij⟩+3⟨vijS2ij⟩−⟨v−1ij(vij⋅Sij)2⟩]}.

## 4 Estimates of the average values

Equations (3.16)–(3.19) express the collisional rates of change of the main quantities as linear combinations of two-body averages of the form (3.15). They are local functions of space and time and functionals of the orientation-averaged pre-collisional distribution . While, thanks to the approximation (3.12), Eqs. (3.16)–(3.19) are much more explicit than the exact results obtained from Eq. (LABEL:3), they still require the full knowledge of .

Suppose, for simplicity, that and define the average angular velocities

 ⟨ωi⟩=\mathchar28938\relaxi,⟨ωj⟩=\mathchar28938\relaxj. (4.1)

Now, let us imagine that, instead of the full knowledge of , we only know the local values of the two densities ( and ), the two average angular velocities ( and ), and the four partial temperatures (, , , and ). The question we want to address in this section is, can we get reasonable estimates of the two-body averages appearing in Eqs. (3.16)–(3.19) by expressing them in terms of , , , , , , , and ? In the absence of further information, the least biased estimates are obtained from the replacement

 ¯f(2)ij(vi,ωi;vj,ωj) → χij⎛⎜⎝mimj4π2TtriTtrj⎞⎟⎠3/2exp⎡⎢⎣−mi(vi−u)22Ttri−mj(vj−u)22Ttrj⎤⎥⎦ (4.2) ×froti(ωi)frotj(ωj),

where is the contact value of the pair correlation function. Equation (4.2) can be justified by maximum-entropy arguments, except that here we do not need to assume a Maxwellian form for the rotational distributions, given that the angular velocities only appear linearly or quadratically in Eqs. (3.16)–(3.19). Like in the approximation (3.12), it is important to stress that we are not making the strong claim that is well approximated by the right-hand side of Eq. (4.2) [see Ref. 41] but only the wekaer one that the two-body averages can be estimated by performing such a replacement. Table I displays those estimates.

Although has been assumed in the results shown in Table I, the generalization to can be carried out following the same steps as done in Ref. 46 for smooth spheres.

## 5 Energy production rates and cooling rate

The most characteristic feature of a granular gas is the energy dissipation taking place after each collision. In the model of inelastic rough hard spheres this is clearly apparent from Eq. (2.14). On the other hand, any of the four partial kinetic energies contributing to in Eq. (2) can either increase or decrease after a given collision. To characterize this effect at a statistical level, it is convenient to introduce the rates of change of the partial temperatures and due to collisions of particles of component with particles of component . More explicitly, we define the (partial) energy production rates and as

 ξtrij≡−13niTtriJij[mi(vi−u)2],ξrotij≡−13niTrotiJij[Iiω2i]. (5.1)

When collisions of particles of component with all the components are considered, we get the (total) energy production rates

 ξtri≡−1Ttri⎛⎝∂Ttri∂t⎞⎠coll=∑jξtrij,ξroti≡−1Troti⎛⎝∂Troti∂t⎞⎠coll=∑jξrotij. (5.2)

Finally, the net cooling rate is

 ζ≡−1T(∂T∂t)coll=∑ini2nT(Ttriξtri+Trotiξroti). (5.3)

In contrast to the energy production rates defined in Eqs. (5.1) and (5.2), the cooling rate is positive definite, i.e., collisions produce a decrease of the total temperature unless and for all pairs .

From Eqs. (3.16)–(3.19) and the expressions of Table I one gets the following estimates:

 n−1iJij[mivi]=0,n−1iJij[Iiωi]=−νij˜βij4σi(σi\mathchar28938\relaxi+σj\mathchar28938\relaxj), (5.4)
 ξtrij = νijmiTtri⎡⎢⎣2(˜αij+˜βij)Ttri−(˜α2ij+˜β2ij)⎛⎜⎝Ttrimi+Ttrjmj⎞⎟⎠ (5.5) −˜β2ij⎛⎜⎝Trotimiκi+Trotjmjκj+16σiσj\mathchar28938\relaxi⋅\mathchar28938\relaxj⎞⎟⎠⎤⎥⎦,
 ξrotij = νi