Heat flux in a vibrated granular gas: the diffusive heat conductivity coefficient

Heat flux in a vibrated granular gas: the diffusive heat conductivity coefficient J. Javier Brey∗and M.J. Ruiz-Montero∗ ∗Física Teórica. Facultad de Física. Apdo. de Correos 1065. 41080-Sevilla. Spain Abstract. The transport coefficient coupling heat flux and density gradient in a granular gas is measured by taking advantage of the existence of a minimum in the temperature profile of an open vibrated granular medium. This temperature inversion is closely related to the existence of the new transport coefficient. Particle simulations using the Direct Simulation Monte Carlo method will be used to compute the transport coefficient, and the results will be compared with theoretical predictions derived from the Boltzmann equation. Finally, the accuracy of a boundary condition requiring the hydrodynamic heat flux to vanish for infinite heights will be discussed. INTRODUCTION Steady states of granular materials are obtained when energy is supplied to the system in order to compensate the energy loss in collisions. One of the easiest ways to add energy to a system in experiments is to do so through a vibrating wall. If the vibration is strong enough, the system will be fluidized, and one can expect its behavior to be described by the extension of the Navier-Stokes equations to granular systems. In this work, we will use the hydrodynamical description to study the steady state of an open, vibrated, granular system in presence of gravity. The granular fluid will be modelled as a system of inelastic hard particles, whose collisions are characterized by a constant coefficient of normal restitution α . Besides, we will be interested in the dilute limit, when the Boltzmann equation applies. A distinctive feature of granular gases as compared to molecular (elastic) ones, is that the expression of the heat flux has to be generalized by including a new term coupling the heat flux and the density gradient. This implies the introduction of a new coefficient, the diffusive heat conductivity µ , that vanishes in the elastic limit. This coupling has been derived by kinetic theory methods [1, 2, 3], and its consequences confirmed in computer simulations [4, 5, 6]. For the particular case of a vibrated granular gas in the presence of gravity, this term implies a peculiar behavior of the temperature profile, that increases with height after a minimum. The existence of the minimum was first derived from the hydrodynamic equations [5], and was confirmed by computer simulations [5, 7, 8], and also in experiments [9]. Here, we will show that the value of the diffusive heat conductivity µ can be obtained from the behavior of the system at the temperature minimum. Then, the value of the transport coefficient will be compared with the theoretical prediction from the Boltzmann equation derived in [2]. Finally, the boundary conditions to be used when solving the hydrodynamic equations will be also discussed. HYDRODYNAMIC DESCRIPTION Let us consider a system of Ninelastic hard spheres (d=3) or disks (d=2) of mass mand diameter σ , in presence of a gravitational field g=−gez, where gis a positive constant and eza unit vector in the Zdirection. In the hydrodynamic description, the state of the system is completely specified by the local number of particles density, n(r,t), the velocity flow u(r,t), and the temperature, T(r,t). For a dilute gas, the evolution of these quantities is given by the extension to the inelastic case of the Navier-Stokes equations [2, 10], ∂ n ∂ t+∇·(nu) = 0,(1) 809 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 Heat ux in a vibrated granular gas: the diffusive heat conductivity coefficient 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) Física Teórica. Facultad de Física. Apdo. de Correos 1065. 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 ∂ u ∂ t+u·∇u+1 nm∇·P−g=0,(2) ∂ T ∂ t+u·∇T+2 dnkB [P:∇u+∇·q]+T ζ =0.(3) Here, kBis the Boltzmann constant, Pis the pressure tensor, P=pI− η (∇u)+(∇u)+−2 dI∇·u,(4) where p=nkBTis the hydrostatic pressure, Ithe identity tensor, and η the shear viscosity coefficient. qappearing in Eq. (3) is the heat flux, q=− κ ∇T− µ ∇n,(5) with κ the heat conductivity coefficient, and µ the diffusive heat conductivity. Finally, ζ is the cooling rate associated to the energy dissipation in collisions. The expression for these coefficients reads: η = η ∗( α ) η 0(T), κ = κ ∗( α ) κ 0(T), µ = µ ∗( α ) µ 0(T),(6) ζ ≃ ζ (0)= ζ ∗( α )p η 0 ,(7) with η 0, κ 0the Boltzmann elastic values of the viscosity and heat conductivity, µ 0=T κ 0/n, while η ∗( α ), κ ∗( α ), µ ∗( α ), and ζ ∗( α )are dimensionless functions of the coefficient of restitution. Although their expressions will not be given here (they can be found in [2, 10]), it is important to remember that in the elastic limit α →1, η ∗and κ ∗tend to unity, while µ ∗and ζ ∗vanish. Besides, both η 0and κ 0are proportional to T1/2. Let us consider that energy is supplied to the system through a vibrating wall of size Slocated at z=0. Although the details of the wall movement are not very important here, we will consider that the vibration amplitude is small enough as to approximate the position of the wall as fixed. Also, the wall moves with a sawtooth profile, so particles find it moving upwards with a characteristic velocity. When the energy supplied by the wall compensates the one lost in collisions, the system reaches a steady state. Besides, and because of the symmetry of the problem, we can expect that, once in this state, there will be gradients only in the Zdirection. When Eqs. (1)–(3) are particularized for this state, they take the form ∂ p ∂ z=−nmg,(8) 2 dnkB dqz dz +T ζ (0)=0.(9) Besides, Eq. (8) implies that the heat flux in this case is qz(z) = −( κ ∗− µ ∗) κ 0 dT dz + µ ∗ κ 0 mg kB .(10) In order to solve these equations, it is convenient to introduce the new dimensionless length scale ξ as: ξ =√aZ∞ z dz01 λ (z0),(11) where λ (z)is the local mean free path for hard disks or spheres, λ (z) = [Cn σ d−1]−1, with C=2√2 for d=2 and C= π √2 for d=3. Let us notice that, when z→∞, ξ →0, while ξ takes its maximum value, ξ 0=√aC σ d−1Nz, with Nz=N/S, at z=0. The coefficient aappearing in Eq. (11) is a function of the coefficient of restitution and reads a( α ) = 32(d−1) π d−1 C2(d+2)3Γ(d/2)2 ζ ∗( α ) κ ∗( α )− µ ∗( α ),(12) vanishing, therefore, in the elastic limit. In terms of the new scale, the solution to the Navier-Stokes equations is [5]: T1/2( ξ ) = ξ − ν [AI ν ( ξ )+BK ν ( ξ )],(13) 810 0 20 40 60 z 0.0 0.5 1.0 1.5 T/TA n/nh FIGURE 1. Temperature and density profiles for a system of hard spheres with α =0.925, ξ 0=1.92. The symbols are from the simulations, while the lines are the solution to the hydrodynamic equations, with the arbitrary constants determined from the temperature minimum. n( ξ ) = mg ξ 1+2 ν CkB σ d−1pa( α )[AI ν ( ξ )+BK ν ( ξ )]2,(14) where I ν and K ν are the modified Bessel functions of first and second kind, and Aand Bare constants to be determined from the boundary conditions. The parameter ν is ν ( α ) = µ ∗( α ) 4[ κ ∗( α )− µ ∗( α )] >0,(15) In the limit ξ →0 (z→∞)I ν →0, while K ν →∞[11]. Nevertheless, this does not imply that the constant Bhas to be taken identically equal to zero, as there is nothing unphysical in the fact that Tdiverges as far as the density decreases fast enough as to guarantee that the local kinetic energy density goes to zero in that limit. Besides, it is clear that the hydrodynamic description will not be valid for very large heights, as the density will have decayed to very small values, and the local Knudsen number will be very small. Then, we will keep the Bconstant different from zero. The presence of the term proportional to K ν in the temperature profile implies that it has a minimum located at ξ = ξ mgiven by AI ν +1( ξ m)−BK ν +1( ξ m) = 0,(16) and the temperature Tmat the minimum is T1/2 m= ξ − ν m[AI ν ( ξ m)+BK ν ( ξ m)].(17) Then, if hydrodynamics is valid in the vicinity of the temperature minimum, the constants Aand Bcan be determined from the measured temperature profiles. It has been also shown [5] that the density has a maximum at ξ = ξ nthat is approximately given by the solution of the equation: I ν ( ξ n)−2 ξ nI ν +1( ξ n) = 0,(18) which can be solved numerically for each value of α . Let us just comment that ξ ntakes values of the order of unity. Of course, in order to observe the density maximum in an experiment the number of particles in the system has to be large enough so that ξ 0> ξ n. If this condition is not fulfilled, the density will decay monotonically with height. In order to check the above hydrodynamic description, computer simulations for two and three dimensional systems by using the Direct Simulation Monte Carlo method (DSMC) [12] have been performed. This method is particularly suited for this problem as we are interested in the low density limit, and as it allows to exploit the symmetry of the system in the transversal direction. In Fig. 1 the steady density and temperature profiles for a system of hard spheres have been plotted. The density has been scaled with the initial, homogeneous density, while the temperature is scaled with some arbitrary value. Space is measured in units of the initial, homogeneous, mean free path. The values of the 811 0.70 0.80 0.90 1.00 α 0.0 0.2 0.4 0.6 µ∗ FIGURE 2. Reduced transport coefficient µ ∗as a function of α for a system of hard spheres. The solid line is the theoretical prediction derived in [2], and the circles the results from the simulation of a vibrated system. The squares are also simulation results but using a Green-Kubo expression for µ . parameters of the simulations were α =0.925, ξ 0=1.92 . The symbols are the results of the simulations: as predicted, the density shows a maximum, while the temperature exhibits a minimum, increasing from there on. The continuous lines are the theoretical prediction, with Aand Bdetermined from the temperature minimum. The agreement is very good, confirming the validity of hydrodynamics in the temperature minimum region. The same qualitative results are achieved for other values of the parameters both in the two and three dimensional cases, as far as the coefficient of restitution is not too low. Due to an intrinsic coupling between gradients and inelasticity in this steady state, retaining only up to the Navier-Stokes order in the Chapman-Enskog expansion may be not enough for small values of α . THE HEAT FLUX In the steady state of a granular material, the heat flux is given by Eq. (11). At the temperature minimum, as the derivative of Tvanishes, we have qz(zm) = µ ∗ κ 0(Tm)mg kB .(19) Then, the computation of the heat flux at the temperature minimum provides a direct measurement of the scaled diffusive heat conductivity coefficient, µ ∗. It must be pointed out that Eq. (19) only requires the validity of the hydrodynamic description, and no additional condition has to be introduced. In Fig. 2 we have plotted the coefficient µ ∗as a function of α for a system of hard spheres. The solid line is the theoretical prediction derived in [2] by using the Chpaman-Enskog expansion from the Boltzmann equation. The circles are the results of the simulation using Eq. (19). The error bars are obtained by giving different values to the parameter ξ 0and the velocity of the vibrating wall. The agreement between theory and simulation is quite good, even at the lowest values of α investigated. Nevertheless, it must be said that, for those values, the shape of the profiles begin to show discrepancies from the theoretical prediction derived here. This is the reason of the large error bars for the lowest values of α , and for not having considered smaller values of this parameter. Finally, we have also included in the figure the results of independent DSMC simulations where µ ∗was computed by means of Green-Kubo expressions derived in Ref. [13] (squares). The agreement is again very good. It could be argued that it is not surprising to obtain a good agreement between the results of the DSMC simulation and a theoretical description derived from the Boltzmann equation. Nevertheless, it must be remembered the the theoretical derivation requires given hypothesis and approximations (validity of hydrodynamics, gradient expansions, Sonine expansion) that are not assumed at all by the simulation method. It is important to stress that the fact that the heat flux does not vanish at the temperature minimum is a direct proof of the existence of the coupling between heat flux and density gradient. Besides, Fig. 2 shows that, although µ ∗→0 in the elastic limit, its contribution cannot be neglected as α goes beyond the quasi-elastic limit. 812 0.70 0.80 0.90 1.00 α 0.0 0.5 1.0 1.5 ξm,n FIGURE 3. Position of the temperature minimum ξ mand of the density maximum ξ nin the three dimensional case. The lines are the theoretical prediction (solid line for ξ m, dot-dashed line for ξ n. The symbols are the results of the DSMC simulations. Eq. (19) provides the heat flux at the temperature minimum. A general expression for qzvalid at any position can be obtained when the solution given by Eq. (14) is substituted in Eq. (11), and it reads qz=A( κ ∗− µ ∗)2mg κ 0 kBT1/2I ν −1( ξ )−B AK ν −1( ξ ) ξ 1− ν .(20) Let us consider now the limit z→∞or, equivalently, ξ →0. Taking into account the asymptotic behavior of the modified Bessel functions, it follows that, for very large heights is given by: qz=A( κ ∗− µ ∗)2mg κ 0 kBT1/22− ν 2 Γ( ν )−Γ(1− ν )B A.(21) If we require that the heat flux vanishes for very large heights, we get the relation B A=2 Γ( ν )Γ(1− ν ),(22) i.e. the ratio of the two constants is given by a function of α alone. Taking into account the behavior of the Γfunction, µ ∗=0 (i.e. ν =0) implies B=0, and the temperature would be constant for large heights. If we introduce the above relation into the equation that determines the position of the temperature minimum, Eq. (16), we get a closed equation for ξ m: I ν +1( ξ m)−2 Γ( ν )Γ(1− ν )K ν +1( ξ m) = 0.(23) Then, the position of the temperature minimum in a vibrated system in the scaled space variable is a function only of the coefficient of normal restitution, independent of the other relevant parameters of the system (number of particles, gravity and vibration velocity). In Fig. 3 the position of the temperature minimum ξ min the three dimensional case is plotted as a function of α . The symbols are the results of the DSMC simulation, while the solid line is the solution of Eq. (23). The agreement is very good, supporting the use of the above mentioned boundary condition to solve the hydrodynamic equations. It must be noticed that the validity of this boundary condition is not at all clear, as it is imposed in the very large heights region, i.e., once the density has decayed to very small values and the hydrodynamic description is not valid. Of course, the heat flux must vanish in that region, but the question is whether this can be translated into an effective boundary condition for the hydrodynamic equations. For instance, a vanishing heat flux implies a constant temperature gradient whose value is directly related to ν . But identifying the asymptotic region where this linear behavior is achieved is 813 very difficult in practical applications, because of the failure of hydrodynamics to describe the upper region of the system [6]. Here we have shown that the condition of vanishing flux at large heights translates into a condition inside the hydrodynamic region, so it can be clearly tested. We have also included in Fig. 3 the position of the density maximum, ξ n. The dot-dashed line is the theoretical prediction given by the solution of Eq. (18), while the triangles are the results of the simulations. The agreement is again quite good, although some discrepancies appear for the lowest values of α studied. The reason for this might be that Eq. (18) is not exact [5], and is in fact obtained for ν 1, while for α =0.7, ν ∼0.14. Nevertheless, it must be noticed the weak dependence on α of ξ nas compared to ξ m. In fact, for the values of the coefficient of restitution considered in the simulations ξ n∼1. In conclusion, hydrodynamics provides a very useful tool to study non-homogeneous steady states of granular systems. Nevertheless, this hydrodynamic description has distinctive characteristics that cannot be guessed from the one of molecular fluids. In particular, a new transport coefficient, the diffusive heat conductivity, has to be introduced. In this work we have shown that this new transport coefficient has relevant consequences in the hydrodynamic profiles, so it cannot be neglected in the description of granular flows. ACKNOWLEDGMENTS We acknowledge financial support from the Ministerio de Ciencia y Tecnología (Spain) through Grant No. BFM200200303 (partially financed by FEDER funds). REFERENCES 1. Brey, J. J., Moreno, F., and Dufty, J. W., Phys. Rev. E 54, 445–456 (1996) 2. Brey,J. J., Dufty, J. W., Kim, C. S., and Santos, A., Phys. Rev. E 58, 4638–4653 (1998). 3. Sela, N., and Goldhirsch, I., J. Fluid Mech., 361, 41–74 (1998). 4. Soto, R., Mareschal, M., and Risso, D., Phys. Rev. Lett. 83, 5003–5006 (1999). 5. Brey, J.J., Ruiz-Montero, M.J., and Moreno, F., Phys. Rev. E 63, 061305 (2001). 6. Brey, J. J., Ruiz-Montero, M. J., Europhys. Lett. 66, 805–811 (2004). 7. Helal, K., Biben, T., and Hansen, J. P., Physica A 240, 361–373 (1997). 8. Ramírez, R., and Soto, R., Physica A 322, 73–80 (2003). 9. Blair, D. L., and Kudrolli, A., Phys. Rev. E 67, 061311 (2001). 10. Brey, J. J., Cubero, D., in Granular Gases, Pöschel, T., and Luding, S. eds., Lectures Notes in Physics, Springer Verlag (Berlín), 59–78 (2001). 11. Handbook of mathematical functions, Abramowitz, M., and Stegun, I. A., Dover, (New York, 1965). 12. Bird, G., Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Clarendon Press (Oxford, 1994). 13. Dufty, J. W., and Brey, J. J., Phys. Rev. E 68, 030302 (R) (2003). 814