Hydrodynamics of an open vibrated granular system

Hydrodynamics of an open vibrated granular system J. Javier Brey, M. J. Ruiz-Montero, and F. Moreno Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, 41080 Sevilla, Spain 共Received 17 January 2001; published 21 May 2001兲 Using the hydrodynamic description and molecular dynamics simulations, the steady state of a fluidized granular system in the presence of gravity is studied. For an open system, the density profile exhibits a maximum, while the temperature profile goes through a minimum at high altitude, beyond that the temperature increases with the height. The existence of the minimum is explained by the hydrodynamic equations if the presence of a collisionless boundary layer is taken into account. The energy dissipated by interparticle collisions is also computed. A good agreement is found between theory and simulation. The relationship with previous works is discussed. DOI: 10.1103/PhysRevE.63.061305 PACS number共s兲: 45.70.Mg, 81.05.Rm, 51.10.⫹y, 05.20.Dd I. INTRODUCTION The aim of this paper is to investigate the hydrodynamic description of a granular fluid in a gravitational field when energy is continuously provided to the system from below through a vibrating plate. Experiments, computer simulations, and also theoretical studies have revealed the existence of a steady fluidized state under these conditions 关1–8兴. One of the main conclusions of all these studies is that vibrofluidized granular media show a fluidlike behavior, in the sense that their state seems to be well characterized by the profiles of the hydrodynamic fields. Of course, a different question is whether hydrodynamic equations, derived as an extension of those for ordinary fluids, provide a quantitatively, or at least qualitatively, accurate description of what is observed. In a previous work 关9兴, we have analyzed a vibrated granular fluid in absence of external fields, paying special attention to the bulk behavior of the fluid, far away from the boundaries. The system was shown to exhibit a normal behavior, independent of the details of the boundaries, characterized by a closed constitutive relationship between the uniform pressure and the temperature gradient. The idea here is to extend the above study to systems submitted to a uniform external field. Such an extension is not at all trivial, since the field generates spatial gradients that couple in an intricate way to those associated with the inelasticity of the system. Most of the experimental studies deal with open systems, i.e., formally with a system of infinite height, and that is the situation we will focus on here. It could be expected that the behavior of the system becomes simple far away from the vibrating surface at the bottom. Nevertheless, the presence of a free surface introduces some additional complications in the description of the fluid. As the density decreases, particles tend to move in a ballistic way restrained by the gravitational force. This was already noted by Haff 关1兴, who concluded that a hydrodynamiclike description cannot account for those effects. In 1991, Clement and Rajchenbach 关2兴carried out an interesting experimental study of a fluidized two-dimensional vertical granular system. They measured the hydrodynamic profiles and established a series of important observations. The point we want to emphasize is presented in Fig. 4 of their paper, where it is observed that the granular temperature increases for large enough heights. Unfortunately, the authors do not comment on the origin and relevance of this finding. The same effect has been found again in experiments and in molecular dynamics simulations by Helal et al. 关10兴. The temperature profile presents a minimum at a high altitude, beyond that it is found to increase. In order to explain this behavior, the authors suggest a model that involves two coupled differential equations for the packing fraction and the temperature. These equations must be numerically solved by using boundary conditions determined from the simulations. Although there is a good qualitative agreement with the experimental and simulation results, the physical origin of the rise in temperature as well as its compatibility with a hydrodynamic description of the system is not clear. In particular, the continuous approach used by Helal et al. does not include a density dependent heat flux term, that plays a relevant role in granular systems 关11,12兴. A detailed discussion of these questions will be one of the main points to be addressed in the following sections. The most direct implication of inelasticity in collisions is the dissipation of energy in the system. The balance between the energy dissipated and the energy supplied through the vibrating wall is often used to derive scaling laws for vibrofluidized granular materials, as well as boundary conditions to solve the hydrodynamic equations. Warr et al. 关13兴modeled the vibrated granular system under gravity as an isothermal fluid, with all the particles having the same velocity. Later on, Kumaran 关14兴, using a kinetic theory description, found that the velocity distribution function is a MaxwellBoltzmann distribution in the limit of small inelasticity. His expression for the dissipated power Donly differs from the result derived by Warr et al. 关13兴by a constant. Extensive molecular dynamics simulations carried out by McNamara and Luding 关5兴showed that the prediction in Ref. 关14兴fitted better the simulation data than the expression in Ref. 关13兴, although significant discrepancies were found. In particular, the simulation values of Dindicated a dependence on the number of particles in the system that was not accounted for by the theory. Trying to understand this dependence is one of the goals of the present work. This paper is organized as follows. In Sec. II, the NavierStokes-like hydrodynamic equations for a granular gas are PHYSICAL REVIEW E, VOLUME 63, 061305 1063-651X/2001/63共6兲/061305共10兲/$20.00 ©2001 The American Physical Society63 061305-1 shortly reviewed and particularized for the steady state of a vibrated system under the influence of an homogeneous external force. By introducing an appropriate scaling for the space coordinate in the direction of the field, explicit expressions for the hydrodynamic profiles are obtained. They involve two constants that must be determined from the boundary conditions. The limit of an infinitely high open system is considered in Sec. III. This limit does not imply that any of the two constants appearing in the general expression of the temperature profile must vanish, contrary to what has been established in some previous works. This is due to the presence of a free-molecule boundary layer in the upper region of the granular gas, so that the validity of the hydrodynamic equations cannot be extrapolated up to an infinite height, for which a divergent behavior of one of the contributions to the theoretical prediction for the temperature profile would show up.A relevant consequence of the above is that the temperature profile exhibits a minimum, becoming an increasing function for large enough heights. Moreover, by using molecular dynamics simulations, we have verified that the region of increasing temperature is accurately described by the hydrodynamic equations. In fact, the values of the minimum of the temperature and its position provide enough information to build up the hydrodynamic profiles in the bulk of the system, as discussed in the last part of Sec. III. The form of the profiles for large heights, but still in the hydrodynamic regime, is analyzed in Sec. IV. Section V is devoted to the study of the dissipation in collisions, again by means of the hydrodynamic description. The results are compared with previous works as well as with molecular dynamics simulations. The origin of the discrepancies between the several theories is clarified, showing, for instance, that the expression for Dderived by Kumaran 关14兴corresponds to the low inelasticity and small system limits of the more general form derived here. Finally, the main conclusions are summarized in Sec. VI. II. HYDRODYNAMIC EQUATIONS We consider a granular gas composed by smooth inelastic hard spheres (d⫽3) or disks (d⫽2) of mass mand diameter ␴ , in presence of a uniform external force f. The particles collide with a constant coefficient of normal restitution ␣ .In the hydrodynamic description, it is assumed that the state of the system is characterized by the local number density n(r,t), velocity flow u(r,t), and temperature T(r,t)关1,15兴. For a dilute gas, the time evolution of these fields is given by the equations 关16,17兴 ⳵ tn⫹“•共nu兲⫽0, 共1兲 ⳵ tui⫹u•“ui⫹共mn兲⫺1ⵜip⫺共nm兲⫺1ⵜj 共2兲 ⫻ 冋 ␩ 冉 ⵜiuj⫹ⵜjui⫺2 d ␦ ij“•u 冊 册 ⫺m⫺1fi⫽0, ⳵ tT⫹u•“T⫹2共dnkB兲⫺1p“•u⫺2共dnkB兲⫺1ⵜiuj ⫻ 冋 ␩ 冉 ⵜiuj⫹ⵜjui⫺2 d ␦ ij“•u 冊 册 ⫺2共dnkB兲⫺1“•共 ␬ “T⫹ ␮ “n兲⫹T ␨ (0)⫽0. 共3兲 In the above equations, p⫽nkBTis the pressure, kBthe Boltzmann constant, ␩ the shear viscosity coefficient, ␬ the heat conductivity coefficient, and ␮ a new transport coefficient that has no analogous in the elastic limit. Although kBis taken as unity in most of the literature of fluidized granular systems, we will keep it here just to stress the analogy with molecular fluids. Finally, ␨ (0) is the cooling rate associated to the energy dissipation in collisions. These quantities have the form ␩ ⫽ ␩ *共 ␣ 兲 ␩ 0共T兲, ␬ ⫽ ␬ *共 ␣ 兲 ␬ 0共T兲, ␮ ⫽ ␮ *共 ␣ 兲 ␮ 0共T兲, 共4兲 ␨ (0)⫽ ␨ *共 ␣ 兲p ␩ 0,共5兲 where ␩ 0and ␬ 0are the Boltzmann elastic values of the shear viscosity and heat conductivity, respectively, ␮ 0 ⫽T ␬ 0/n, and ␬ *, ␩ *, ␮ *, and ␨ *are dimensionless functions of the coefficient of normal restitution ␣ . For ␣ →1, ␩ *and ␬ *tend to unity, while ␮ *and ␨ *vanish. The explicit expressions of these quantities are given in Appendix A. Let us note that the existence of a transport coefficient ␮ giving a contribution of the density gradients to the heat flux is a peculiarity of granular fluids that has been confirmed by molecular dynamics 关11兴and by Monte Carlo simulations of the Boltzmann equation 关12兴. We will consider a force of the gravitational type, namely, f⫽⫺mge ˆz,共6兲 with ga positive constant and e ˆzthe unit vector in the positive direction of the zaxis. The system we will study is a granular medium of Nparticles contained in a box of section Sand height L. The quantity Sis an area for d⫽3 and a length for d⫽2. Energy is added to the system by the bottom of the box which vibrates in a given way. Since many of the results we will derive in the following are independent of the specific way in which the wall is vibrated, we delay the discussion of the details of its motion until they are needed. Moreover, we are not interested in the boundary effects associated to the side walls and, therefore, we will consider a region of the system far away from them. In fact, in the molecular dynamics simulations to be reported later on, periodic boundary conditions were used in the directions perpendicular to the external field. Then, because of symmetry considerations, in the steady state only gradients in the z direction are expected and the hydrodynamic equations reduce to ⳵ p ⳵ z⫽⫺nmg,共7兲 J. JAVIER BREY, M. J. RUIZ-MONTERO, AND F. MORENO PHYSICAL REVIEW E 63 061305 061305-2 2 dnkB ⳵ ⳵ z 冉 ␬ ⳵ T ⳵ z⫹ ␮ ⳵ n ⳵ z 冊 ⫺T ␨ (0)⫽0. 共8兲 Equation 共7兲implies T n ⳵ n ⳵ z⫽⫺ mg kB ⫺ ⳵ T ⳵ z,共9兲 and use of this expression into Eq. 共8兲together with Eqs. 共4兲 and 共5兲, leads to, 2 dnkB关 ␬ *共 ␣ 兲⫺ ␮ *共 ␣ 兲兴 ⳵ ⳵ z 冋 ␬ 0共T兲 ⳵ T ⳵ z 册 ⫺2mg dnkB 2 ␮ *共 ␣ 兲 ⳵ ␬ 0共T兲 ⳵ z⫺ ␨ *共 ␣ 兲nkBT2 ␩ 0 ⫽0. 共10兲 In order to analyze the above equation, it is convenient to introduce a dimensionless length scale lby l⫽ 冕 z Ldz⬘1 ␭共z⬘兲,共11兲 where ␭(z) is the local mean free path for hard disks or spheres, ␭共z兲⫽关Cn ␴ d⫺1兴⫺1,共12兲 with C⫽2 冑 2 for d⫽2 and C⫽ ␲ 冑 2 for d⫽3. The variable lmeasures the number of mean free paths from the wall located at z⫽Lto the parallel plane located at height z. For z⫽0itis l共z⫽0兲⬅l0⫽C ␴ d⫺1Nz,共13兲 Nz⫽N/Sbeing the number of particles in the system per unit of section. In terms of l, the solution of Eq. 共7兲can be written as p⫽mgl C ␴ d⫺1⫹pL,共14兲 where pLis the pressure of the gas next to the upper wall. In particular, at z⫽0, p共z⫽0兲⬅p0⫽mgNz⫹pL.共15兲 Equation 共10兲is equivalent to ⳵ 2T1/2 ⳵ l2⫹b共 ␣ 兲 l⫹pL * ⳵ T1/2 ⳵ l⫺a共 ␣ 兲T1/2⫽0, 共16兲 with pL *⫽C ␴ d⫺1 mg pL共17兲 and a共 ␣ 兲⫽32共d⫺1兲 ␲ d⫺1 C2共d⫹2兲3⌫共d/2兲2 ␨ *共 ␣ 兲 ␬ *共 ␣ 兲⫺ ␮ *共 ␣ 兲,共18兲 b共 ␣ 兲⫽2 ␬ *共 ␣ 兲⫺ ␮ *共 ␣ 兲 2关 ␬ *共 ␣ 兲⫺ ␮ *共 ␣ 兲兴.共19兲 Inspection of Eq. 共16兲indicates that 冑 a( ␣ ) determines the coupling between gradients and dissipation that is intrinsic in granular flows, while b( ␣ ) is a scaling factor of the inhomogeneities associated to the external field. It is still possible to express Eq. 共16兲in a more familiar way by defining a new variable ␰ by ␰ ⫽ 冑 a共 ␣ 兲共l⫹pL *兲⫽C ␴ d⫺1 冑 a共 ␣ 兲 冋 冕 z Ldz⬘n共z⬘兲⫹pL mg 册 . 共20兲 This variable, as well as l, is a decreasing function of the original coordinate z. Note that ␰ cannot be defined in the elastic limit ␣ →1, in which a( ␣ ) vanishes. Performing the change, Eq. 共16兲becomes ␰⳵ 2T1/2 ⳵␰ 2⫹b共 ␣ 兲 ⳵ T1/2 ⳵␰ ⫺ ␰ T1/2⫽0, 共21兲 whose general solution is T1/2共 ␰ 兲⫽A ␰ ⫺ ␯ I ␯ 共 ␰ 兲⫹B ␰ ⫺ ␯ K ␯ 共 ␰ 兲,共22兲 where ␯ 共 ␣ 兲⫽ ␮ *共 ␣ 兲 4关 ␬ *共 ␣ 兲⫺ ␮ *共 ␣ 兲兴 ⬎0, 共23兲 I ␯ and K ␯ are the modified Bessel functions of first and second kind, respectively, and Aand Bare constants that must be determined from the boundary conditions. Since the behavior of the functions 冑 a( ␣ ) and ␯ ( ␣ ) will play an important role in the discussions in the following sections, we have plotted them in Fig. 1 for d⫽2. Although both quantities vanish in the elastic limit, 冑 a( ␣ ) grows much faster than ␯ ( ␣ )as ␣ decreases in the vicinity of ␣ ⫽1. On the other hand, when ␣ becomes smaller, the behavior inverts and ␯ ( ␣ ) presents a much larger slope. The pressure profile in the ␰ scale follows directly from Eq. 共14兲and 共20兲, p共 ␰ 兲⫽mg ␰ C ␴ d⫺1 冑 a共 ␣ 兲.共24兲 In this way, we have formally solved the hydrodynamic equations for the system under consideration. The expression for the density follows from Eqs. 共22兲,共24兲, and the equation of state. Afterwards, the relationship between ␰ and the original coordinate zis obtained by solving the equation HYDRODYNAMICS OF AN OPEN VIBRATED GRANULAR SYSTEM PHYSICAL REVIEW E 63 061305 061305-3 d ␰ n共 ␰ 兲⫽⫺ 冑 a共 ␣ 兲C ␴ d⫺1dz,共25兲 that is the differential form of the definition of ␰ given in Eq. 共20兲. Of course, the solution of this equation involves the pressure of the gas next to the upper wall pL, which is unknown up to now. III. OPEN SYSTEMS In order to particularize the general results obtained in the previous section for a given physical situation, i.e., for specific forms of the walls at z⫽0 and z⫽L, the crucial and nontrivial point is the introduction of the hydrodynamic boundary conditions needed to determine the constants Aand Bappearing in Eq. 共22兲, as well as the value of pLentering in the definition of ␰ , Eq. 共20兲. We are interested in an open system, i.e., in the limit of infinite height L. In this limit, it is obvious that pL⫽0, so that the lowest value of ␰ is now ␰ ⫽0, corresponding to the limit z→⬁. More explicitly, Eq. 共20兲becomes ␰ ⫽C ␴ d⫺1 冑 a共 ␣ 兲 冕 z ⬁dz⬘n共z⬘兲.共26兲 The behavior of the modified Bessel functions in the limit ␰ →0is关18兴 I ␯ 共 ␰ 兲⬃1 ⌫共1⫹ ␯ 兲 冉 ␰ 2 冊 ␯ ,共27兲 K ␯ 共 ␰ 兲⬃⌫共 ␯ 兲 2 冉 ␰ 2 冊 ⫺ ␯ .共28兲 Therefore, it follows from Eq. 共22兲that in the same limit, T1/2共 ␰ 兲⬃B⌫共 ␯ 兲 2 冉 ␰ 2 2 冊 ⫺ ␯ ,共29兲 indicating a divergent behavior of the temperature as the limit ␰ →0(z→⬁) is approached. Nevertheless, it cannot be concluded from here that the constant Bmust identically vanish, contrary to what has been inferred in other works 关1,4兴. There are two main reasons for that. First, the fact that Tformally diverges does not imply anything unphysical, as long as the density decreases fast enough as to guarantee that the local kinetic energy goes to zero as zgoes to infinity. Second, Eq. 共22兲is based on a continuous hydrodynamic description of the granular flow, and such a description is not valid in the region in which ␰ is very small, so that the local Knudsen number, defined as the ratio of the mean free path to the length scale of the macroscopic gradients, is very large. There, the arguments leading from a microscopic description to a continuous approach fail, and the gas has to be described as a free-molecule flow, with a so-called transition regime between the hydrodynamic region and the collisionless one 关19,20兴. The detailed analysis of the gas in these regimes is an interesting but very complex problem that will be addressed elsewhere. So, we will keep the constant Bin Eq. 共22兲different from zero, although the question is still whether its contribution is relevant within the hydrodynamic region. The presence of the term proportional to K ␯ ( ␰ ) in the expression of T1/2 implies that the temperature profile exhibits a minimum. Using the expressions of the derivatives of the modified Bessel functions 关18兴, it is obtained that the minimum is located at a value ␰ ⫽ ␰ Tgiven by the solution of the equation AI ␯ ⫹1共 ␰ T兲⫺BK ␯ ⫹1共 ␰ T兲⫽0, 共30兲 and the temperature Tmat the minimum is Tm 1/2⫽ ␰ T ⫺ ␯ 关AI ␯ 共 ␰ T兲⫹BK ␯ 共 ␰ T兲兴.共31兲 Therefore, if the hydrodynamic description is valid in the vicinity of ␰ ⫽ ␰ T, the values of ␰ Tand Tmallow us to determine the constants Aand Bcharacterizing the hydrodynamic profiles everywhere in the system, except in the boundary layers next to ␰ ⫽ ␰ 0and ␰ ⫽0, where a more microscopic description, such as that provided by kinetic theory, is needed. As an example, in Fig. 2 we plot the temperature profile in the ␰ variable obtained by molecular dynamics simulation in a two-dimensional system with ␣ ⫽0.95 and Nz⫽6. The wall at the bottom is vibrated with a sawtooth velocity profile having a velocity vW⫽6. This means that all the particles colliding with the wall find it with that velocity 关5,21兴. Moreover, the amplitude of the wall motion is much smaller than the mean free path of the particles next to it, so that the position of the wall can be taken as fixed at z⫽0. Periodic boundary conditions are employed in the direction perpendicular to the field. The units are defined by m⫽1, and ␴ ⫽1. We take kB⫽0.5, and the value of the external field is g⫽1. For the above value of ␣ ,itis ␯ ⯝0.021. From the simulation data it is estimated that ␰ T⯝0.21 and Tm⯝162.6, and using Eqs. 共30兲and 共31兲one gets A⯝12.2 and B ⯝0.76. The dashed line in the figure is the temperature proFIG. 1. Functions 冑 a( ␣ )共solid line兲and ␯ ( ␣ )共dashed line兲 defined in the main text, for d⫽2. J. JAVIER BREY, M. J. RUIZ-MONTERO, AND F. MORENO PHYSICAL REVIEW E 63 061305 061305-4 file obtained with these values of the constants. A quite good agreement is observed between the constructed profile and the simulation data outside the boundary layers. In particular, the agreement extends well inside the region of small values of ␰ . This confirms that the hydrodynamic regime includes the minimum of the temperature and also a part of the system where the temperature increases with z, i.e., Tincreases as ␰ decreases. In Fig. 3 the same temperature profile as in Fig. 2 is shown as a function of the original variable z. It is seen that the increase of the temperature with the height is not just a theoretical artifact, but in practice it is observed over a wide region in real space. Similar results have been obtained for other values of ␣ in the interval 0.85⭐ ␣ ⭐0.99. Details of the practical limitations in the molecular dynamics simulations will be given in the last section of the paper. Since the existence of a region where the temperature increases is associated with the presence of the collisionless boundary layer and the transition regime, the location of the temperature minimum and, consequently, of the hydrodynamic region beyond it, are expected to correspond to small values of ␰ . Therefore, to study this region we can approximate in Eq. 共22兲the modified Bessel functions by their expressions in the limit of small arguments. This yields A B⬃21⫺2 ␯ 共 ␯ ⫹1兲⌫共 ␯ ⫹1兲2 ␰ T ⫺2(1⫹ ␯ )共32兲 and, since ␰ Tis small, it follows that A/BⰇ1. The conclusion reached in this way is that the term involving K ␯ ( ␰ )in Eq. 共22兲is only relevant in the region in which ␰ is small and K ␯ ( ␰ ) is large. More precisely, a detailed asymptotic analysis of Eq. 共22兲indicates that the relevant increasing temperature region corresponds to ␰ 2Ⰶ1 but ␰ 2ⲏ ␯ 2. This range of values of ␰ exists as long as ␯ is small enough, i.e., the system be not too inelastic. Moreover, the analysis shows that in this ␰ window one can approximate K ␯ 共 ␰ 兲⬃⫺ln ␰ .共33兲 When ␰ takes values of the order of unity, the term B ␰ ⫺ ␯ ln ␰ is negligible as compared with A ␰ ⫺ ␯ I ␯ ( ␰ ). Then, we propose, as an accurate approximation to Eq. 共22兲, the expression T1/2共 ␰ 兲⯝A ␰ ⫺ ␯ I ␯ 共 ␰ 兲⫺B ␰ ⫺ ␯ ln ␰ .共34兲 In Fig. 4 we compare Eqs. 共22兲and 共34兲for A⫽12.2 and B⫽0.26, which are the values of the constants found from the molecular dynamics data corresponding to the situation described in Fig. 2. The agreement is fairly good in the plotted interval 10⫺2⬍ ␰ ⬍2. The accuracy is even better for larger values of ␰ , as expected from the above discussion. Let us next analyze the density profile. Equations 共24兲and 共34兲yield FIG. 2. Temperature profile in units defined in the main text, in the ␰ variable of a vibrated system with ␣ ⫽0.95, Nz⫽6( ␰ 0 ⫽1.722) and vw⫽6. The solid line is from the molecular dynamics simulation, while the dashed line is the theoretical prediction discussed in the text. FIG. 3. The same as Fig. 2, but in terms of the real space variable z, measured in units of ␴ . FIG. 4. Comparison of the exact theoretical expression for the temperature profile 关Eq. 共22兲兴 and the approximated expression 关Eq. 共34兲兴 for values of the parameters corresponding to the situation of Fig. 2. HYDRODYNAMICS OF AN OPEN VIBRATED GRANULAR SYSTEM PHYSICAL REVIEW E 63 061305 061305-5 n共 ␰ 兲⫽p共 ␰ 兲 kBT共 ␰ 兲⫽mg ␰ 1⫹2 ␯ CkB ␴ d⫺1 冑 a共 ␣ 兲关AI ␯ 共 ␰ 兲⫺Bln ␰ 兴2. 共35兲 The comparison of this expression with the simulation data for the same system as in Fig. 2 is shown in Fig. 5. The dashed line, corresponding to the theoretical prediction, has been plotted by using the values of Aand Bobtained by fitting the temperature minimum. Again, the agreement is quite good, outside the boundary layer next to the vibrating wall. The density profile exhibits a maximum at ␰ ⫽ ␰ n, that is approximately given by the solution of the equation I ␯ 共 ␰ n兲⫺2 ␰ nI ␯ ⫹1共 ␰ n兲⫽0, 共36兲 which can be numerically solved for each value of ␯ , i.e., of ␣ . A sketch of the derivation of Eq. 共36兲is given in Appendix B. In Fig. 6, ␰ nis shown as a function of ␣ for 0.5⭐ ␣ ⭐0.99. Let us remark that contrary to the position of the temperature minimum ␰ T, ␰ ntakes values of the order of unity. Moreover, the dependence of ␰ non ␣ is rather weak since it grows roughly from 1.07 to 1.46 as ␣ decreases from 0.99 to 0.5. Nevertheless, when the values of ␰ nare translated into the lscale by means of Eq. 共20兲with pL *⫽0, the position of the maximum of the density turns out to be strongly influenced by the inelasticity of the system, increasing very fast as ␣ approaches unity. One of the main implications of Eq. 共36兲is that the position of the density maximum, measured in the scale ␰ , does not depend on the boundary conditions, i.e., on the values of the constants Aand B, or on the total number of particles in the system, as measured for instance by Nz. Of course, in order to actually observe the maximum in an experiment or in a computer simulation, it must be ␰ n⬍ ␰ 0⬅ 冑 a共 ␣ 兲l0⫽C 冑 a共 ␣ 兲 ␴ d⫺1Nz.共37兲 Taking into account Eq. 共11兲, the total number of particles Nz (⫹)per unit of length or area of the vibrating wall above the position of the density maximum is Nz (⫹)⫽ln C ␴ d⫺1共38兲 where lndenotes the position of the density maximum in the lscale. This number increases as ␣ increases. Therefore, the more elastic the system is the larger is the number of particles needed in order to see the maximum of the density. This explains why in some molecular dynamics simulations the density profile shows an almost monotonic decay with an apparent maximum next to the vibrating wall. On the other hand, the position of the density maximum in the actual variable zdoes depend on the boundary conditions, since the conversion from ␰ to zinvolves the constants Aand B,asit follows from Eqs. 共26兲and 共35兲. This is clearly seen, for instance, in the experimental results shown in Fig. 2共b兲of Ref. 关13兴. The prediction of our theory is that the area bellow the density profile to the right of the maximum is the same for the several vibration amplitudes. This seems to be qualiFIG. 5. Density profile in the scaled ␰ variable 共a兲, and in the real variable 共b兲, for the same values of the parameters as in Fig. 2. The solid line is from the molecular dynamics simulation, while the dashed line is the theoretical prediction discussed in the text. FIG. 6. Theoretical value of the scaled position of the density maximum as a function of the coefficient of restitution, ␣ . J. JAVIER BREY, M. J. RUIZ-MONTERO, AND F. MORENO PHYSICAL REVIEW E 63 061305 061305-6 tatively true in the reported case, although the vibration amplitude can affect the degree of fluidization of the granular system, so that for small amplitudes the theory developed here may not apply. IV. THE UPPER REGION OF THE GRANULAR GAS For large z共small ␰ ), but still inside the region where hydrodynamics holds, the temperature and density profiles can be approximated by T1/2共 ␰ 兲⬃A 冉 1 2 冊 ␯ 1 ⌫共1⫹ ␯ 兲⫺B ␰ ⫺ ␯ ln ␰ ,共39兲 n共 ␰ 兲⬃mg ␰ CkB ␴ d⫺1 冑 a共 ␣ 兲关A2⫺ ␯ ⌫共1⫹ ␯ 兲⫺1⫺B ␰ ⫺ ␯ ln ␰ 兴2, 共40兲 respectively. Note that the own structure of Eq. 共39兲, that does not present any minimum, implies that this approximation only holds for values of ␰ smaller than the position ␰ Tof the temperature minimum. On the other hand, ␰ should be large enough as to the hydrodynamic description be accurate. In the previous section we have discussed the existence of a relevant region verifying both conditions. Substitution of these expressions into Eq. 共25兲yields d ␰ ⬃⫺mg ␰ dz kB 冋 A 冉 1 2 冊 ␯ ⌫共1⫹ ␯ 兲⫺1⫺B ␰ ⫺ ␯ ln ␰ 册 2共41兲 and by differentiation of Eq. 共39兲one gets T1/2 dT dz ⬃2mgB kB ␰ ⫺ ␯ 共1⫺ ␯ ln ␰ 兲.共42兲 For ␯ Ⰶ1, that means not very inelastic systems, the above equation can be approximated by dT3/2 dz ⬃3mgB kB.共43兲 This is compatible with what is seen in the molecular dynamics simulations, although due to the relatively small variation of the temperature with zin the region with positive slope, the behavior predicted by Eq. 共43兲is hard to discern from a simple linear in zprofile. We have fitted the temperature profiles obtained by molecular dynamics simulations for large zto the behavior predicted by Eq. 共43兲. The values of B obtained in this way were compared with those determined from the minimum of the temperature, as discussed in Sec. III, and a good agreement was found. Since the constant Bis small, a good estimation of the temperature profile in this upper region of the vibrated granular system is obtained in some cases by considering that the temperature reaches a constant plateau 关13,14,22兴. Also for ␯ Ⰶ1, Eq. 共40兲leads to dlnn dz ⬃⫺mg kBT.共44兲 Therefore, if the temperature is approximately constant in the zinterval considered, an apparently exponential behavior of the density can be observed. This is equivalent to saying that n( ␰ )⬀ ␰ , as easily seen from Eq. 共40兲. Again, the simulations have confirmed these predictions. It is important to stress that the region in which the temperature shows a minimum followed by an apparently linear profile, and the density seems to decay exponentially, can only be explained correctly if the contribution of the term B ␰ ⫺ ␯ K ␯ ( ␰ ) to the temperature profile in Eq. 共22兲is taken into account. Moreover, the above discussion only explains an approximately exponential decay of the density for values of zlying well inside the region of increasing temperature, but not where the temperature decreases with the height, due to the different variation rates of the temperature in both regions. V. DISSIPATED POWER An expression for the total power Ddissipated in the system is directly obtained from the hydrodynamic equation for the temperature, Eq. 共3兲, D⫽ 冕 drdnkB 2T ␨ (0) ⫽dkBS ␨ *共 ␣ 兲T1/2 2C ␴ d⫺1 冑 a共 ␣ 兲 ␩ 0 冕 0 ␰ 0d ␰ p共 ␰ 兲T共 ␰ 兲1/2.共45兲 Upon writing the above expression we have taken into account that ␩ 0is proportional to T1/2. Using now the expressions for the pressure and temperature profiles, Eqs. 共22兲and 共24兲, respectively, one gets D⫽dkBSmg ␨ *共 ␣ 兲T1/2 2关C ␴ d⫺1 冑 a共 ␣ 兲兴2 ␩ 0 冕 0 ␰ 0d ␰␰ 1⫺ ␯ 关AI ␯ 共 ␰ 兲⫹BK ␯ 共 ␰ 兲兴. 共46兲 The integral on the right hand side of this expression can be easily evaluated by employing Eqs. 共B2兲with the result D⫽dkBS ␨ *共 ␣ 兲mgT1/2 2关C ␴ d⫺1 冑 a共 ␣ 兲兴2 ␩ 0 再 A 冋 ␰ 0 1⫺ ␯ I ␯ ⫺1共 ␰ 0兲⫺21⫺ ␯ ⌫共 ␯ 兲 册 ⫺B关 ␰ 0 1⫺ ␯ K1⫺ ␯ 共 ␰ 0兲⫺2⫺ ␯ ⌫共1⫺ ␯ 兲兴 冎 .共47兲 This expression is an exact consequence of the hydrodynamic equations for a granular gas. In particular, no assumption has been made about the values of ␰ 0,A,orB. HYDRODYNAMICS OF AN OPEN VIBRATED GRANULAR SYSTEM PHYSICAL REVIEW E 63 061305 061305-7 Let us define a dimensionless quantity Fby F⫽D g共NmE ¯ ␰ 0兲1/2 ,共48兲 where E ¯ is the total kinetic energy of the system, E ¯ ⫽ 冕 drd 2nkBT⫽SdkB 2C ␴ d⫺1 冑 a共 ␣ 兲 冕 0 ␰ 0d ␰␰ ⫺2 ␯ 关AI ␯ 共 ␰ 兲 ⫹BK ␯ 共 ␰ 兲兴2.共49兲 Substitution of Eqs. 共47兲and 共49兲into Eq. 共48兲leads to an explicit expression for F. It is a rather complicated and not very illuminating expression. Therefore, we will not write it here explicitly, although it will be referred to in the following as Fexact . Let us now suppose that we neglect the part of the profiles that are responsible for the increase of the temperature, i.e., we formally take B⫽0. Because of the discussion in Sec. III, this could be expected a priori to be a good approximation as long as ␰ 0is not small. Then, using the explicit expression of the elastic shear viscosity ␩ 0,itis obtained Fapprox⫽4共2d兲1/2 ␲ d⫺1/2 ␨ *共 ␣ 兲 共d⫹2兲⌫共d/2兲C 冑 a共 ␣ 兲 ␰ 0 ␰ 0 1⫺ ␯ I ␯ ⫺1共 ␰ 0兲⫺21⫺ ␯ ⌫共 ␯ 兲 冋 冕 0 ␰ 0d ␰␰ ⫺2 ␯ I ␯ 共 ␰ 兲2 册 1/2 . 共50兲 Although this expression is not at all simple, it only depends on the values of ␣ and ␰ 0, but not on the boundary conditions that determine the constants Aand B, then representing a scaling law prediction. In the limit of a large system, in the sense that ␰ 0Ⰷ1, the asymptotic behavior of Eq. 共51兲is given by Fapprox⬃8d1/2 ␲ (d⫺1)/2 ␨ *共 ␣ 兲 共d⫹2兲C⌫共d/2兲 冑 a共 ␣ 兲,共51兲 while for small ␰ 0it is Fapprox⬃2共2d兲1/2 ␲ (d⫺1)/2 ␨ *共 ␣ 兲 共d⫹2兲⌫共d/2兲C 冑 a共 ␣ 兲 ␰ 0 1/2.共52兲 Kumaran 关5,14兴modeled the vibrated granular media as an isothermal fluid with a Maxwellian velocity distribution and derived an expression for the dissipated power Din a two-dimensional system. His expression leads to a value of the quantity Fgiven by F⫽共1⫺ ␣ 兲 冋 ␲␴ Nz 2 冑 2a共 ␣ 兲 册 1/2 .共53兲 In the limit of quasielasticity, i.e., for ␣ very close to unity, this result is equivalent to Eq. 共52兲. Consequently, as derived here, its applicability is restricted to small inelasticity and small systems. In Fig. 7 we have plotted the function F( ␰ ) in the interval 0⭐ ␰ 0⭐10 for ␣ ⫽0.95 ( ␯ ⯝0.021). The solid line is Fapprox as given by Eq. 共50兲, while the dashed line is the expression derived by Kumaran, Eq. 共53兲. The dotted line shows the asymptotic constant value predicted by Eq. 共51兲,F⯝0.85. Similar behaviors are obtained for other values of ␣ . The symbols are molecular dynamics simulation results. While the circles are from simulations carried out by us, the squares are from Fig. 2 in Ref. 关5兴by taking into account that the quantity Cpp defined there is related to Fby Cpp⫽共1⫺ ␣ 兲⫺1 冋 C 冑 a共 ␣ 兲 ␴ Nz 册 1/2 .共54兲 Equation 共53兲predicts Cpp⫽ 冑 2 ␲ . All the reported simulation data in the figure correspond to the dilute fluidized regime, i.e., for high enough velocities of the vibrating wall, so that Fhas already reached a steady value that does not depend any more on vw. For smaller values of the velocity, F is an increasing function of it 关5兴. The value of ␰ 0has been varied by modifying the number of particles Nz. It is seen that Eq. 共50兲reproduces fairly well the simulation results, providing a definitely better approximation than Eq. 共53兲.In fact, the dependence of Cpp on the size of the system, measured by Nz, was already realized by McNamara and Luding 关5兴. Let us also stress that the asymptotic behavior for large systems, Eq. 共51兲, is only accurate for quite large values of ␰ 0and, therefore, it is not very useful in practice. We have also computed Fexact by using the values of A and Bresulting from the fitting of the minimum of the temperature profile obtained in the simulations, as discussed in Sec. III. The results 共not shown兲always lie between the simulation symbols and the curve Fapprox . The discrepancies between Fapprox and the simulation results increase as the value of the coefficient of normal restitution ␣ decreases. FIG. 7. Dimensionless quantity Fdefined in the text for ␣ ⫽0.95 as a function of the parameter ␰ 0. The continuous line is the approximated expression derived in the text 关Eq. 共50兲兴, the dashed line the prediction by Kumaran 关Eq. 共53兲兴, and the horizontal dotted line, the asymptotic value for ␰ 0→⬁. The symbols are from molecular dynamics simulations, as discussed in the main text. J. JAVIER BREY, M. J. RUIZ-MONTERO, AND F. MORENO PHYSICAL REVIEW E 63 061305 061305-8 The reason is that the energy dissipated in the region with a positive temperature slope, which is neglected in Eq. 共50兲, becomes more relevant as the inelasticity of the system increases. This is confirmed by the fact that a much better agreement is found if the expression of Fexact , with Aand B obtained from the simulation, is used. On the other hand, it is still true that, for large vw,Freaches a value that only depends on ␰ 0and ␣ , then indicating the existence of a scaling law. This scaling is not trivially seen in the expression of Fexact that depends on both constant Aand Bin a nontrivial manner. Nevertheless, it must be realized that Aand Bare not in fact independent. They must be determined from the same boundary conditions specifying the vibrating wall, although their calculation actually requires considering the free particle region. If there is a proportionality relationship between Aand B, it is easily seen that the expression of Fexact turns out to be independent of them. VI. CONCLUSIONS In this paper, we have studied a fluidized granular system submitted to an external force of the gravitational type. The general conclusion we have reached is that the hydrodynamic description provided by the 共inelastic兲Navier-Stokes equations is able to explain what is observed in molecular dynamics simulations and also in experiments. We have focused on open systems, and showed that the presence of a collisionless regime in the very high region of the gas must be taken into account when introducing the matching conditions between the bulk of the granular medium and the boundaries. Now we summarize the most important results. 共a兲The temperature profile as a function of the height presents a minimum, increasing monotonically afterwards. The minimum lies in the hydrodynamic region, but when interpreting this result it must be realized that the hydrodynamic description is not valid when the density becomes too small. 共b兲The density profile presents a maximum when the system has a large enough number of particles. This maximum is not associated, in principle, to any clustering hydrodynamic instability, but follows directly from the NavierStokes equations. 共c兲The position of the density maximum is quite accurately only determined by the coefficient of restitution of the system, being independent of the number of particles and the way in which the system is being vibrated. 共d兲An accurate description of the energy dissipated in collisions requires considering the nonuniformity of the hydrodynamic fields. The results obtained by using the exact hydrodynamic profiles derived from the Navier-Stokes equations are in better agreement with molecular dynamics simulations than those using a uniform temperature and an exponentially decreasing density. 共e兲The approximations used in some previous works have been obtained as limiting approximations of the more general results obtained here. This also applies to the scaling behavior predicted by some authors. A sound justification of scaling laws can only follow from a detailed analysis of both boundary layers, the one next to the vibrating wall and that associated with the transition to the free-particle flow. Nevertheless, we have observed in the simulations that the hydrodynamic fields seem to scale with the velocity of the vibrating wall, as already found in Ref. 关4兴 The range of validity of the analysis we have carried out deserves some comments. We have verified that there is a reasonable good agreement between the theoretical predictions derived here and the molecular dynamics results for ␣ ⬎0.9. For smaller values of the coefficient of restitution, the discrepancies become important and they increase very rapidly as ␣ decreases. There are two main related reasons that restrict a priori the applicability of our theory to the small inelasticity range. For large inelasticity, the gradients become very large and the Navier-Stokes approximation fails. Moreover, the density in the vicinity of its maximum becomes very high so that the low density hydrodynamic equations should be substituted by equations more accurate for dense granular fluids. ACKNOWLEDGMENT This research has been partially supported by the Direccio ´n General de Investigacio ´n Cientı ´ficayTe ´cnica 共Spain兲 through Grant No. PB98-1124, APPENDIX A In this Appendix, the explicit expressions for the quantities appearing in Eqs. 共4兲and 共5兲are given 关16,17兴. The Boltzmann elastic values for the shear viscosity and thermal conductivity are ␩ 0⫽2⫹d 8⌫共d/2兲 ␲ ⫺(d⫺1)/2共mkBT兲1/2 ␴ ⫺(d⫺1),共A1兲 ␬ 0⫽d共d⫹2兲2 16共d⫺1兲⌫共d/2兲 ␲ ⫺(d⫺1)/2kB 冉 kBT m 冊 1/2 ␴ ⫺(d⫺1), 共A2兲 while the dimensionless functions have the form ␩ *共 ␣ 兲⫽ 冋 ␯ 1 *共 ␣ 兲⫺ ␨ *共 ␣ 兲 2 册 ⫺1 ,共A3兲 ␬ *共 ␣ 兲⫽ 冋 ␯ 2 *共 ␣ 兲⫺2d d⫺1 ␨ *共 ␣ 兲 册 ⫺1 关1⫹c*共 ␣ 兲兴,共A4兲 ␮ *共 ␣ 兲⫽2 ␨ *共 ␣ 兲 冋 ␬ *共 ␣ 兲⫹共d⫺1兲c*共 ␣ 兲 2d ␨ *共 ␣ 兲 册 ⫻ 冋 2共d⫺1兲 d ␯ 2 *共 ␣ 兲⫺3 ␨ *共 ␣ 兲 册 ⫺1 ,共A5兲 ␨ *共 ␣ 兲⫽2⫹d 4d共1⫺ ␣ 2兲 冋 1⫹3 32c*共 ␣ 兲 册 .共A6兲 HYDRODYNAMICS OF AN OPEN VIBRATED GRANULAR SYSTEM PHYSICAL REVIEW E 63 061305 061305-9