Green-Kubo representation of the viscosity of granular gases

Green-Kubo representation of the viscosity of granular gases J. Javier Brey Área de Física Teórica, Universidad de Sevilla. Apartado de Correos 1065, E-41080, Sevilla, Spain Abstract. The Green-Kubo representation of the shear viscosity of a dilute gas of inelastic hard disks derived from the Boltzmann equation is evaluated by means of the direct simulation Monte Carlo method. The relationship between the oneparticle dynamics of the original expression and the N-particle dynamics needed for the simulation is analyzed. The presence of velocity correlations in the latter is identified and their possible implications are discussed. INTRODUCTION The prototypical model for granular gases is an ensemble of smooth inelastic hard spheres or disks with a constant coefficient of normal restitution α . For this model, the methods of kinetic theory have been applied and, in particular, the Boltzmann equation has been extended [1, 2]. Then, by using a modification of the Chapmann-Enskog procedure, hydrodynamicequations have been derived [3, 4]. As the reference state, the isotropic and homogeneouscooling state (HCS) of a granular system is used. Due to the energy dissipation, the HCS is not stationary, but its energy decreases monotonically. Hydrodynamics has been employed with success to describe the behavior of granular flows. As for molecular systems, the Chapmann-Enskogmethod leads to complicated differential equations that have to be solved in some approximations.It has been shown that the aforementionedequations are equivalentto a representation of the transport coefficients that has a similar structure to the Green-Kubo relations for elastic systems, although with significant differences [5]. Equivalent expressions have been derived by identifying the hydrodynamic part of the spectrum of the linearized Boltzmann equation and assuming it dominates for long long times and wavelengths [6, 7]. Here, an evaluation of the Green-Kubo relation for the shear viscosity by means of N-particle simulations, namely the direct simulation Monte Carlo (DSMC) method [8], is presented. This requires going from the one-particle effective dynamics, as defined by a linearized Boltzmann operator, to a supposed equivalent N-particle dynamics. The relationship between both descriptions will turn out to be important for the analysis of the results. The existence an exact mapping of the HCS onto a steady state will be exploited as the basis of the simulation method [9, 10]. THE GREEN-KUBO EXPRESSION FOR THE SHEAR VISCOSITY By using the Chapman-Enskog method or eigenfunctions expansions [5, 6], the shear viscosity η of a dilute granular gas of Nsmooth inelastic hard disks (d=2) or spheres (d=3) of mass mand diameter σ can be expressed in the form η (T) = nm`v0(T)Z∞ 0dse−s ζ 0/2Zdc1 χ HCS(c1)∆xy(c1,s)Φ2,xy(c1).(1) Here nis the density, Tthe temperature, v0(T) = (2kBT/m)1/2with kBthe Boltzmann constant, `= (n σ d−1)−1, and ∆xy(c) = cxcy,Φ2,xy(c) = −cx ∂ ln χ HCS(c) ∂ cy.(2) The function χ HCS(c1)is defined from the solution of the Boltzmann equation describing the HCS, fHCS(v1,t), as [1] fHCS(v1,t) = nv−d 0[T(t)] χ HCS(c1),c1=v1 v0[T(t)] ,(3) 815 CP762, Rarefied Gas Dynamics: 24th International Symposium, edited by M. Capitelli © 2005 American Institute of Physics 0-7354-0247-7/05/$22.50 Report Documentation Page Form Approved OMB No. 0704-0188 Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 1. REPORT DATE 13 JUL 2005 2. REPORT TYPE N/A 3. DATES COVERED - 4. TITLE AND SUBTITLE Green-Kubo representation of the viscosity of granular gases 5a. CONTRACT NUMBER 5b. GRANT NUMBER 5c. PROGRAM ELEMENT NUMBER 6. AUTHOR(S) 5d. PROJECT NUMBER 5e. TASK NUMBER 5f. WORK UNIT NUMBER 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) Área de Física Teórica, Universidad de Sevilla. Apartado de Correos 1065, E-41080, Sevilla, Spain 8. PERFORMING ORGANIZATION REPORT NUMBER 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 11. SPONSOR/MONITOR’S REPORT NUMBER(S) 12. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release, distribution unlimited 13. SUPPLEMENTARY NOTES See also ADM001792, International Symposium on Rarefied Gas Dynamics (24th) Held in Monopoli (Bari), Italy on 10-16 July 2004. 14. ABSTRACT 15. SUBJECT TERMS 16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF ABSTRACT UU 18. NUMBER OF PAGES 6 19a. NAME OF RESPONSIBLE PERSON a. REPORT unclassified b. ABSTRACT unclassified c. THIS PAGE unclassified Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std Z39-18 and it is an isotropic function of the vector c1. Finally, the time dependence of ∆xy is given by ∆xy(c1,s) = es¯ Λc(c1)∆xy(c1),(4) ¯ Λc(c1) = Zdc2 χ HCS(c2)Zdb σ Θ(c12 ·b σ )c12 ·b σ [b σ (c1,c2)−1](1+P12)+ ζ 0 2c1· ∂ ∂ c1,(5) where c12 =c1−c2,Θis the Heaviside step function, b σ is the unit vector pointing from the center of particle 2 to the center of particle 1 at contact, the operator P12 interchanges the subindexes 1 and 2 to its right, and b σ (c1,c2)is an operator replacing all the velocities c1and c2appearing to its right by the postcollisional values c0 1and c0 2given by c0 1≡b σ c1=c1−1+ α 2(b σ ·c12)b σ ,c0 2≡b σ c2=c1+1+ α 2(b σ ·c12)b σ .(6) The last term on the right hand side of Eq. (5) contains the time-independent reduced cooling rate of the HCS, ζ 0, defining the evolution of the temperature in that state, ∂ tT(t) = − ζ HCS(T)T, ζ 0=` ζ HCS v0(T).(7) Equation (1) differs from the standard Green-Kubo expression for dilute molecular systems in several aspects: •The average is taken over the velocity distribution of the HCS and not over the Maxwellian characterizing the equilibrium state. Both differ for α <1. •The time correlation function couples the transversal momentum flux ∆xy and a “modified flux” Φ2,xy. •The dynamics includes an accelerating streaming between collisions. •The time integral contains, in addition to the correlation function, an exponentially decreasing in time factor. To render Eq. (1) suitable for evaluation by means of N-particle simulations, some issues must be addressed. First, it presents the technical difficulty that Λcinvolves the cooling rate ζ 0that, therefore, must be known a priori. Consequently, it is useful to make a change of scale [9, 10], τ = ζ 0 2 ω 0s,w=2 ω 0` ζ 0c(8) with ω 0being an arbitrary time-independent dimensionless frequency. Then, Eq. (1) transforms into η (T) = nmv0(T) ev0,stNZ∞ 0dte− ω 0th∆xy(v,t)Φ2,xy (v/ev0,st)ist .(9) Here, we have renamed the time and velocity variables, ev0,st ≡(2kBe Tst/m)1/2where e Tst =2m kB ω 0` ζ 02 ,(10) and the angular brackets denote average defined by ha(v,t)b(v)ist ≡ZdrZdve fst(v)a(v,t)b(v),e fst(v) = nev−d 0,st χ HCS v ev0,st .(11) Moreover, the time dependence is now given by a(v,t) = et¯ Λst(v)a(v)(12) with ¯ Λst(v1) = σ d−1Zdv2e fst(v2)Zdb σ Θ(v12 ·b σ )v12 ·b σ [b σ (v1,v2)−1](1+P12)+ ω 0v1· ∂ ∂ v1.(13) 816 N-PARTICLE REPRESENTATION In Eq. (9), the dynamics is defined by means of the linear operator e Λst. Before computating it by N-particle simulation methods, it must be rewritten in terms of a well defined particle dynamics. The structure of e Λst suggests to introduce an acceleration streaming between collisions, ∂ ∂ tRi(t) = Vi(t), ∂ ∂ tVi(t) = ω 0Vi(t),(14) while the effect of a collision between particles iand jis to instantaneously modify their velocities accordingly with the same rules as in the original dynamics and given in Eqs. (6), i.e. Vi→V0 i≡b σ Vi=Vi−1+ α 2(b σ ·Vij)b σ ,Vj→V0j≡b σ Vj=Vj+1+ α 2(b σ ·Vij)b σ .(15) For this dynamics, the Liouville equation can be obtained. Then, by means of any of the standard procedures, the Boltzmann equation can be derived in the low density limit. It has the form  ∂ ∂ t+ ω 0 ∂ ∂ v·v+v· ∂ ∂ rf(r,v,t) = J[f,f],(16) where J[f,f]is the same Boltzmann collision operator as in the original dynamics. It is now easily verified that e fst(v) is a steady solution of this equation [10]. Let us assume that the Liouville equation has a stationary solution ρ st(Γ), where Γdenotes a point in phase space. A sufficient condition for this is that the distribution function of the HCS in the actual, original dynamics has the scaling property ρ HCS(Γ,t) = (v0[T(t)])−dN ρ ∗ HCS ({Rij,Vi/v0[T(t)]}).(17) Consider time correlation functions of the form CAB,st(t) = hA(t)BiN,st −hAiN,sthBiN,st,(18) where hAiN,st =ZdΓ ρ st(Γ)A(Γ),hBiN,st =ZdΓ ρ st(Γ)B(Γ),(19) hA(t)BiN,st =ZdΓ ρ st(Γ)A(Γ,t)B(Γ),(20) with A(Γ) = N ∑ i=1a(Vi),B(Γ) = N ∑ i=1b(Vi).(21) The dynamical variable A(Γ,t)≡A[Γ(t)] is generated from A(Γ)through the dynamics defined above.The low density limit of CAB,st can be obtained by using the same procedures as for molecular systems [10]. The hypothesis needed are the same as those required to derive the Boltzmann equation [11], plus the additional assumption that velocity correlations are negligible in the HCS at low densities. Then, it is found that CAB,st is given by Eq. (11). This provides a direct way of evaluating the latter expression, and therefore the shear viscosity, by means of N-particle simulations, since the dynamics defined by Eqs. (14) and (15) is easily implemented. Of course, the equality CAB,st =ha(v,t)b(v)ist (22) only holds at low enough densities and as long as velocity correlations effects are negligible in the HCS. The DSMC method is a very useful tool to generate the low density dynamics of a system [8, 12]. This method is designed as a real N-particle dynamics simulation of a low density gas, consequently providing its complete dynamical description. 817 012345 v −6 −4 −2 0 (ln Χ)’ FIGURE 1. Plot of (lnX)0≡ ∂ ln χ HCS/ ∂ cas a function of c=v/ev0,st for α =0.6. The circles are the numerical derivative of the simulation results and the solid line the fitted function. The velocity is measured in the units defined in the text. SIMULATION RESULTS We have performedDSMC simulations of a system of N=104inelastic hard disks. Since we are dealing with position independent quantities in a homogeneous state, it is enough to consider one cell in configurationspace. This allows to increase thestatistics and alsoavoids that the system spontaneouslydevelopsspatial inhomogeneities.This is important since the HCS is unstable with respect to spatial long wavelength perturbations [13, 14]. We introduce a dimensionless reduced shear viscosity η ∗by scaling the viscosity with its elastic limit in the first Sonine approximation η 0, η ∗= η (T) η 0(T), η 0=1 2 σ mkBT π 1/2 .(23) Then, Eq. (9) yields η ∗=2√2 π `ev0,st Z∞ 0dtJ η (t)e− ω 0t,(24) with J η (t) = 1 Nh∆xy(v,t)Φ2,xy(v/ev0,st)ist .(25) Let us remind that what is actually computed in the N-particle description is J0 η (t) = 1 N N ∑ i N ∑ jh∆xy(vi,t)Φ2,xy(vj/ev0,st)iN,st .(26) The simulations show that the system reaches, after a transient period, a stationary state with a temperature e Tst. Then, we measured the velocity function ∂ ln χ HCS(c)/ ∂ c,c=v/ev0,st, to determine the modified flux Φ2,xy(c), defined in Eq. (2). For this, the range of vwas partitioned into small non-overlapping bins and the frequency distribution was built. Afterwards, the numerical derivative was computed. An example of the results, for α =0.6, is given in Fig. 1. Here and in the following, the unit of mass is m, the unit of length is `= (n σ )−1, and the unit of time is `[2kBe T(0)/m]−1/2, where e T(0)is the initial temperature. Moreover, we set kB=1, implying that in our units it is e T(0) = 1/2. For large velocities, ∂ ln χ HCS(c)/ ∂ ctends to a constant, reflecting the exponential decay of χ HCS(c)[15]. Once the function Φ2,xy has been determined for each of the values of α considered, the correlation function J0 η (t) has been measured. It is seen to decay exponentially in time, at least over more than two decades, for all the values of the restitution coefficient [16]. Then, J0 η (t)was fitted to an exponential, J0 η (t) = J0 η (0)e− λη t, and the fitted parameters 818 0.1 0.3 0.5 0.7 0.9 α 1.0 1.2 1.4 1.6 η∗ FIGURE 2. The reduced coefficient of shear viscosity η ∗for a system of inelastic disks as a function of the coefficient of normal restitution α . The symbols are from the DSMC method and the solid line is the first Sonine approximation. J0 η (0)and λη were determined. In terms of them, Eq. (24) reads η ∗=2√2 π J0 η (0) λη + ω 0.(27) The values of η ∗obtained in this way are shown in Fig. 2. A good agreement is observed with the results obtained in the first Sonine approximation [17], although a slightly different behavior is also clearly identified. It is interesting to consider the initial value of J0 η (t)and decompose it in the form J0 η (0) = J0(1) η (0)+J0(2) η (0),(28) where J0(1) η (0) = 1 N N ∑ ih∆xy(vi)Φ2,xy(vi/ev0,st)iN,st ,J0(2) η (0) = 1 N N ∑ i N ∑ j6=ih∆xy(vi)Φ2,xy(vj/ev0,st)iN,st .(29) The first component, J0(1) η (0)is identically the same as J η (0), where J η (t)was defined in Eq. (25). It is independent of the velocity correlations present in the system. In fact, its exact value can be easily obtained, and for the choice of ω 0used in the simulations and the units we are employing, it is J0(1) η (0) = J η (0) = 1/2 for all α [16]. The other part, J(2) η (0), contains the effect of velocity correlations and was assumed to be negligible. The simulation results for J0 η (0)are plotted in Fig. 3. For α <0.6, contributionsfrom J0(2) η become significant, indicating the presence of velocity correlations in the system. This result raises some fundamental questions. Are the observed correlations an artifact of the DSMC method? We believe they are not, since velocity correlations in the HCS also have been recently observed and measured in molecular dynamics simulations [18]. Why there is such a good agreement between the first Sonine approximation predictions and the simulation results for the shear viscosity, even from strong dissipation? A possible answer is that the Sonine expansion is introduced, and truncated, both in the initial conditions and in the dynamics of the fluxes, leading to some kind of compensation. For small values of α , which Green-Kubo expression is the right one, that with J η (t)or that with J0 η (t)?, or are both wrong? And, finally, can it be concluded from these results that the Boltzmann equation should be somehow revised for strong inelasticity? All these points clearly deserve further analysis. Although in this paper we have restricted ourselves to the shear viscosity, the Green-Kubo expressions for the other two transport coefficients, associated with the heat flux, can be computed in a similar way. For not too strong dissipation, the simulation results agree quite well with the first Sonine approximation predictions. Nevertheless, the discrepancies grow rather fast as α decreases below α ≈0.7. This might be, at least partially, due to the presence 819 0 0.2 0.4 0.6 0.8 1 α 0.45 0.5 0.55 0.6 FIGURE 3. DSMC results for the correlation function J0 η (0)(filled circles) and its diagonal part J0(1) η (0)(empty circles). of velocity correlations in the DSMC results. A detailed discussion of these transport coefficients will be published elsewhere [16]. ACKNOWLEDGMENTS We acknowledge financial support from the Ministerio de Ciencia y Tecnología (Spain) through Grant No. BFM2002- 00303 (partially financed by FEDER funds) REFERENCES 1. Goldsthein M. and Shapiro M., J. Fluid Mech.,282, 75-114 (1995). 2. Brey J.J., Dufty J.W., and Santos A., Physica,87, 1051-1066 (1997). 3. Brey J.J., Dufty J.W., Kim C.S., and Santos A., Phys. Rev. E,58, 4638-4653 (1998). 4. Sela N. and Goldhirsch I., J. Fluid Mech.,361, 41-74 (1998). 5. Dufty J.W. and Brey J.J., Physica,109, 433-448 (2002). 6. Brey J.J., Dufty J.W., and Ruiz-Montero M.J., “Linearized Boltzmann Equation and Hydrodynamics for Granular gases”, in Granular Gas Dynamics, edited by T. Pöschel and N. Brilliantov, Springer-Verlag, Berlin, 2003, pp. 225-246. 7. Dufty J.W. and Brey J.J., Phys. Rev. E,68, 030302-1-4 (2003). 8. Bird G., Molecular gas Dynamics and the Direct Simulation of Gas Flows, Clarendon Press, Oxford (1994). 9. Lutsko J., Phys. Rev. E,63, 061211-1-14 (2001) 10. Brey J.J., Ruiz Montero M.J., and Moreno F., Phys. Rev. E,69, 051303-1-13 (2004). 11. Ernst M.H. and Cohen E.G.D., J. Stat. Phys.,25, 153-180 (1981). 12. García A., Numerical Methods for Physics, Prentice Hall, Englewood Cliffs, NJ (2000). 13. Goldhirsch I. and Zanetti G. Phys. Rev. Lett.,70, 1619-1622 (1993). 14. Brey J.J., Ruiz-Montero M.J., and Cubero D. Phys. Rev. E,54, 3664-3671 (1996). 15. Esipov S.E. and Pöschel T., J. Stat. Phys.,86, 1385-1395 (1997). 16. Brey J.J. and Ruiz-Montero M.J., unpublished. 17. Brey J.J. and Cubero D., “Hydrodynamic Tranport Coefficients of Granular Gases”, in Granular Gases, edite by T. Pöschel and S. Luding, Springer-Verlag, Berlin (2001), pp. 59-78. 18. Brey J.J., García de Soria M.I., Maynar P., and Ruiz-Montero M.J., Phys. Rev. E, to appear (2004). 820