Vibrated granular gas confined by a piston

Vibrated granular gas confined by a piston J. Javier Brey and M. J. Ruiz-Montero Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain 共Received 26 January 2009; published 23 March 2009兲 The steady state of a vibrated granular gas confined by a movable piston on top is discussed. Particular attention is given to the hydrodynamic boundary conditions to be used when solving the inelastic Navier- Stokes equations. The relevance of an exact general condition relating the grain fluxes approaching and moving away from each of the walls is emphasized. It is shown how it can be used to get a consistent hydrodynamic description of the boundaries. The obtained expressions for the fields do not contain any undetermined parameter. Comparison of the theoretical predictions with molecular-dynamics simulation results is carried out, and a good agreement is observed for low density and not too large inelasticity. A practical way of introducing small finite density corrections to the dilute limit theory is proposed to improve the accuracy of the theory. DOI: 10.1103/PhysRevE.79.031305 PACS number共s兲: 45.70.⫺n, 47.70.Nd, 51.10.⫹y I. INTRODUCTION The behavior of fluidized granular systems resembles in many cases that of ordinary molecular fluids. Actually, it is by now well established that generalized Navier-Stokes equations describe quite accurately many of the experimental and numerical features of granular flows, especially at low density and small inelasticity 关1–3兴. The justification for this fluidlike description, and derivation of theoretical predictions for the transport coefficients appearing in the equations, has been intensively studied for some time. Quite often, idealized systems of inelastic hard spheres or disks have been considered. For monodisperse models of this kind, the studies carried out include the derivation of the hydrodynamic equations to Navier-Stokes order by using kinetic theory methods 关4–9兴, direct Monte Carlo simulation of 共inelastic兲 kinetic equations 关10,11兴, molecular-dynamics 共MD兲simulations 关12兴, and, very recently, linear-response theory 关13,14兴. Interactions between grains are inherently inelastic. As a consequence, the kinetic energy of isolated granular systems decreases monotonically in time, and in order to keep them fluidized, it is necessary to continuously supply energy to them. A prototypical way of doing it is by vibrating one of the walls of the container, usually the one at the bottom. Also often, the interest focuses on the bulk properties of the system, i.e., on the behavior of the system far enough from the walls, where the governing laws are expected to be independent of the details of the boundaries. Then, the most appropriate possible way of vibration for this purpose is chosen. In many situations of interest, the criterion for this choice is twofold: simplicity and avoiding undesired effects such as the induction of propagating waves into the system. These goals are formally achieved in the limits of very high frequency and very small amplitude; the former as compared with the typical relaxation frequency of the granular fluid next to the wall and the latter with its mean free path. Additional simplifications occur if collisions of the grains with the wall are taken as elastic and if the wall is assumed to move with a sawtooth profile 关15,16兴. A kind of idealized walls often used in kinetic theory and particle simulations are the so-called thermal walls. By definition, the flux of after-collision particles leaving it corresponds to a Maxwellian flux with the temperature parameter characterizing the wall 关17兴. Therefore, the shape of the velocity distribution of the particles moving away from the thermal wall is independent of the distribution of the ingoing particles. It is evident that thermal walls also provide a mechanism to compensate the energy dissipated in inelastic collisions once the granular gas in the vicinity of the wall tends to have a temperature smaller than the one of the wall. It is then not surprising that thermal walls have also been extensively used in the literature of fluidized granular gases since a decade ago 关18,19兴. Nevertheless, it must be stressed that it is not at all evident that thermal walls correspond to any limit of a rapidly vibrating plate. Actually, it has been shown that in some cases the stability of systems driven by a thermal wall and by a wall vibrating in a sawtooth way can differ 关20兴. The granular system considered in this paper is confined between two parallel plates in the presence of a gravity force. The mission of the one at the bottom of the system is to fluidize the granular medium, as in the previous studies mentioned above. On the other hand, the wall on the top is floating, in the sense that it can move in the vertical direction, being supported by the granular fluid below it. As a consequence, the position and motion of the upper wall is interrelated with the state of the granular media below it, and the boundary conditions to be imposed to the hydrodynamic equations, following from the interaction between the gas and the piston, must be determined in a self-consistent way. Here, the steady state eventually reached by the granular gas between the two plates will be investigated using the hydrodynamic description provided by the Navier-Stokes equations for inelastic hard spheres. Collisions of the grains with the piston will be modeled as hard inelastic collisions. This defines in a deterministic way the mechanical interaction between the particles and the movable wall. The question addressed afterward is how to translate it into an appropriate boundary condition to be used in the context of hydrodynamics. Here this will be done by means of an intermediate stage in which an exact boundary condition for the kinetic theory description of the system is formulated. This condition relates the distribution functions of particles leaving the wall and approaching it by expressing the conservation of the particles flux at the piston. PHYSICAL REVIEW E 79, 031305 共2009兲 1539-3755/2009/79共3兲/031305共13兲©2009 The American Physical Society031305-1 An additional boundary condition is required to determine the hydrodynamic fields. It can be obtained, in equivalent ways, from the energy flux at the vibrating bottom wall or from the global balance of energy in the system. In this way, the theoretical prediction is completed, and explicit expressions for the fields with no adjustable parameters are derived. In order to verify the accuracy of this description, the predictions are compared with molecular-dynamics simulation results. As in other steady states of granular gases, the range of applicability of the theory is restricted to values of the restitution coefficient for gas particle collisions close to unity because of the coupling between inelasticity and gradients. Under these conditions, reasonable agreement between theory and simulations is observed in the bulk of the system, i.e., outside the kinetic boundary layers next to the walls. This confirms the validity of the hydrodynamic description, including the needed boundary conditions, to describe vibrated granular gases in quite realistic situations. The plan of the paper is as follows. In Sec. II, a previously derived 关21兴stationary solution of the inelastic Navier- Stokes equations for a vibrated dilute granular system in presence of gravity is reviewed. The main results are explicit expressions for the hydrodynamic fields of the system having two arbitrary parameters. Additionally, they involve the height of the system. The results are particularized in Sec. III for the granular gas between two plates described above. Also, the boundary effects following from the interaction between the gas and the piston on the top are formulated as a condition for the gas-piston distribution function at contact. In the same section and in Appendix B, it is discussed why this is the appropriate starting point to derive the hydrodynamic boundary condition and to introduce self-consistent approximations, such as those used in Sec. IV. An analysis along the same lines of the hydrodynamic boundary effects due to the vibrating wall at the bottom is presented in Appendix C. Also in Sec. IV, the derived boundary conditions are used to identify the arbitrary constants in the hydrodynamic profiles derived in Sec. II. There are no adjustable parameters in these expressions. The comparison of the obtained theoretical predictions with molecular-dynamics simulation results is carried out in Sec. Vfor two-dimensional systems. It includes both the detailed description provided by the hydrodynamic profiles and also some global properties, such as the average position of the piston and the balance of the total energy of the system. A fairly good agreement is observed, especially if some 共small兲finite density effects are partially incorporated into the hydrodynamic description, through the equation of state of the gas. The paper ends with a short summary and some general comments. II. GENERAL ONE-DIMENSIONAL SOLUTION In this section, some of the results already discussed in Ref. 关21兴will be shortly reviewed and summarized for the sake of completeness. The system considered is a dilute granular gas composed of Nequal smooth inelastic hard spheres 共d=3兲or disks 共d=2兲of mass mand diameter ␴ . The position and velocity of grain iwill be denoted by riand vi, respectively. The inelasticity of collisions between grains is modeled by means of a constant velocity-independent coefficient of normal restitution ␣ , defined in the interval 0⬍ ␣ ⱕ1. There is an external gravitational field acting on the system so that each particle is submitted to a force −mg0e ˆz, where g0is a positive constant and e ˆzis the unit vector in the positive direction of the zaxis. For steady states with vanishing macroscopic flow and gradients only in the direction of the external field, i.e., the z axis, the inelastic hydrodynamic Navier-Stokes equations of this system reduce to 关21兴 ⳵ p ⳵ z=−nmg0,共1兲 2 nd ⳵ ⳵ z 冉 ␬ ⳵ T ⳵ z+ ␮ ⳵ n ⳵ z 冊 −T ␨ 共0兲=0. 共2兲 Here, n共z兲is the local number of particles density, T共z兲is the local granular temperature, and p共z兲=n共z兲T共z兲is the pressure. The temperature is defined from the kinetic energy in the usual way but with the Boltzmann constant set equal to unity. Moreover, ␬ is the thermal heat conductivity and ␮ is the diffusive heat conductivity that is peculiar of granular systems. More specifically, the generalized Fourier law giving the heat flux qzin the system is qz=− ␬ ⳵ T ⳵ z− ␮ ⳵ n ⳵ z.共3兲 Finally, ␨ 共0兲共z兲is the cooling rate accounting for the energy dissipated in collisions. Upon deriving Eq. 共1兲, use has been made of the local equation p共z兲=n共z兲T共z兲, where pis the hydrodynamic pressure, valid in the low-density limit. The expressions of the transport coefficients and the cooling rate appearing in the above expressions can be written in the form 关7,11兴 ␬ 共 ␣ ,T兲= ␬ ⴱ共 ␣ 兲 ␬ 0共T兲,共4兲 ␮ 共 ␣ ,T兲= ␮ ⴱ共 ␣ 兲 ␮ 0共T兲,共5兲 ␨ 共0兲共 ␣ ,T兲= ␨ ⴱ共 ␣ 兲p ␩ 0共T兲,共6兲 with ␬ 0共T兲and ␩ 0共T兲being the elastic 共 ␣ =1兲values of the thermal heat conductivity and the shear viscosity, ␬ 0共T兲=d共d+2兲2⌫共d/2兲 16共d−1兲 ␲ d−1 2 冉 T m 冊 1/2 ␴ −共d−1兲,共7兲 ␩ 0共T兲=共d+2兲⌫共d/2兲 8 ␲ d−1 2 共mT兲1/2 ␴ −共d−1兲,共8兲 respectively, and ␮ 0共T兲=T ␬ 0共T兲 n.共9兲 The dimensionless quantities ␬ ⴱ, ␮ ⴱ, and ␨ ⴱare given in Appendix A. They only depend on the restitution coefficient ␣ . J. JAVIER BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 79, 031305 共2009兲 031305-2 Equation 共2兲shows the physically evident feature that, in this state, hydrodynamic gradients are induced by the inelasticity through the cooling rate. Consequently, a restriction to small gradients, as it is the case in the Navier-Stokes approximation used above, also implies a limitation on the value of ␣ for which the theory can be expected to apply. The interval of values of this parameter for which the theory actually provides an accurate description is very hard to determine a priori. The system is supposed to be confined between two parallel walls located at z=0 and z=L, respectively. The nature of this two walls will be specified and discussed later on. It is convenient to introduce a dimensionless length scale ␰ by 关21兴 ␰ =冑a共 ␣ 兲 冋 冕 z L dz⬘1 ␭共z⬘兲+ ␴ d−1pL mg0 册 ,共10兲 where ␭共z兲⬅关 ␴ d−1n共z兲兴−1 共11兲 is proportional to the local mean free path, pL⬅p共z=L兲is the pressure of the gas next to the wall located at z=L, and a共 ␣ 兲⬅ 32共d−1兲 ␲ d−1 ␨ ⴱ共 ␣ 兲 共d+2兲3⌫2共d/2兲关 ␬ ⴱ共 ␣ 兲− ␮ ⴱ共 ␣ 兲兴.共12兲 The ␰ coordinate is a monotonic decreasing function of z varying between ␰ M⬅ ␰ 共z=0兲=冑a共 ␣ 兲 ␴ d−1 冉 Nz+pL mg0 冊 共13兲 and ␰ m⬅ ␰ 共z=L兲=冑a共 ␣ 兲 ␴ d−1pL mg0 .共14兲 In Eq. 共13兲,Nzdenotes the number of particles in the system per unit of section W共length or area兲perpendicular to the external field, Nz⬅N/W. It must be noted that the variation interval of ␰ depends on pLand Nzbut not on the value of L. The physical meaning of ␰ can be illustrated by realizing that it is proportional to the local pressure at the corresponding height z. This follows from Eq. 共1兲that leads to p共z兲=mg0 冕 z L dz⬘n共z⬘兲+pL=mg0 冑a共 ␣ 兲 ␴ d−1 ␰ .共15兲 By substituting Eq. 共1兲into Eq. 共2兲and doing the change of variable defined in Eq. 共10兲, a closed equation for the temperature is obtained, ␰ ⳵ 2T1/2 ⳵␰ 2+b共 ␣ 兲 ⳵ T1/2 ⳵␰ − ␰ T1/2=0, 共16兲 with b共 ␣ 兲⬅ 2 ␬ ⴱ− ␮ ⴱ 2共 ␬ ⴱ− ␮ ⴱ兲.共17兲 The general solution of the above differential equation reads as 关22兴 T1/2共 ␰ 兲=A ␰ − ␯ I ␯ 共 ␰ 兲+B ␰ − ␯ K ␯ 共 ␰ 兲,共18兲 where Aand Bare constants to be identified from the boundary conditions, ␯ 共 ␣ 兲⬅b共 ␣ 兲−1 2= ␮ ⴱ 4共 ␬ ⴱ− ␮ ⴱ兲,共19兲 and I ␯ and K ␯ are the modified Bessel functions of first and second kind, respectively 关22兴. Therefore, for the system being considered, the pressure and temperature profiles are given by Eqs. 共15兲and 共18兲, respectively. Consequently, the density profile is n共 ␰ 兲=p共 ␰ 兲 T共 ␰ 兲=mg0 ␰ 冑a共 ␣ 兲 ␴ d−1关A ␰ − ␯ I ␯ 共 ␰ 兲+B ␰ − ␯ K ␯ 共 ␰ 兲兴2. 共20兲 Finally, the transformation from the ␰ coordinate to the zone is given by z=1 冑a共 ␣ 兲 ␴ d−1 冕 ␰ ␰ Md ␰ ⬘ n共 ␰ ⬘兲共21兲 that follows from Eq. 共10兲. It is worth to mention that the presence of the diffusive heat conductivity ␮ in the above expressions is not at all irrelevant. Predictions implied by its existence have been checked both by particle simulations and experimentally 关23–26兴. To proceed any further, the boundary conditions of the system at the top 共z=L兲and the bottom 共z=0兲must be specified. This will be done in Sec. III. III. CLOSED SYSTEM WITH A PISTON: THE KINETIC BOUNDARY CONDITION In order to maintain the system fluidized, it will be assumed that energy is being continuously supplied to it through the wall located at the bottom and that this is achieved by vibrating it. The simplest possible way of vibration will be considered here, namely, with a sawtooth velocity profile. This means that all the particles colliding with the wall find it with the same upward velocity vW关15,16兴. Moreover, the amplitude of vibration of this wall is taken much smaller than the mean free path of the grains in its vicinity. As a consequence, the position of the wall can be taken in practice as fixed at z=0 with very good accuracy. Finally, since the main reason to introduce this vibrating wall is to keep the granular matter fluidized, collisions of particles with it will be considered as elastic for the sake of simplicity. Of course, all the above corresponds to a very idealized wall that cannot be fully implemented in actual experiments. Next, the upper boundary condition must be specified. The case of an open system 共L→⬁兲was studied in Ref. 关21兴. Here, a different physical situation will be investigated. It will be considered that there is a movable lid or piston on top of the gas, as illustrated in Fig. 1. The piston has a finite mass Mcontrary to the vibrating wall at the bottom that is taken infinitely massive. The piston can only move in the z direction, remaining always perpendicular to it, i.e., parallel to the bottom wall. Its position and velocity will be denoted by Zand Vz, respectively, so that Lcorresponds to the aver- VIBRATED GRANULAR GAS CONFINED BY A PISTON PHYSICAL REVIEW E 79, 031305 共2009兲 031305-3 age value of Zin the steady state. There is no friction between the piston and the lateral walls of the container. Collisions of particles with the piston on the top are smooth and inelastic, with a velocity-independent coefficient of normal restitution ␣ P, where 0⬍ ␣ Pⱕ1. Therefore, the vector component, v⬜, of the velocity of the particle perpendicular to the zaxis remains unchanged in a collision with the piston, v⬜ ⬘=v⬜,共22兲 where the prime is used here and henceforth to denote aftercollision quantities. On the other hand, when a particle with a component vzof the velocity collides with the piston being Vzthe velocity of the latter, these values change instantaneously to vz ⬘=vz−M m+M共1+ ␣ P兲共vz−Vz兲,共23兲 Vz ⬘=Vz+m m+M共1+ ␣ P兲共vz−Vz兲.共24兲 Therefore, in the collision the total momentum is conserved, the relative velocity gz⬅vz−Vzchanges to gz ⬘=− ␣ Pgz, and there is a variation in the total kinetic energy given by ⌬E=− mM 2共m+M兲共1− ␣ P 2兲gz 2.共25兲 For ␣ P⬍1 there is a loss of energy in every collision between a particle and the piston. The change in the zcomponent of the momentum Pzof the piston in a collision is ⌬Pz⬅M共Vz ⬘−Vz兲=mM m+M共1+ ␣ P兲gz.共26兲 Since a collision of a particle with the piston is only possible if vz⬎Vz,itis⌬Pz⬎0 for all the collisions, indicating that momentum is continuously transferred from the gas to the piston. A relevant relationship to be used in the following is dv⬘dVz ⬘= ␣ PdvdVz.共27兲 For later use, it is convenient to consider also the so-called restituting collision corresponding to the velocities vand Vz. It is defined by the velocities vⴱand Vz ⴱleading as a consequence of a collision to vand Vz. Their expressions are obtained directly by inverting Eqs. 共22兲–共24兲, v⬜ ⴱ=v⬜,共28兲 vz ⴱ=vz−M m+M 1+ ␣ P ␣ P gz,共29兲 Vz ⴱ=Vz+m m+M 1+ ␣ P ␣ P gz.共30兲 Also, it is gz ⴱ= ␣ P −1gzand dvⴱdVz ⴱ= ␣ P −1dvdVz.共31兲 The question now is how to translate the above collision rules into one or more boundary conditions for the description of the granular gas below the piston. To discuss this point in some detail, let us introduce the two-body distribution function for the piston and the gas, ⌽共x,Z,Vz,t兲, defined for an arbitrary state as ⌽共x,Z,Vz,t兲=N 冕 dx2¯ 冕 dxN ␳ 共x,x2, ...,xN,Z,VZ,t兲, 共32兲 where ␳ 共x1,x2,...,xN,Z,VZ,t兲, with xi⬅兵ri,vi其, is the distribution function for the system composed by the piston and the Ngrains at time t, and x⬅兵r,v其. Therefore, ⌽共x,Z,Vz,t兲 is proportional to the probability density of finding the piston at height Zwith velocity Vzand a grain at position rwith velocity vat time t. It is normalized as 冕 dx 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVz⌽共x,Z,Vz,t兲=N.共33兲 The one-particle distribution function of the gas, f共x,t兲, can be obtained from ⌽by integration over the piston position and velocity, f共x,t兲= 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVz⌽共x,Z,Vz,t兲.共34兲 Similarly, the probability distribution for the piston, F共Z,Vz,t兲, is given by F共Z,Vz,t兲=1 N 冕 dx⌽共x,Z,Vz,t兲,共35兲 and it is normalized to unity. No reference to any particular state is involved in the above definitions. For initial conditions in which all the particles are located below the piston, the two-body distribution ⌽at arbitrary later times can be expressed in the form M vw Z g o W FIG. 1. Sketch of the system studied in this paper. J. JAVIER BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 79, 031305 共2009兲 031305-4 ⌽共x,Z,Vz,t兲=⌰共Z−z兲⌽0共x,Z,Vz,t兲,共36兲 where ⌰共x兲is the Heaviside step function defined by ⌰共x兲 =1 for xⱖ0 and ⌰共x兲=0 for x⬍0. The function ⌽0共x,Z,Vz,t兲is not defined by Eq. 共36兲for z⬎Z, and it can be considered as being regular everywhere as well as its derivatives without restriction. In order to discuss the form of ⌽0when the position of the particle is taken next to the piston, it is useful to decompose it in the form ⌽0共x,Z,Vz,t兲=⌽+共x,Z,Vz,t兲+⌽−共x,Z,Vz,t兲,共37兲 with ⌽+共x,Z,Vz,t兲⬅⌰共gz兲⌽0共x,Z,Vz,t兲共38兲 and ⌽−共x,Z,Vz,t兲⬅⌰共−gz兲⌽0共x,Z,Vz,t兲.共39兲 In the last expression, it is ⌰共−x兲⬅1−⌰共x兲. Conservation of the flux of particles at the piston implies that ⌽+共r,vⴱ,Z,Vz ⴱ,t兲dvⴱdVz ⴱ兩gz ⴱ兩 ␦ 共Z−z兲 =⌽−共r,v,Z,Vz,t兲dvdVz兩gz兩 ␦ 共Z−z兲,共40兲 where vⴱand Vz ⴱare the restituting velocities defined by Eqs. 共28兲–共30兲. This is an exact relationship, valid for arbitrary density, that can be derived also by starting from the evolution equation for ⌽共x,Z,Vz,t兲following from the Liouville equation for ␳ 共x1,x2,...,xN,Z,VZ,t兲. Equation 共40兲is obtained by isolating the singular terms at z=Z关27兴. When Eqs. 共28兲–共31兲are employed into Eq. 共40兲, it can be reduced to ⌽−共r,v,Z,Vz,t兲 ␦ 共Z−z兲= ␣ P −2⌽+共r,vⴱ,Z,Vz ⴱ,t兲 ␦ 共Z−z兲. 共41兲 Combination of Eqs. 共37兲and 共41兲yields ⌽0共x,Z,Vz,t兲 ␦ 共Z−z兲=共1+ ␣ P −2bP −1兲⌽+共x,Z,Vz,t兲 ␦ 共Z−z兲. 共42兲 Here bP −1 is an operator acting on the velocities vand Vzto its right, replacing them by the precollisional values given by Eqs. 共28兲–共30兲. Equation 共42兲shows that ⌽0is fully determined at z=Zif ⌽+is known at the same position. Consequently, upon introducing simplifications or approximations on the value of ⌽0 at the piston, i.e., on ⌽0 ␦ 共Z−z兲, they must refer only to either ⌽+ ␦ 共Z−z兲or ⌽− ␦ 共Z−z兲. Otherwise, the exact relationship given by Eq. 共42兲may be violated. A natural question in this context is how relevant is that relation in practice. In other words, is any fundamental physical property possibly lost if Eq. 共42兲is not satisfied in a given approximate description? The answer to this question is affirmative as it will be seen in Sec. IV. In the limit of a very dilute gas, a simple approximation similar to the one leading to the Boltzmann equation 关27兴is to assume that ⌽+共x,Z,Vz,t兲 ␦ 共Z−z兲=f共x,t兲F共Z,Vz,t兲⌰共gz兲 ␦ 共Z−z兲, 共43兲 therefore neglecting all the correlations between the piston and the particles colliding with it before the collision. Using the above approximation into Eq. 共42兲, it is found that ⌽0共x,Z,Vz,t兲 ␦ 共Z−z兲=f共x,t兲F共Z,Vz,t兲 ␦ 共Z−z兲 +⌰共−gz兲共 ␣ P −2bP −1 −1兲 ⫻f共x,t兲F共Z,Vz,t兲 ␦ 共Z−z兲.共44兲 Upon deriving this equation use has been made of the relations bP −1⌰共gz兲=⌰共gz ⴱ兲bP −1 =⌰共−gz兲bP −1,共45兲 following from the definition of bP −1 and the relationship between gzand gz ⴱ. The physical meaning of Eq. 共44兲is evident: correlations between the piston and the particles in its neighborhood are created by the collisions. In this section, only the kinetic boundary condition for the piston at the top of the system has been considered. The analysis of the vibrating wall located at the bottom is much simpler, and it will be discussed later on. IV. HYDRODYNAMIC BOUNDARY CONDITIONS In the following, attention will be restricted to the macroscopic steady state described in Secs. II and III. Then, the height of the system Lthere corresponds to the average position of the piston, Z ¯ , in the kinetic theory description, while the average velocity of the piston is required to vanish; V ¯ z⬅1 N 冕 dx 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVzVz⌽st共x,Z,Vz兲 = 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVzVzFst共Z,Vz兲=0, 共46兲 where the indices st are used to refer to properties of the system in the steady state. Although the above property can be accomplished in other ways, here the simplifying assumption that ⌽st共x,Z,Vz兲and, therefore, Fst共Z,Vz兲are even functions of Vzis made. This is consistent with the results from molecular-dynamics simulations to be reported later on. Also, it must be kept in mind that for the state we are considering, ⌽st depends on the position of the particles only through the coordinate zalthough it is not made explicit in the notation. It is interesting and illuminating to compute the force F that the granular gas makes on the piston in this state. Because of symmetry reasons, the force only has zcomponent, that is computed in Appendix B with the result Fz=WnLTL,z,共47兲 where nLis the number density of the granular gas next to the piston, i.e., nL=n共z→L兲, and TL,zis a temperature parameter of the gas in the same region defined from the zcomponent of the velocity or, equivalently, proportional to the zz component of the pressure tensor, VIBRATED GRANULAR GAS CONFINED BY A PISTON PHYSICAL REVIEW E 79, 031305 共2009兲 031305-5 TL,z⬅1 nL 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVzmvz 2⌽0,st共x,Z,Vz兲 ␦ 共Z−z兲. 共48兲 In the Navier-Stokes approximation used in Sec. I,TL,zcoincides with the temperature of the gas at the piston, TL, and Eq. 共47兲is the expected result since it agrees with the hydrodynamic interpretation of the pressure tensor. In its derivation, a crucial role is played by the kinetic boundary condition given in Eq. 共41兲, as it is shown in Appendix B. Still more, the specific form of ⌽+,st is not relevant as long as ⌽−,st be consistently derived from it. Otherwise, the pressure of the dilute gas defined as the product of the local density times the local temperature would not agree with the scalar defined from the hydrodynamic pressure tensor, characterizing the internal forces in the fluid. In order to develop a consistent theory, it is then convenient to express the hydrodynamic fields of the granular gas in the vicinity of the piston in terms of only ⌽+,st instead of the complete distribution ⌽0,st. The expression of the number density nLintroduced above in terms of ⌽0,st reads as nL⬅ 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVz⌽0,st共x,Z,Vz兲 ␦ 共Z−z兲. 共49兲 By means of the decomposition formulated in Eq. 共37兲and using boundary condition 共41兲, the above expression can be put in the form nL=1+ ␣ P ␣ P nL 共+兲,共50兲 with nL 共+兲⬅ 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVz⌽+,st共x,Z,Vz兲 ␦ 共Z−z兲. 共51兲 Also, it is easily verified that the local velocity flow vanishes at the piston in the steady state as it should. Finally, consider the temperature of the gas next to the piston, TL, given by d 2nLTL⬅ 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVz mv2 2⌽0,st共x,Z,Vz兲 ␦ 共Z−z兲. 共52兲 It is convenient to distinguish between the perpendicular and zcontributions to this temperature, TL=1 dTL,z+d−1 dTL,⬜,共53兲 where TL,zis defined in Eq. 共48兲and TL,⬜⬅1 共d−1兲nL 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVzmv⬜ 2 ⫻⌽0,st共x,Z,Vz兲 ␦ 共Z−z兲.共54兲 Boundary condition 共41兲leads directly to TL,⬜=TL,⬜ 共+兲,共55兲 with TL,⬜ 共+兲⬅1 共d−1兲nL 共+兲 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVzmv⬜ 2 ⫻⌽+,st共x,Z,Vz兲 ␦ 共Z−z兲.共56兲 A more involved calculation gives TL,z= ␣ P共M−m兲 M− ␣ PmTL,z 共+兲−Hz 共+兲,共57兲 where TL,z 共+兲⬅1 nL 共+兲 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVzmvz 2⌽+,st共x,Z,Vz兲 ␦ 共Z −z兲,共58兲 Hz 共+兲⬅2 ␣ PMm 共M− ␣ Pm兲nL 共+兲 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVzvzVz ⫻⌽+,st共x,Z,Vz兲 ␦ 共Z−z兲.共59兲 At this point, a simplifying hypothesis on the precollisional two-body distribution at contact is made. The associated marginal velocity distribution is approximated by a product of Gaussian distributions, namely, it is assumed that 冕 dz 冕 0 ⬁ dZ⌽+,st共x,Z,Vz兲 ␦ 共Z−z兲=fst 共+兲共v兲Pst共VZ兲⌰共gz兲, 共60兲 where fst 共+兲共v兲=2nL 共+兲 ␸ MB共v⬜兲 冉 m 2 ␲ TL,z 共+兲 冊 1/2 e−mvz 2/2TL,z 共+兲,共61兲 ␸ MB共v⬜兲= 冉 m 2 ␲ TL,⬜ 冊 共d−1兲/2 e−mv⬜ 2/2TL,⬜,共62兲 Pst共Vz兲= 冉 M 2 ␲ TP 冊 1/2 e−MVz 2/2TP.共63兲 Although in the same spirit, Eq. 共60兲is a somewhat stronger assumption than the particularization of Eq. 共43兲for the steady state under consideration. The presence of nL 共+兲and TL,z 共+兲in Eq. 共61兲and of TL,⬜in Eq. 共62兲is required by consistency with the previous results in this section. Moreover, it will be assumed in the following that TL,⬜=TLfor the sake of simplicity. Because of Eq. 共53兲this implies that also TL,z =TLand, therefore, all the diagonal components of the pressure tensor of the gas in the vicinity of the piston are the same, consistently with the Navier-Stokes approximation. Using Eqs. 共60兲–共63兲it is obtained that HZ 共+兲=0 and, therefore, Eq. 共57兲reduces to J. JAVIER BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 79, 031305 共2009兲 031305-6 TL=TL,z= ␣ P共M−m兲 M− ␣ PmTL,z 共+兲.共64兲 Equations 共50兲and 共64兲relate the properties of the flux of grains reaching the wall with the local hydrodynamic fields of the granular gas next to the wall. They will be employed now to derive an expression for the heat flux at the piston, QL. Using standard kinetic theory arguments, this quantity can be computed from the variation in kinetic energy of the grains colliding with it, namely, it is given by the energy of the flux leaving the piston minus that of the flux reaching it. Then, it can be written as QL= 冕 dz 冕 dv 冕 0 ⬁ dZ 冕 −⬁ ⬁ dVz m 2共vz 2 −vz ⬘2兲gz⌽+,st共x,Z,Vz兲 ␦ 共Z−z兲.共65兲 This expression is easily evaluated in the Gaussian approximation given by Eqs. 共60兲–共63兲with the result QL=2 冉 2 ␲ m 冊 1/2M M+m 冉 1+ ␾ ␣ P 冊 1/2 冉 M− ␣ Pm M−m 冊 3/2 ⫻ 冋 1−共1+ ␣ P兲共1+ ␾ 兲M 2共M+m兲 册 nLTL 3/2,共66兲 where ␾ ⬅mTP MTL,z 共+兲= ␣ P共M−m兲mTP 共M− ␣ Pm兲MTL .共67兲 The calculation of the power injected into the system through the vibrating wall at the bottom is much more direct, and it is outlined in Appendix C. There, it is shown that the heat flux at this wall is given by Q0=vWp0,共68兲 with p0⬅p共z=0兲=p共 ␰ = ␰ M兲is the pressure of the granular gas in the region just above the vibrating wall. Equation 共68兲 was proposed in Refs. 关15,16兴and used many times in the literature since then. To get the value of p0, first note that Eq. 共15兲gives p0=pL+mg0Nz.共69兲 The value of pLcan be obtained by means of Eq. 共47兲after substituting TL,zby TLaccordingly with Eq. 共64兲. Mechanical equilibrium of the piston in the steady state requires that Fz=Mg0and, therefore, pL⬅nLTL=Mg0 W.共70兲 Once pLis known, the values of ␰ Mand ␰ mdefined in Eqs. 共13兲and 共14兲, respectively, can be determined. The next step is to identify the constants Aand Bintroduced in Sec. II and appearing in the expressions of the density and temperature profiles 关Eqs. 共18兲and 共20兲兴. To do so, two boundary conditions are needed. They are provided by requiring the values of QLand Q0given above to agree with the hydrodynamic expression of the heat flux qzgiven in Eq. 共3兲, particularized for each of the wall boundaries, i.e., QL= lim z→Z qz共z兲,共71兲 Q0= lim z→0 qz共z兲,共72兲 with qz共z兲=− ␬ 0 冋 共 ␬ ⴱ− ␮ ⴱ兲 ⳵ T ⳵ z− ␮ ⴱmg0 册 .共73兲 By using Eqs. 共10兲,共18兲, and 共19兲, as well as properties of the modified Bessel functions 关22兴, the above expression becomes qz共 ␰ 兲=2mg0 ␬ 0共T兲 T1/2共 ␬ ⴱ− ␮ ⴱ兲 ␰ 1− ␯ 关AI ␯ −1共 ␰ 兲−BK ␯ −1共 ␰ 兲兴. 共74兲 Making in this equation ␰ = ␰ mand equating the result to the expression of QLin Eq. 共66兲a relationship between the constants Aand Bis obtained, B A=I ␯ −1共 ␰ m兲−e共 ␣ , ␣ P兲I ␯ 共 ␰ m兲 K ␯ −1共 ␰ m兲+e共 ␣ , ␣ P兲K ␯ 共 ␰ m兲,共75兲 where e共 ␣ , ␣ P兲=16冑2共d−1兲 ␲ 共d−2兲/2 d共d+2兲2⌫共d/2兲 M M+m 冉 1+ ␾ ␣ P 冊 1/2 ⫻ 冉 M− ␣ Pm M−m 冊 3/2 冋 1−共1+ ␣ P兲共1+ ␾ 兲M 2共M+m兲 册 ⫻1 冑a共 ␣ 兲共 ␬ ⴱ− ␮ ⴱ兲.共76兲 Next, Eq. 共74兲is particularized for ␰ = ␰ Mand afterward put equal to the right-hand side of Eq. 共68兲. This leads to AI ␮ −1共 ␰ M兲−BK ␮ −1共 ␰ M兲= ␰ M ␯ 8共d−1兲 ␲ 共d−1兲/2 d共d+2兲2⌫共d/2兲 m1/2vW 冑a共 ␣ 兲共 ␬ ⴱ− ␮ ⴱ兲. 共77兲 An alternative to one of the two boundary conditions used above would be to require that the total power dissipated in the system due to the inelasticity of collisions be the same as the net heat flux injected in the system through the boundaries. This leads to a relation between Aand Bthat is a combination of Eqs. 共75兲and 共77兲, consistently with the fact that the condition follows directly from the hydrodynamic equations. This can be realized by noting that the energy balance equation in the steady state reads as 关see Eqs. 共2兲and 共3兲兴 ⳵ qz ⳵ z=−nd 2T ␨ 共0兲,共78兲 whose integration between z=0 and z=Lgives VIBRATED GRANULAR GAS CONFINED BY A PISTON PHYSICAL REVIEW E 79, 031305 共2009兲 031305-7 qz,L−qz,0 =−d 2 冕 0 L dzp共z兲 ␨ 共z兲.共79兲 The right-hand side of this equation is clearly identified as the energy dissipated in the system per unit of time and section due to the inelasticity of collisions. To close the identification of the hydrodynamic profiles in the system, a prescription to determine the parameter ␾ defined in Eq. 共67兲is needed. This in turn calls for an expression for the temperature of the piston TP. It is well known that a peculiar feature of granular gases is the violation of energy equipartition 关8,28–30兴that in the present case manifests itself by the difference between TPand TL. This has been discussed with some detail in 关31兴. For values of the restitution coefficient ␣ close to unity, a good approximation to the value of ␾ is obtained by taking TP/TLequal to unity. A relevant quantity characterizing the macroscopic state of the system is the average height of the piston, L=Z ¯ . The theoretical prediction for it is obtained from Eq. 共21兲once the values of ␰ Mand ␰ mhave been obtained and the density profile is known, L=1 冑a共 ␣ 兲 ␴ d−1 冕 ␰ m ␰ Md ␰ n共 ␰ 兲.共80兲 V. MOLECULAR-DYNAMICS SIMULATIONS To test the theoretical predictions discussed in Secs. II–IV, MD simulations of a two-dimensional 共d=2兲system of inelastic hard disks have been performed. To avoid undesired boundary effects, periodic boundary conditions with periodicity Wwere used in the direction perpendicular to the zaxis. In all the simulations to be reported in the following, the grains were initially located in a square lattice with a Gaussian velocity distribution. The initial position of the piston was slightly above the highest layer of grains. Then, the system was evolved in time accordingly with the mechanical rules governing the dynamics of the particles and the piston, and it was observed that, for values of the restitution coefficient ␣ between 0.9 and 1, a steady state with only gradients in the vertical direction and no flow field was reached. Actually, in order for the system to reach this state, its width W has to be not too large since otherwise transversal inhomogeneities of the type discussed in Refs. 关20,32兴develop into the system. The results presented below have been time averaged once the system is in the steady state and also over several trajectories. More precisely, once the stationary state was reached, its trajectory was followed for 4000 collisions per particle, and 10 different trajectories were used in each case. To carry out a systematic analysis of the theoretical predictions, it is necessary to reduce the number of parameters of the system by fixing some of them. In the present study, the values N=420 and W=70 ␴ 共Nz=6 ␴ −1兲have been used in all the simulations. The reason is that the dependence on Nz of the theoretical results follows trivially once the hydrodynamic description is assumed to hold, and analysis of the simulation results indicates that the above values are appropriate for the purposes here, in the sense that a hydrodynamic behavior is observed. Moreover, only results for values of ␣ in the interval 0.9ⱕ ␣ ⬍1 will be reported, the reason being that for smaller values of the restitution coefficient of the gas, qualitative and quantitative strong deviations from the theoretical predictions were found. This is not surprising since the Navier-Stokes approximation is expected to fail beyond the weak dissipation limit as a consequence of the coupling between gradients and inelasticity, as pointed out in Sec. II. A. Hydrodynamic profiles In Figs. 2and 3, the dimensionless pressure 共pⴱ兲and temperature 共Tⴱ兲profiles, respectively, are plotted as a function of the length scale ␰ defined in Eq. 共10兲, for two different combinations of the parameters M/m, ␣ , and ␣ P, as indicated in the figures themselves. Remind that ␰ is a decreasing function of z. The dimensionless fields are defined by Tⴱ⬅T mvW 2,pⴱ⬅pm ␴ g0 .共81兲 It is easily seen that in these units, the theoretical predictions for the hydrodynamic fields become independent of the velocity vWof the vibrating wall and also of the acceleration g0. Of course, this applies as long as the system is fluidized and its density low everywhere. Consider first the pressure field. The theoretical prediction is given by Eq. 共15兲and does not involve the parameters A and Bappearing in the temperature and density profiles. It is represented by the solid lines in Fig. 2. In the simulations, the ␰ coordinate has been measured directly by using the discrete version of Eq. 共10兲. On the other hand, the pressure at each value of zhas been obtained from the data for the density and the temperature at the same value. In the low- 0.3 0.8 1. 3 ξ 0 4 8 p *α=0.98 αP=0.8 M=92m 0.3 0.8 1.3 1.8 ξ 0 4 8 p *α=0.94 αP=0.99 M=48m FIG. 2. Dimensionless pressure 共pⴱ兲profiles for two steady states of the system of inelastic hard disks depicted in Fig. 1. The solid lines are the theoretical predictions given in the main text, while the symbols are MD simulation results. The empty circles have been computed using pⴱ=nⴱTⴱ, while the triangles have been calculated by means of Eq. 共82兲. The values of the parameters of the system are M=48m, ␣ =0.94, ␣ P=0.99, and vW=7冑g0 ␴ for the top figure, and M=92m, ␣ =0.98, ␣ P=0.8, and vW=3冑g0 ␴ for the bottom one. J. JAVIER BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 79, 031305 共2009兲 031305-8 density limit considered in the theory, it is pⴱ=nⴱTⴱ, where nⴱ⬅n ␴ vW 2/g0. The results for the pressure represented by empty circles in the figures have been computed in this way. A good agreement is observed between theory and simulation in the bulk of the system, i.e., outside the boundary layers next to the upper and lower walls. The results can be considered as satisfactory, especially taking into account that there are no adjustable parameters. Nevertheless, it is true that a systematic although small deviation is observed. Equation 共15兲is a quite general result, which only requires for its derivation the restriction to the Navier-Stokes order approximation without any particular expression of the transport coefficients or of the 共local兲equation of state. Thus it seems sensible to check whether the origin of the observed discrepancy lies in the way in which the pressure is computed from the simulation data, namely, by using the ideal gas equation of state. For this reason, the local pressure has also been calculated from the MD data by using the equation of state proposed by Grossman et al. 关18兴for hard disks at finite density, p=nTnc+n nc−n,共82兲 where nc=2 ␴ −2/冑3 is the maximum packing number density. The pressure values obtained by this procedure are represented by empty triangles in Fig. 2. Now, a much better agreement is obtained. It is concluded that density corrections are identifiable in the pressure in spite of the density being quite small. Actually, the same values of the pressure are obtained if instead of Eq. 共82兲, the second virial approximation for the equation of state of a gas of hard disks is used. The scaled temperature profiles shown in Fig. 3also exhibit a good agreement between the theoretical predictions and the simulation results. Moreover, the best fits obtained varying the parameters Aand Bin Eq. 共18兲are also plotted. Introducing density corrections to the expression for the temperature as discussed above for the pressure is far from trivial. The arguments in Sec. II leading to a closed separated equation for the temperature do not apply anymore, and the resulting differential equations do not seem to have an elementary solution. In any case, the above results clearly confirm the accuracy of both the hydrodynamic equations and the used boundary conditions. Notice that the boundary layer next to the movable piston is much narrower for the temperature than for the pressure. Similar degree of agreement has been found for all the studied combinations of parameters in the intervals 20 ⱕM/mⱕ120, 0.9ⱕ ␣ ⬍1, and 0.6ⱕ ␣ Pⱕ1. In particular, the inelasticity of the collisions between the movable piston and the grains, measured by the coefficient ␣ P, seems to affect very weakly the accuracy of the theory. It is interesting to plot also the hydrodynamic profiles in the original scale zor, equivalently, as a function of zⴱ =zg0/vW 2. In Fig. 4, the temperature profiles in the zⴱscale are shown for the same two steady states as in the previous figures. Also included are the profiles obtained with the best fitting parameters Aand B. It is observed that the agreement between theory and simulation is now worse than in the previous figures, where the spatial scale ␰ was employed. This is not at all surprising since in the transformation of the theoretical prediction from the scale ␰ to the scale z关Eq. 共10兲兴, the density profile is involved, and the ideal gas equation of state has been used. Consequently, the density effects discussed above manifest themselves in the transformation. Moreover, the scale zis defined from ␰ in a cumulative way, as an integral over the density profile so that the discrepancies increase as zdecreases. This effect is also clearly identified in Fig. 5, where the profiles for the reduced density nⴱ are given. It is also seen that the simulation results extend to larger values of zthan the theoretical predictions for the profiles. It is worth to emphasize that the maximum value of 0.3 0.8 1. 3 ξ 5 15 25 T * α=0.98 αP=0.8 M=92m 0.3 0.8 1.3 1.8 ξ 0 3 6 T * α=0.94 αP=0.99 M=48m FIG. 3. 共Color online兲Dimensionless temperature profile Tⴱfor the same two steady states as in Fig. 2. The solid line is the theoretical prediction given by Eq. 共18兲with the values of Aand B determined in Sec. IV. The empty circles are MD results and the dashed line 共red online兲is the best fit to Eq. 共18兲by varying Aand B. 01020 z* 5 15 25 T * α=0.98 αP=0.8 M=92m 0123 z* 0 2 4 6 T * α=0.94 αP=0.99 M=48m FIG. 4. 共Color online兲Temperature profiles for the same systems as in Fig. 2, but now as a function of the reduced original coordinate zⴱdefined in the text. The solid line is again the theoretical prediction with no adjustable parameters, and the symbols are the MD simulation results. The dotted-dashed lines 共green online兲 have been obtained by fitting the two parameters Aand Bappearing in the expression of the temperature field using, moreover, the equation of state 关Eq. 共82兲兴 as discussed in the main text. VIBRATED GRANULAR GAS CONFINED BY A PISTON PHYSICAL REVIEW E 79, 031305 共2009兲 031305-9