Impact of roughness on the instability of a free-cooling granular gas
aa r X i v : . [ c ond - m a t . s o f t ] M a y Impact of roughness on the instability of a free-cooling granular gas
Vicente Garz´o ∗ and Andr´es Santos † Departamento de F´ısica and Instituto de Computaci´on Cient´ıfica Avanzada (ICCAEx),Universidad de Extremadura, E-06006 Badajoz, Spain
Gilberto M. Kremer ‡ Departamento de F´ısica, Universidade Federal do Paran´a, 81531-980 Curitiba, Brazil (Dated: November 13, 2018)A linear stability analysis of the hydrodynamic equations with respect to the homogeneous coolingstate is carried out to identify the conditions for stability of a granular gas of rough hard spheres. Thedescription is based on the results for the transport coefficients derived from the Boltzmann equationfor inelastic rough hard spheres [Phys. Rev. E , 022205 (2014)], which take into account thecomplete nonlinear dependence of the transport coefficients and the cooling rate on the coefficientsof normal and tangential restitution. As expected, linear stability analysis shows that a doublydegenerate transversal (shear) mode and a longitudinal (“heat”) mode are unstable with respect tolong enough wavelength excitations. The instability is driven by the shear mode above a certaininelasticity threshold; at larger inelasticity, however, the instability is driven by the heat modefor an inelasticity-dependent range of medium roughness. Comparison with the case of a granulargas of inelastic smooth spheres confirms previous simulation results about the dual role playedby surface friction: while small and large levels of roughness make the system less unstable thanthe frictionless system, the opposite happens at medium roughness. On the other hand, such anintermediate window of roughness values shrinks as inelasticity increases and eventually disappearsat a certain value, beyond which the rough-sphere gas is always less unstable than the smooth-sphere gas. A comparison with some preliminary simulation results shows a very good agreementfor conditions of practical interest. I. INTRODUCTION
One of the most characteristic features of granular flu-ids, as compared with ordinary fluids, is the spontaneousformation of velocity vortices and density clusters infreely evolving flows [homogenous cooling state (HCS)].This kind of flow instability originates by the dissipativenature of particle-particle collisions [1]. In the case ofsmooth and frictionless inelastic hard spheres, those in-stabilities were first detected independently, by means ofcomputer simulations, by Goldhirsch and Zanetti [2] andMcNamara [3]. An important property of the cluster-ing instability is that it is restricted to long-wavelengthexcitations, and, hence, it can be suppressed for systemsthat are small enough. In addition, simulations also showthat vortices generally preempt clusters, so that the onsetof instability is generally associated with the transversalshear modes.Instabilities of freely cooling flows of smooth spherescan be well described by means of a linear stability anal-ysis of the Navier–Stokes (NS) hydrodynamic equations.This analysis gives a critical wave number k c and an asso-ciated critical wavelength L c = 2 π/k c , such that the sys-tem becomes unstable if its linear size is larger than L c .The evaluation of L c requires knowledge of the depen-dence of the NS transport coefficients on the coefficient ∗ † ‡ kremer@fisica.ufpr.br of normal restitution α , and therefore the determinationof L c is perhaps one of the most interesting applicationsof the NS equations. In addition, comparison betweenkinetic theory and computer simulations for the criticalsize can be considered a clean way of assessing the degreeof accuracy of the former, since theoretical results are ob-tained in most of the cases under certain approximations(e.g., by considering the leading terms in a Sonine poly-nomial expansion).The accuracy of the kinetic-theory prediction of L c for a low-density monodisperse granular gas of inelasticsmooth hard spheres [4, 5] has been verified by means ofthe direct simulation Monte Carlo DSMC method [6] and,more recently, by molecular dynamics (MD) simulationsfor a granular fluid at moderate density [7, 8]. Similarcomparisons have also been made for polydisperse sys-tems [9–12] in the low-density regime [13] and moderatedensities [14]. In general, the theoretical predictions forthe critical size compare quite well with computer simu-lations, even for strong dissipation.On the other hand, all the above works ignore the influ-ence of collisional friction on instability. Needless to say,grains in nature are typically frictional, and, hence, en-ergy transfer between the translational and rotational de-grees of freedom occurs upon particle collisions. The sim-plest way of accounting for friction is, perhaps, througha constant coefficient of tangential restitution β , rang-ing from − α .
1) and ei-ther nearly smooth particles ( β & −
1) [19–21] or nearlyperfectly rough particles ( β .
1) [19, 22]. An exten-sion of the previous works to arbitrary values of α and β has been recently carried out by Kremer et al. [23]for a dilute granular gas. Explicit expressions for the NStransport coefficients and the cooling rate were obtainedas nonlinear functions of both coefficients of restitutionand the moment of inertia. The knowledge of the trans-port coefficients opens up the possibility of performinga linear stability analysis of the hydrodynamic equationswith respect to the HCS state to identify the conditionsfor stability as functions of the wave vector and the co-efficients of restitution.A previous interesting work on flow instabilities for adense granular gas of inelastic and rough hard sphereswas carried out by Mitrano et al. [24]. In that pa-per, the authors compare their MD simulations againsttheoretical predictions obtained from a simple stabilityanalysis where friction is only accounted for through itsimpact on the cooling rate, since the NS transport coeffi-cients are otherwise assumed to be formally the same asthose obtained for frictionless particles [16]. In spite ofthe simplicity of this theoretical approach, the obtainedpredictions compare well with MD simulations. In par-ticular, they observe that, paradoxically, high levels ofroughness can actually attenuate instabilities relative tothe frictionless case.Even though, as said above, the predictions made byMitrano et al. [24] agree reasonably well with simulationresults, they lack a sounder basis as they are essentiallybased on the hydrodynamic description of smooth parti-cles and only the cooling rate incorporates the completenonlinear dependence on β . As a matter of fact, thederivation of the transport coefficients for rough spheres[23] was not known before Ref. [24] was published. There-fore, it is worthwhile assessing to what extent the pre-vious results are indicative to what happens when thecomplete nonlinear dependence of both α and β on thetransport coefficients is included in the stability analysis.The main aim of this paper is to fill this gap by revis-iting the problem of the HCS stability of a granular gasdescribed by the Boltzmann kinetic equation, this timemaking use of the knowledge of the NS transport coef-ficients of a granular gas of inelastic rough hard spheres[23].As expected, linear stability analysis shows twotransversal (shear) modes and a longitudinal (“heat”)mode to be unstable with respect to long enough wave-length excitations. The corresponding critical values forthe shear and heat modes are explicitly determined asfunctions of the (reduced) moment of inertia and the co-efficients of restitution α and β . The results show thatthe instability is mainly driven by the transversal shear mode, except for values of the coefficient of normal resti-tution α smaller than a certain threshold value α th . Be-low that value, the instability is dominated by the heatmode, provided that the coefficient of tangential restitu-tion β lies inside an α -dependent window of values around β ≈
0. Moreover, below an even smaller value α D , anyvalue of β is enough to attenuate the instability effectwith respect to the frictionless case.The plan of the paper is as follows. The NS hydrody-namic equations for a dilute granular gas of inelastic andrough hard spheres are shown in Sec. II, the explicit ex-pressions for the transport coefficients being given in theAppendix. Next, the linear stability analysis of the NSequations is worked out in Sec. III, where it is found thatthe two transversal shear modes are decoupled from thethree longitudinal ones. Section IV deals with a detaileddiscussion of the results, with a special emphasis on theimpact of roughness on the instability of the HCS. Thepaper ends in Sec. V with a summary and concludingremarks. II. NAVIER–STOKES HYDRODYNAMICEQUATIONS
We consider a dilute granular gas composed of inelas-tic and rough hard spheres of mass m , diameter σ , andmoment of inertia I = κmσ / κ characterizes the mass distribution ina sphere). Since, in general, particle rotations and sur-face friction are relevant in the description of granularflows, collisions are characterized by constant coefficientsof normal restitution ( α ) and tangential restitution ( β ).The coefficient α measures the postcollisional shrinkingin the magnitude of the normal component of the rela-tive velocity of the two points at contact; the coefficient β does the same but for the tangential component [25, 26].As said in Sec. I, the coefficient α ranges from 0 (per-fectly inelastic particles) to 1 (perfectly elastic particles)while the coefficient β ranges from − n ( r , t ), the local tempera-ture T ( r , t ), and the flow velocity u ( r , t ): D t n + n ∇ · u = 0 , (2.1a) D t T + 13 n ( ∇ · q + P : ∇ u ) + ζT = 0 , (2.1b) mnD t u + ∇ · P = 0 . (2.1c)In Eqs. (2.1), D t = ∂ t + u · ∇ is the material derivative, P is the pressure tensor, q is the heat flux, and ζ is thecooling rate due to the collisional energy dissipation.The balance equations (2.1) do not constitute a closedset of equations for the hydrodynamic fields { n, T, u } ,unless the fluxes and the cooling rate are further speci-fied in terms of the hydrodynamic fields and their gra-dients. As mentioned in Sec. I, the detailed form of theconstitutive equations and the transport coefficients ap-pearing in them have recently been obtained by applyingthe Chapman–Enskog method to the Boltzmann equa-tion. To first order in the gradients, the correspondingconstitutive equations are [23] P ij = nT τ t δ ij − η (cid:18) ∇ j u i + ∇ i u j − δ ij ∇ · u (cid:19) − η b δ ij ∇· u , (2.2a) q = − λ ∇ T − µ ∇ n, (2.2b) ζ = ζ (0) − ξ ∇ · u . (2.2c)Here, τ t ≡ T (0) t /T is the ratio between the (partial) trans-lational temperature T (0) t and the granular temperature T in the reference HCS, η is the shear viscosity, η b isthe bulk viscosity, λ is the thermal conductivity, µ is anew transport coefficient (heat diffusivity coefficient orDufour-like coefficient), not present in the conservativecase, and ξ characterizes the NS correction to the HCScooling rate ζ (0) . The last three coefficients have bothtranslational ( λ t , µ t , ξ t ) and rotational ( λ r , µ r , ξ r ) con-tributions. Moreover, the (partial) rotational temper-ature T (0) r in the HCS is given by T (0) r = τ r T , where τ r = 2 − τ t .Insertion of the NS constitutive equations (2.2) intothe balance equations (2.1) yields the corresponding NShydrodynamic equations for n , T , and u : D t n = − n ∇ · u , (2.3a) (cid:16) D t + ζ (0) (cid:17) T = (cid:16) ξ − τ t (cid:17) T ∇ · u + 13 n ∇ · ( λ ∇ T + µ ∇ n )+ 13 n h η ( ∇ i u j + ∇ j u i ) − (cid:18) η − η b (cid:19) × δ ij ∇ · u i ∇ i u j , (2.3b) mnD t u i = − τ t ∇ i ( nT ) + ∇ j h η ( ∇ i u j + ∇ j u i ) − (cid:18) η − η b (cid:19) δ ij ∇ · u i . (2.3c)As noted in previous stability studies [5, 10], in princi-ple one should consider terms up to second order in thegradients in Eq. (2.2c) for the cooling rate, since this isthe order of the terms of Eq. (2.3b) coming from the pres-sure tensor and the heat flux. On the other hand, it hasbeen shown for dilute granular gases of smooth spheresthat the second order contributions to the cooling rate arenegligible as compared with their corresponding zeroth- and first-order contributions [4]. It is assumed here thatthe same happens for inelastic rough hard spheres.To close the hydrodynamic problem posed by Eqs.(2.3), one needs to express the coefficients, τ t , ζ (0) , η , η b , λ , µ , and ξ as functions of the hydrodynamic fields( n and T ) and of the mechanical parameters ( α , β , and κ ). The former dependence is dictated by dimensionalanalysis, while the dependence on α , β , and κ requiresresorting to approximations.In the two-temperature Maxwellian approximation,the HCS parameters τ t , τ r , and ζ (0) are given by[25, 26, 29, 30] τ t = 21 + θ , τ r = 2 θ θ , θ = p h + h, (2.4a) ζ (0) = ζ ∗ ν, ν = 165 nσ χ ( φ ) r πT τ t m , (2.4b)where θ = T (0) r /T (0) t is the rotational-to-translationaltemperature ratio and h = (1 + κ ) κ (1 + β ) (cid:20) − α − (1 − β ) 1 − κ κ (cid:21) , (2.5a) ζ ∗ = 512 11 + θ (cid:20) − α + (1 − β ) κ + θ κ (cid:21) . (2.5b)In the second equality of Eq. (2.4b), χ ( φ ) is the con-tact value of the pair correlation function (Enskog fac-tor), which is a function of the solid volume fraction φ ≡ π nσ . A simple and accurate prescription is χ ( φ ) =(1 − φ/ / (1 − φ ) [31]. In principle, the results of thispaper apply to the Boltzmann limit φ →
0, in which case χ ( φ ) →
1. Nevertheless, keeping the Enskog factor χ ( φ )in the collision frequency ν allows us to account for basicexcluded volume effects which are present if φ is smallbut not strictly zero.From Eqs. (2.4a) and (2.5a), we note that θ divergesin the quasismooth limit β → − θ = (1 − α )(1 + κ ) κ (1 + β ) − . (2.6)Consequently, in that limit the reduced cooling rate van-ishes as ζ ∗ = 56(1 + κ ) (1 + β ) . (2.7)This contrasts with the purely smooth case [4], where ζ ∗ sm = 512 (1 − α ) , (2.8)which can be obtained from Eq. (2.5b) by setting β = − θ = 0. The difference between Eqs.(2.7) and (2.8) illustrates the strong singular character ofthe limit β → − η = η η ∗ , η b = η η ∗ b , λ = λ λ ∗ , µ = T λ n µ ∗ , (2.9)where η = nT τ t ν , λ = 154 η m (2.10)are the shear viscosity and the thermal conductivity co-efficients, respectively, of a gas of elastic and smoothspheres at the (translational) temperature T τ t . The di-mensionless quantities η ∗ , η ∗ b , λ ∗ , and µ ∗ are nonlinearfunctions of κ , α , and β given by η ∗ = 1 ν ∗ η − ζ ∗ , η ∗ b = τ r γ E , (2.11a) λ ∗ = 23 τ t γ A t + 25 τ r γ A r , µ ∗ = 23 τ t γ B t + 25 τ r γ B r , (2.11b)where the explicit forms of the quantities ν ∗ η , γ E , γ A t , γ A r , γ B t , and γ B r are displayed in the Appendix [seeEqs. (A.1) and (A.3)].Finally, the NS cooling-rate coefficient ξ is ξ = γ E Ξ , (2.12)where Ξ is given by Eq. (A.5).In the purely smooth case ( β = − τ t = 2, τ r = 0, ξ = 0, η ∗ b = 0, ζ ∗ = ζ ∗ sm , η ∗ = η ∗ sm , λ ∗ = 2 λ ∗ sm ,and µ ∗ = 2 µ ∗ sm , with η ∗ sm = 24(1 + α )(13 − α ) , (2.13a) λ ∗ sm = 32(1 + α )(9 + 7 α ) , (2.13b) µ ∗ sm = 640(1 − α )(1 + α )(9 + 7 α )(19 − α ) . (2.13c) III. LINEAR STABILITY ANALYSIS OF THENS EQUATIONS
It is well known that the NS equations admit a simplesolution corresponding to the so-called HCS. It describesa uniform state with vanishing flow field and a tempera-ture T decreasing monotonically in time, ∇ n H = ∇ T H = 0 , ∂ t ln T H = − ζ (0) H , u H = , (3.1)where henceforth the subscript H denotes quantities eval-uated in the HCS. On the other hand, the HCS is known to become un-stable under sufficiently long wavelength excitations [2–14, 24]. Our aim here is to study the stability conditionsof the HCS by means of a linear stability analysis of Eqs.(2.3) for small initial excitations with respect to the ho-mogeneous state. Thus, we assume that the deviations δy a ( r , t ) = y a ( r , t ) − y Ha ( t ) of the hydrodynamic fields { y a ; a = 1 , . . . , } ≡ { n, T, u } from their values in theHCS are small, and neglect terms of second and higherorder. The resulting equations for δn , δT , and δu i are ∂ t δnn H = − v H ∇ · δ u v H , (3.2a) ∂ t δTT H = − ζ (0) H (cid:18) δnn H + 12 δTT H (cid:19) + (cid:16) ξ − τ t (cid:17) v H ∇ · δ u v H + λ H n H ∇ (cid:18) λ ∗ δTT H + µ ∗ δnn H (cid:19) , (3.2b) ∂ t δu i v H = ζ (0) H δu i v H − v H ∇ i (cid:18) δnn H + δTT H (cid:19) + η H mn H (cid:20)(cid:18) η ∗ η ∗ b (cid:19) ∇ i ∇ · δ u v H + η ∗ ∇ δu i v H (cid:21) , (3.2c)Note that in Eqs. (3.2) the deviations { δn, δT, δ u } have been scaled with respect to the HCS quantities { n H , T H ( t ) , v H ( t ) } , where v H ( t ) = r T H ( t ) τ t m (3.3)is the thermal velocity in the HCS. This scaling is non-trivial in the cases of the flow velocity and the temper-ature since the reference state is cooling down and thus v H and T H are time-dependent quantities.As in the purely smooth case ( β = −
1) [4, 5], it isadvisable to introduce the scaled variables s ( t ) = 12 Z t dt ′ ν H ( t ′ ) , ℓ = 12 ν H ( t ) v H ( t ) r , (3.4)where ν H ( t ) is defined by Eq. (2.4b) with T → T H eval-uated in the HCS. The quantity s ( t ) is a measure of thenumber of collisions per particle up to time t , while ℓ measures position in units of the mean free path (no-tice that the ratio ν H /v H is independent of time). A setof Fourier transformed dimensionless variables are thendefined by ρ k ( s ) = δ e n k ( s ) n H , Θ k ( s ) = δ e T k ( s ) T H ( s ) , w k ( s ) = δ e u k ( s ) v H ( s ) , (3.5)where the Fourier transforms { δ e y k ,a ( s ); a = 1 , . . . , } ≡{ δ e n k ( s ) , δ e T k ( s ) , δ e u k ( s ) } are δ e y k ,a ( s ) = Z d ℓ e − ı k · ℓ δy a ( ℓ , s ) . (3.6)Note that in Eq. (3.6) the wave vector k is dimensionless.In terms of the variables (3.4) and (3.5), Eqs. (3.2)become ∂ s ρ k = − ı k · w k , (3.7a) ∂ s Θ k = − ζ ∗ (2 ρ k + Θ k ) + (cid:16) ξ − τ t (cid:17) ı k · w k − k ( λ ∗ Θ k + µ ∗ ρ k ) , (3.7b) ∂ s w k ,i = ζ ∗ w k ,i − ık i ( ρ k + Θ k ) − (cid:20)(cid:18) η ∗ η ∗ b (cid:19) k i k · w k + η ∗ k w k ,i (cid:21) . (3.7c)As expected, the two transversal velocity components w k , ⊥ = w k − ( w k · b k ) b k (orthogonal to the wave vector k ) decouple from the other three modes and hence canbe analyzed independently. Their evolution equation is (cid:18) ∂ s − ζ ∗ + η ∗ k (cid:19) w k , ⊥ = 0 , (3.8)whose solution is w k , ⊥ ( s ) = w k , ⊥ (0) e ω ⊥ ( k ) s , (3.9)where ω ⊥ ( k ) = ζ ∗ − η ∗ k . (3.10)This identifies two shear (transversal) modes analogousto the ones for molecular gases. According to Eq. (3.10),there exists a critical wave number k ⊥ , given by k ⊥ = s ζ ∗ η ∗ , (3.11)separating two regimes: shear modes with k > k ⊥ alwaysdecay, while those with k < k ⊥ grow exponentially. Forpurely smooth spheres, the critical wave number k ⊥ , sm isgiven by Eq. (3.11) with ζ ∗ → ζ ∗ sm and η ∗ → η ∗ sm .The remaining (longitudinal) modes correspond to ρ k ,Θ k , and the longitudinal velocity component of the ve-locity field, w k , || = w k · b k (parallel to k ). These modesare coupled and obey the matrix equation ∂ s δy k ,a ( s ) = M ab ( k ) δy k ,b ( s ) , (3.12)where δy k ,a ( s ) denotes now the set (cid:8) ρ k , Θ k , w k , || (cid:9) and M ( k ) is the square matrix M ( k ) = − ık ζ ∗ + µ ∗ k ζ ∗ + λ ∗ k τ t − ξ ıkık ık η ∗ +3 η ∗ b k − ζ ∗ . (3.13)The longitudinal three modes have the formexp[ ω k ,a ( k ) s ] for a = 1 , ,
3, where { ω k ,a ( k ) } are the eigenvalues of the matrix M ( k ), i.e., they are the solu-tions of the cubic equation ω + A ( k ) ω + B ( k ) ω + C ( k ) = 0 , (3.14)where A ( k ) = (cid:18) λ ∗ η ∗ η ∗ b (cid:19) k , (3.15a) B ( k ) = 5 λ ∗ (cid:18) η ∗ η ∗ b (cid:19) k + h − ξ + τ t ζ ∗ (cid:18) η ∗ η ∗ b − λ ∗ (cid:19) i k − ζ ∗ , (3.15b) C ( k ) = 58 ( λ ∗ − µ ∗ ) k − ζ ∗ k . (3.15c)In the purely smooth case ( β = − ζ → ζ (0) [4]. IV. DISCUSSIONA. Long-wavelength limit
Before considering the general case ( k = 0), it is con-venient to consider first the solutions to Eq. (3.14) inthe long-wavelength limit k → ω ⊥ (0) = ζ ∗ , { ω || ,a (0); a = 1 , , } = {− ζ ∗ , , ζ ∗ } . (4.1)Since two of the eigenvalues are positive (correspondingto growth of the initial perturbation in time), this meansthat some of the solutions are unstable, as expected. Thezero eigenvalue represents a marginal stability solution,while the negative eigenvalue provides a stable solution.The unstable (positive) modes simply emerge from thenormalization of the flow velocity by the time-dependentthermal velocity v H ( t ), which is required to get time-independent coefficients in the linearized hydrodynamicequations after scaling the hydrodynamic fields with re-spect to their homogeneous values.The solution to the cubic equation (3.14) for smallwave numbers k can be written in the form of the se-ries expansion ω k ,a ( k ) = ω (0) k ,a + k ω (2) k ,a + · · · , a = 1 , , , (4.2)where the Euler eigenvalues ω (0) k ,a = ω k ,a (0) are given bythe second equality in Eq. (4.1). Note that only evenpowers in k appear in Eq. (4.2). Substitution of the ex-pansion (4.2) into Eq. (3.14) yields ω (2) k , = 6 − ξ + τ t ζ ∗ − λ ∗ , (4.3a) -0.3-0.2-0.10.00.10.20.30.0 0.2 0.4 0.6-0.3-0.2-0.10.00.10.2 0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6 0.8 ( k ) =0.8, =-1 (a) (b)=0.8, ! -1 (c)=0.8, =-0.5(d) ( k ) =0.8, =0k (e)=0.8, =0.5k (f)=0.8, =1k FIG. 1. Dispersion relations ω ( k ) for the hydrodynamic modes vs the reduced wave number k . The curves correspond to thedegenerate shear mode ω ⊥ (– – –), the heat mode ω k , (– · – · –), and the sound modes ω k , and ω k , (—). Note that when ω k , and ω k , become a complex conjugate pair, only the (common) real part is plotted. The circles denote the critical wave numbers k ⊥ and k k . The coefficient of normal restitution is α = 0 . κ = , while the coefficientsof tangential restitution are (a) β = − β → − β = − .
5, (d) β = 0, (e) β = 0 .
5, and (f) β = 1. ω (2) k , = − ζ ∗ , (4.3b) ω (2) k , = − η ∗ − η ∗ b − τ t − ξ ζ ∗ . (4.3c)Notice that the small wave number solutions (4.3) arerelevant since they are of the same order O ( k ) as theNS hydrodynamic equations. B. Finite wavelength
By extending the usual nomenclature in normal fluids[37] to the case of granular gases, the three longitudinalmodes can be referred to as two “sound” modes ( ω k , and ω k , ) and a “heat” mode ( ω k , ). As happens for smooth-sphere granular gases [4, 5], the shear and heat modes( ω ⊥ and ω k , ) are real for all k . On the other hand,the two sound modes ( ω k , and ω k , ) become a complexconjugate pair of propagating modes in a certain interval k ∗ < k < k ∗∗ , where k ∗ and k ∗∗ depend on α , β , and κ .In particular, k ∗ → β → − ω ( k ) are illustrated in Fig.1 for the case of uniform spheres ( κ = ) at a repre-sentative value of the coefficient of normal restitution( α = 0 . β = −
1) [4, 5], while Fig. 1(b) represents the qua-sismooth limit β → −
1. In Figs. 1(c)–1(f) roughness in-creases from β = − . β = 1. The respective values of ( k ∗ , k ∗∗ ) are (a) (0 . , . , . . , . . , . . , . . , . k < k k , where k k = s ζ ∗ λ ∗ − µ ∗ ) (4.4)can be obtained from Eq. (3.14) by setting ω = 0. Theonly exception is the quasismooth limit ( β → −
1) sincein that case ζ ∗ → k ⊥ and k k tend tozero.For purely smooth spheres ( β = − k k , sm can be obtained from Eq. (4.4) with thereplacements ζ ∗ = ζ ∗ sm , λ ∗ = 2 λ ∗ sm , and µ ∗ = 2 µ ∗ sm . Inthat case, one has k ⊥ , sm > k k , sm if α is higher than athreshold value α th , sm = 0 .
345 [5, 8], as illustrated in Fig.1(a) for α = 0 .
8. Thus, in the purely smooth case, theinstability is driven by the transversal shear mode (vortexinstability) in the inelasticity domain α > α th , sm . Figure1 shows that this is also the case for rough spheres with α = 0 .
8. The property k ⊥ > k k for α = 0 . k ⊥ and k k as functions of β .On the other hand, if α is small enough, there exists anintermediate window of β -values where k k > k ⊥ , so theinstability is driven in that case by the longitudinal heatmode (clustering instability). This is the case shown inFig. 2(b) for α = 0 .
4. If the spheres are uniform ( κ = ), c r i t i c a l w a v e nu m b e r (a) =0.8=0.4 (b) c r i t i c a l w a v e nu m b e r FIG. 2. Plot of the critical wave numbers k ⊥ (– – –) and k k (– · – · –) as functions of β for a reduced moment of inertia κ = with (a) α = 0 . α = 0 .
4. The horizontallines represent the respective values ( k ⊥ , sm and k k , sm ) in thepurely smooth case. the threshold value of α below which it is possible to have k k > k ⊥ is α th = 0 . α th , sm = 0 . k ⊥ and k k for rough spheres, Fig. 2 alsoincludes as horizontal lines those quantities ( k ⊥ , sm and k k , sm ) for purely smooth spheres. As can be seen, one has k ⊥ < k ⊥ , sm and k k < k k , sm in the regions of small andlarge roughness, while the opposite happens for interme-diate roughness. This explains the dual role of roughnesspreviously observed in MD simulations [24]: “high levelsof friction actually attenuate instabilities relative to thefrictionless case, whereas moderate levels enhance insta-bilities compared to frictionless systems.” However, aswill be seen below, the attenuation effect dominates forany level of friction if the inelasticity is high enough.The full combined dependence of k ⊥ , k k , and k ⊥ − k k on both α and β is illustrated in Fig. 3, again for κ = .Figures 3(a) and 3(b) confirm that, in agreement withFig. 2, the critical wave numbers k ⊥ and k k reach theirmaximum values at intermediate roughness ( β ≈
0) fora given value of α . On the other hand, while at β ≈ k k monotonicallyincreases with increasing inelasticity, the transversal crit- - - FIG. 3. Density plot of (a) the critical wave number k ⊥ , (b)the critical wave number k k , and (c) the difference k ⊥ − k k , allin the case of uniform spheres ( κ = ). Two adjacent contourlines correspond to a step (a) ∆ k ⊥ = 0 .
05, (b) ∆ k k = 0 .
1, and(c) ∆( k ⊥ − k k ) = 0 .
05. In panel (c), the thick contour line k ⊥ − k k = 0 divides the plane ( α, β ) into the (outer) regionwhere the transversal shear mode is the most unstable one( k ⊥ > k k ) and the (inner) region where the longitudinal heatmode is the most unstable one ( k k > k ⊥ ). -1.0 -0.5 0.0 0.5 1.0020406080100 L c / =0.5 =0.95 FIG. 4. Plot of the reduced critical length L c /σ as a functionof β for a reduced moment of inertia κ = and a solid volumefraction φ = 0 .
05 with α = 0 .
95 (– – –) and α = 0 . · – · –)[see Eq. (4.6)]. The horizontal lines represent the respectivevalues ( L c, sm /σ ) in the purely smooth case, while the symbolsare MD results at ( α, β ) = (0 . ,
0) (circle) and (0 . , . -1.0 -0.5 0.0 0.5 1.00.00.20.40.60.81.0 A’B’C’DCBA L c
0. As a consequence, and as alreadyanticipated in Fig. 2, the vortex instability (associatedwith k ⊥ ) is preempted by the clustering instability (as-sociated with k k ) in a dome-shaped region with an apexat ( α, β ) = ( α th , β th ) = (0 . , − . -6 -5 -4 -3 -2 -1 T r T ( s ) / T ( ) s T t smooth FIG. 6. Translational and rotational temperatures versus thescaled time variable s [see Eq. (3.4)] for α = 0 . κ = ,and β = − .
99. The initial condition is T r (0) = T t (0). Thedashed line represents the evolution of T t ( s ) in the purelysmooth case ( β = − -1.0 -0.5 0.0 0.5 1.00.00.20.40.60.81.00.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.01.2 L c,M
D th (b)
FIG. 8. Plot of (a) α th and α D , and (b) β th and β D as func-tions of the reduced moment of inertia κ . The vertical dottedlines mark the values κ = 0 .
277 (above which k k is well de-fined for all α and β ), κ = (uniform mass distribution), and κ = (mass concentrated on the spherical surface). C. Critical length
In a system with periodic boundary conditions, thesmallest allowed wave number is 2 π/L , where L is thelargest system length. Thus, we can identify a criticallength L c , such that the system becomes unstable when L > L c . According to the scaling (3.4), the value of L c is given by L c = 4 πv H ν H min { k − ⊥ , k − k } . (4.5)Making use of Eqs. (2.4b) and (3.3), this can be rewrittenas L c σ = 5 π √ π φχ ( φ ) min { k − ⊥ , k − k } . (4.6)Figure 4 shows the β -dependence of L c /σ for α = 0 . α = 0 . κ = ) at a small (but nonzero)solid volume fraction φ = 0 .
05, in which case χ ( φ ) =1 . L c, sm /σ , are represented by horizontal lines. The dualrole of roughness (or friction) is again quite apparent. For small and large roughness the system is less unstable(instability attenuation) than its smooth-sphere counter-part, while the opposite happens (instability enhance-ment) for intermediate roughness. It can be observedthat the central enhancement region is much shorter with α = 0 . α = 0 .
95. Figure 4 also includes val-ues obtained by MD simulations [36] for a representativecase of instability enhancement ( α = 0 . β = 0) anda representative case of instability attenuation ( α = 0 . β = 0 . L c < L c, sm shrinks as inelasticity increases(i.e., as α decreases). Is a threshold value of α reachedbelow which L c > L c, sm for all β ? To address this ques-tion, Fig. 5 displays a phase diagram (assuming uni-form spheres, i.e., κ = ) where the enhancement re-gion (or “phase”) L c < L c, sm is separated from the at-tenuation region (or “phase”) L c > L c, sm by the lo-cus line L c = L c, sm . The peculiar shape of the lo-cus is due to the change from the shear mode to theheat mode as the most unstable one as inelasticity in-creases. More specifically, in the segments A–B andA’–B’ one has ( L c , L c, sm ) ∝ ( k − ⊥ , k − ⊥ , sm ), while in thesegments B–C and B’–C’ the situation is ( L c , L c, sm ) ∝ ( k − k , k − ⊥ , sm ) and ( L c , L c, sm ) ∝ ( k − ⊥ , k − k , sm ), respec-tively. Finally, in the segment C–D–C’, ( L c , L c, sm ) ∝ ( k − k , k − k , sm ). The coordinates of the relevant points are( α B , β B ) = (0 . , − . α B ′ , β B ′ ) = (0 . , . α C , β C ) = (0 . , − . α C ′ , β C ′ ) = (0 . , . α D , β D ) = (0 . , . α B < α th ,while α B ′ = α C = α th , sm . As Fig. 5 clearly shows, thelocus L c = L c, sm presents a minimum vertex at point D.Therefore, if α < α D , the frictional system is always lessunstable than the frictionless system. D. On the quasismooth limit
It seems paradoxical that, as illustrated in Figs. 1(a),1(b), 2, 4, and 5, the HCS is stable in the limit β → − L > L c, sm ) in the purelysmooth case. The explanation lies in the fact that theHCS with β = − β & −
1. In the former case, the rotational degrees offreedom do not play any role at all and the translationaltemperature decays with the cooling rate (2.8). On theother hand, in the quasismooth case the rotational andtranslational degrees of freedom are coupled, the tem-perature ratio is approximately given by Eq. (2.6), andboth temperatures decay with a much smaller coolingrate given by Eq. (2.7).The interesting question is whether the transientregime prior to the HCS for β & − et al. [30], the initial de-cay of the translational temperature T t of nearly smoothparticles is dominated by the coefficient of normal resti-0tution α and reaches a very small value (relative to T r )before the asymptotic HCS is attained. This is illustratedby Fig. 6, which shows the evolution of T t and T r for α = 0 . κ = , and β = − .
99, as obtained by numer-ically solving the coupled set of equations ˙ T t = − ζ t T t and ˙ T r = − ζ r T r , starting from an initial condition ofenergy equipartition, i.e., T r (0) = T t (0). Here, ζ t and ζ r , with ζ = ( ζ t T t + ζ r T r ) / ( T t + T r ), are the collisionalrates of change of T t and T r , respectively [23]. It canbe observed that up to s ≈
30 collisions per particle theevolution of T t is practically indistinguishable from thatof the smooth-sphere system. Therefore, a perturbationwith a wavelength larger than L c, sm during this transientregime will develop vortex or cluster instabilities. In or-der to prevent instabilities in the quasismooth limit, onewould need to fine-tune the initial state with a tempera-ture ratio T r (0) /T t (0) ∼ θ . E. Comparison with the approach of Mitrano et al. [24]
According to the heuristic approach of Mitrano et al. [24], the critical wave numbers are obtained by startingfrom the smooth-sphere values k ⊥ , sm and k k , sm , and sim-ply replacing the smooth-sphere cooling rate ζ ∗ sm [see Eq.(2.8)] by the rough-sphere cooling rate ζ ∗ [see Eq. (2.5b)]: k ⊥ , M = s ζ ∗ η ∗ sm , k k , M = s ζ ∗ λ ∗ sm − µ ∗ sm ) . (4.7)The associated critical length is then L c, M σ = 5 π √ π φχ ( φ ) min { k − ⊥ , M , k − k , M } . (4.8)Figure 7(a) shows that, at a given value of β , L c, M 345 andin the lower region ( L c, M < L c ). This is illustrated inFig. 7(b) for the quasismooth limit β → − 1, the mediumroughness case β = 0, and the completely rough case β = 1. F. Influence of the moment of inertia All the results displayed in Figs. 1–7 correspond tospheres with a uniform mass distribution ( κ = ). Nev-ertheless, at least at a quantitative level, the results areexpected to be influenced by the value of the reduced mo-ment of inertia κ . In this respect, it must be mentioned that the critical longitudinal wave number (4.4) turns outto be well defined for all α and β only if κ > . λ ∗ < µ ∗ below a certain κ -dependent value of α (with amaximum α = 0 . κ → 0) and for a certain in-terval of values of β . A similar situation takes place inthe purely smooth case, where λ ∗ sm < µ ∗ sm if α < . α and/or κ . In order to assess the influence of κ , the threshold val-ues α th (below which one may have k k > k ⊥ ) and α D (be-low which L c > L c, sm for all β ) are plotted as functionsof κ ≥ . 25 in Fig. 8(a). Figure 8(b) does the same butfor the other coordinates β th and β D . As the mass of thespheres becomes more concentrated in the inner layers(i.e., as κ decreases), α th tends to increase (although itpresents a maximum value α th = 0 . κ = 0 . α D monotonically decreases. As for the associatedvalues β th and β D of the coefficient of tangential restitu-tion, both decrease as κ decreases (especially in the caseof β D ), but their magnitudes are smaller than about 0 . β ≈ V. CONCLUSIONS In this paper we have undertaken a rather completestudy of the linear stability conditions of the HCS of adilute granular gas modeled as a system of (identical) in-elastic and frictional hard spheres. Inelasticity and sur-face friction are characterized by constant coefficients ofnormal ( α ) and tangential ( β ) restitution. The analysisis based on the NS hydrodynamic equations, linearizedaround the time-dependent HCS solution. The most rel-evant outcome is the determination of the critical length L c , such that the system is linearly unstable for sizeslarger than L c . The novel aspect of our study, not ac-counted for in previous ones [24], is the use of the detailednonlinear dependence of the NS transport coefficients onboth α and β , as well as on the reduced moment of in-ertia κ [23]. This allows one to explore the impact ofroughness on the critical length L c without a priori anyrestriction on α , β , and κ .As in the purely smooth case [4, 5], two of the five hy-drodynamic modes (the two longitudinal “sound” modes)are stable. On the other hand, a doubly degeneratetransversal (shear) mode and a longitudinal “heat” modebecome unstable for (reduced) wave numbers smallerthan certain critical values k ⊥ and k k , respectively. Ingeneral, the instability is driven by the transversal mode,i.e., k ⊥ > k k . On the other hand, if α is small enough( α < α th = 0 . 536 in the case of uniform spheres) and β lies inside an α -dependent interval around β ≈ 0, thesituation is reversed, i.e., the heat mode is the mostunstable one. As a consequence, the critical length L c ∝ min { k − ⊥ , k − k } exhibits a nontrivial dependence on1 α , β , and κ . Comparison of the theoretical predictionsfor L c against preliminary MD simulations [36] shows anexcellent agreement.An interesting point is the comparison between thecritical length L c of a gas of rough spheres and its cor-responding counterpart, L c, sm , of a gas of purely smoothspheres. One could naively expect that the existence offriction would enhance the instability of the HCS, namely L c < L c, sm , for a common value of the inelasticity pa-rameter α . Nevertheless, this is not the general case.As illustrated in Fig. 5 for uniform spheres, an attenua-tion effect is present at sufficiently low or sufficiently highlevels of friction. This dual role of friction was alreadyobserved by Mitrano et al. [24] in MD simulations andin a simple kinetic theory description for dense gases.On the other hand, the enhancement middle region ofroughness disappears if the inelasticity is large enough( α < α D = 0 . 218 for uniform spheres). With respect tothe influence of the moment of inertia, our results showthat, as the mass of each sphere concentrates nearer thesurface (i.e., as κ increases), the values of α th and α D decrease and increase, respectively.It is worthwhile mentioning that the cooling rate of thegas of rough spheres, as compared with that of purelysmooth spheres, already exhibits a dual behavior: dissi-pation of energy is enhanced by friction only if α is largerthan a certain threshold value ( α > . 401 for uniformspheres) and β lies in an α -dependent interval around β ≈ 0. Otherwise, friction attenuates energy dissipation.This explains the good performance of the simple kinetictheory proposed by Mitrano et al. [24], where the criticalwave numbers k ⊥ and k k are obtained by assuming thesame expressions as for the purely smooth system, exceptfor the replacement of the cooling rate.A subtle point is the suppression of clustering and vor-tex formation in the quasismooth limit ( β → − α tends to 1as the impact velocity decreases show that structure for-mation occurs in free granular gases only as a transientphenomenon, whose duration increases with the systemsize [38]. However, the experimental measurement of α at very small impact velocities is very challenging [39].Some independent experiments [40, 41] provide evidenceon a sharp decrease of α at small impact velocities, pos-sibly due to van der Waals attraction at relatively lowsurface energies for typical grain materials.Finally, we hope that the results presented in this workwill stimulate the performance of computer simulationsto further assess their practical usefulness. ACKNOWLEDGMENTS We want to thank Peter P. Mitrano for making the sim-ulation data included in Fig. 4 available to us. V.G. andA.S. acknowledge the financial support of the Ministeriode Econom´ıa y Competitividad (Spain) through GrantNo. FIS2016-76359-P, partially financed by “Fondo Eu-ropeo de Desarrollo Regional” funds. The research ofG.M.K is supported by Conselho Nacional de Desenvolvi-mento Cient´ıfico e Tecnol´ogico (CNPq), Brazil, throughGrant No. 303251/2015-8. Appendix: Explicit expressions for the NS transportand cooling-rate coefficients In this Appendix, the expressions for the NS trans-port coefficients ( η ∗ , η ∗ b , λ ∗ , µ ∗ ) and the NS cooling ratecoefficient ( ξ ) are explicitly given.The reduced coefficients η ∗ and η ∗ b associated with thepressure tensor are given by Eq. (2.11a) with ν ∗ η = ( e α + e β )(2 − e α − e β ) + e β θ κ , (A.1a) γ E = 23 1Ξ t − Ξ r − ζ ∗ , (A.1b)where e α = 1 + α , e β = 1 + β κ κ , (A.2a)Ξ t = 58 τ r h − α + (1 − β ) κ κ − κ θ − (cid:18) β κ (cid:19) i , (A.2b)Ξ r = 58 τ t β κ (cid:20) θ − 23 (1 − β ) + κ θ − 5) 1 + β κ (cid:21) . (A.2c)The reduced coefficients λ ∗ and µ ∗ associated with theheat flux are given by Eq. (2.11b) with γ A t = Z r − Z t − ζ ∗ ( Y t − ζ ∗ ) ( Z r − ζ ∗ ) − Y r Z t , (A.3a) γ A r = Y t − Y r − ζ ∗ ( Y t − ζ ∗ ) ( Z r − ζ ∗ ) − Y r Z t , (A.3b) γ B t = ζ ∗ γ A t (cid:0) Z r − ζ ∗ (cid:1) − γ A r Z t (cid:0) Y t − ζ ∗ (cid:1) (cid:0) Z r − ζ ∗ (cid:1) − Y r Z t , (A.3c) γ B r = ζ ∗ γ A r (cid:0) Y t − ζ ∗ (cid:1) − γ A t Y r (cid:0) Y t − ζ ∗ (cid:1) (cid:0) Z r − ζ ∗ (cid:1) − Y r Z t . (A.3d)2In Eqs. (A.3), we have introduced the quantities Y t = 4112 (cid:16)e α + e β (cid:17) − (cid:16)e α + e β (cid:17) − e α e β − θ e β κ , (A.4a) Y r = 2536 e βκ − e βθ − e βκ ! , (A.4b) Z t = − θ e β κ , (A.4c) Z r = 56 (cid:16)e α + e β (cid:17) + 518 e βκ − e βκ − e β − e α ! . (A.4d)Finally, the first-order contribution ξ to the coolingrate is given by Eq. (2.12), where γ E is given by Eq.(A.1b) and Ξ isΞ = 516 τ t τ r (cid:20) − α + (1 − β ) (cid:18) θ − 51 + κ (cid:19)(cid:21) . (A.5) [1] W. D. Fullmer and C. M. Hrenya, “The clustering insta-bility in rapid granular and gas-solid flows,” Annu. Rev.Fluid Mech. , 485–510 (2017).[2] I. Goldhirsch and G. Zanetti, “Clustering instabilityin dissipative gases,” Phys. Rev. Lett. , 1619–1622(1993).[3] S. McNamara, “Hydrodynamic modes of a uniform gran-ular medium,” Phys. Fluids A , 3056–3070 (1993).[4] J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, “Hy-drodynamics for granular flow at low density,” Phys. Rev.E , 4638–4653 (1998).[5] V. Garz´o, “Instabilities in a free granular fluid describedby the Enskog equation,” Phys. Rev. E , 021106(2005).[6] J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, “Insta-bility and spatial correlations in a dilute granular gas,”Phys. Fluids , 2976–2982 (1998).[7] P. P. Mitrano, S. R. Dahl, D. J. Cromer, M. S. Pacella,and C. M. Hrenya, “Instabilities in the homogeneouscooling of a granular gas: A quantitative assessmentof kinetic-theory predictions,” Phys. Fluids , 093303(2011).[8] P. P. Mitrano, V. Garz´o, A. M. Hilger, C. J. Ewasko, andC. M. Hrenya, “Assessing a modified-Sonine kinetic the-ory for instabilities in highly dissipative, cooling granulargases,” Phys. Rev. E , 041303 (2012).[9] V. Garz´o and J. W. Dufty, “Hydrodynamics for a granu-lar mixture at low density,” Phys. Fluids , 1476–1490(2002).[10] V. Garz´o, J. M. Montanero, and J. W. Dufty, “Mass andheat fluxes for a binary granular mixture at low density,”Phys. Fluids , 083305 (2006).[11] V. Garz´o, J. W. Dufty, and C. M. Hrenya, “Enskog the-ory for polydisperse granular mixtures. I. Navier-Stokesorder transport,” Phys. Rev. E , 031303 (2007).[12] V. Garz´o, C. M. Hrenya, and J. W. Dufty, “Enskog the-ory for polydisperse granular mixtures. II. Sonine poly-nomial approximation,” Phys. Rev. E , 031304 (2007).[13] J. J. Brey and M. J. Ruiz-Montero, “Shearing instabilityof a dilute granular mixture,” Phys. Rev. E , 022210(2013).[14] P. P. Mitrano, V. Garz´o, and C. M. Hrenya, “Instabilitiesin granular binary mixtures at moderate densities,” Phys.Rev. E , 020201(R) (2014).[15] A. Bodrova and N. Brilliantov, “Self-diffusion in granular gases: an impact of particles’ roughness,” Granul. Matter , 85–90 (2012).[16] V. Garz´o and J. W. Dufty, “Dense fluid transport for in-elastic hard spheres,” Phys. Rev. E , 5895–5911 (1999).[17] J. F. Lutsko, “Transport properties of dense dissipa-tive hard-sphere fluids for arbitrary energy loss models,”Phys. Rev. E , 021306 (2005).[18] V. K. Gupta, P. Shukla, and M. Torrilhon, “Higher-ordermoment theories for dilute granular gases of smooth hardspheres,” J. Fluid Mech. , 451–501 (2018).[19] J. T. Jenkins and M. W. Richman, “Kinetic theory forplane flows of a dense gas of identical, rough, inelastic,circular disks,” Phys. Fluids , 3485–3494 (1985).[20] I. Goldhirsch, S. H. Noskowicz, and O. Bar-Lev, “Nearlysmooth granular gases,” Phys. Rev. Lett. , 068002(2005).[21] I. Goldhirsch, S. H. Noskowicz, and O. Bar-Lev, “Hy-drodynamics of nearly smooth granular gases,” J. Phys.Chem. B , 21449–21470 (2005).[22] C. K. K. Lun, “Kinetic theory for granular flow of dense,slightly inelastic, slightly rough spheres,” J. Fluid Mech. , 539 (1991).[23] G. M. Kremer, A. Santos, and V. Garz´o, “Transportcoefficients of a granular gas of inelastic rough hardspheres,” Phys. Rev. E , 022205 (2014).[24] P. P. Mitrano, S. R. Dahl, A. M. Hilger, C. J. Ewasko,and C. M. Hrenya, “Dual role of friction in granular flows:attenuation versus enhancement of instabilities,” J. FluidMech. , 484–495 (2013).[25] A. Santos, G. M. Kremer, and V. Garz´o, “Energy pro-duction rates in fluid mixtures of inelastic rough hardspheres,” Prog. Theor. Phys. Suppl. , 31–48 (2010).[26] F. Vega Reyes, A. Lasanta, A. Santos, and V. Garz´o,“Energy nonequipartition in gas mixtures of inelasticrough hard spheres: The tracer limit,” Phys. Rev. E ,052901 (2017).[27] N. V. Brilliantov and T. P¨oschel, Kinetic Theory of Gran-ular Gases (Oxford University Press, Oxford, 2004).[28] K. K. Rao and P. R. Nott, An Introduction to GranularFlow (Cambridge University Press, Cambridge, 2008).[29] A. Goldshtein and M. Shapiro, “Mechanics of collisionalmotion of granular materials. Part 1. General hydrody-namic equations,” J. Fluid Mech. , 75–114 (1995).[30] S. Luding, M. Huthmann, S. McNamara, and A. Zip-pelius, “Homogeneous cooling of rough, dissipative par- ticles: Theory and simulations,” Phys. Rev. E ,3416–3425 (1998).[31] N. F. Carnahan and K. E. Starling, “Equation of state fornonattracting rigid spheres,” J. Chem. Phys. , 635–636(1969).[32] N. V. Brilliantov, T. P¨oschel, W. T. Kranz, and A. Zip-pelius, “Translations and rotations are correlated in gran-ular gases,” Phys. Rev. Lett. , 128001 (2007).[33] W. T. Kranz, N. V. Brilliantov, T. P¨oschel, and A. Zip-pelius, “Correlation of spin and velocity in the homoge-neous cooling state of a granular gas of rough particles,”Eur. Phys. J. Spec. Top. , 91–111 (2009).[34] A. Santos, G. M. Kremer, and M. dos Santos, “So-nine approximation for collisional moments of granulargases of inelastic rough spheres,” Phys. Fluids , 030604(2011).[35] A. Santos, “Homogeneous free cooling state in binarygranular fluids of inelastic rough hard spheres,” AIPConf. Proc. , 128–133 (2011). [36] P. P. Mitrano, Private communication.[37] P. R´esibois and M. de Leener, Classical Kinetic Theoryof Fluids (John Wiley and Sons, New York, 1977).[38] N. Brilliantov, C. Salue˜na, T. Schwager, and T. P¨oschel,“Transient structures in a granular gas,” Phys. Rev. Lett. , 134301 (2004).[39] Y. Duan and Z.-G. Feng, “Incorporation of velocity-dependent restitution coefficient and particle surface fric-tion into kinetic theory for modeling granular flow cool-ing,” Phys. Rev. E , 062907 (2017).[40] C. M. Sorace, M. Y. Louge, M. D. Crozier, and V. H. C.Law, “High apparent adhesion energy in the breakdownof normal restitution for binary impacts of small spheresat low speed,” Mech. Res. Commun. , 364–368 (2009).[41] Y. Grasselli, G. Bossis, and G. Goutallier, “Velocity-dependent restitution coefficient and granular cooling inmicrogravity,” Europhys. Lett.86