Steady state of a fluidized granular medium between two walls at the same temperature

Steady state of a fluidized granular medium between two walls at the same temperature J. J. Brey and D. Cubero Fı ´sica Teo ´rica, Facultad de Fı ´sica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain ~Received 30 June 1997; revised manuscript received 17 October 1997! The steady state of a low-density gas of inelastic hard spheres confined between two parallel walls at the same temperature is studied. Because of the dissipation in collisions, the state is not uniform but highly inhomogeneous with a nonlinear temperature profile. Direct Monte Carlo simulations show that in the nearly elastic limit the pressure is uniform, but the state exhibits anisotropy of the diagonal terms of the pressure tensor, contrary to the predictions of the Navier-Stokes equations. For larger inelasticity, the pressure becomes nonuniform. These rheological effects, peculiar to granular systems, are explained by means of a model kinetic equation based on the Boltzmann equation. The equation is solved by constructing a systematic perturbative expansion in the square root of the degree of inelasticity. The theoretical predictions compare well with the simulation results for small inelasticity, but they are in conflict for larger values of the degree of inelasticity. The analysis provides strong evidence that this is due to the asymptotic but divergent character of the expansion, similarly to what happens when the usual Chapman-Enskog expansion is applied to molecular fluids. @S1063-651X~98!10602-5# PACS number~s!: 81.05.Rm, 05.20.Dd, 51.10.1y, 47.20.2k I. INTRODUCTION The rapid flow regime of granular media is characterized by the free motion of the particles between collisions. This has led to the consideration of the ~macroscopic!grains as analogous to the particles of a molecular fluid, trying to extend to rapid granular flows the well founded theories developed for ordinary fluids out of equilibrium. In this context, a system of smooth inelastic hard spheres is considered an ideal model to study transport and relaxation in rapid granular flows, and hydrodynamiclike equations have been proposed for this system @1,2#. The equations are similar to the conventional Navier-Stokes equations, modified to account for dissipation in collisions. Nevertheless, the correct inclusion of the effects following from the lack of energy conservation is far from being trivial and a detailed derivation of the equations starting from a fundamental basis is required. Kinetic theory provides a level of description from which the validity of a hydrodynamic description, its form, and explicit expressions for the transport coefficients can be determined. The extension to a system of inelastic hard spheres of the Boltzmann equation and also of the Enskog equation is by now well established, both by using heuristic arguments @3– 5#and by starting from the Liouville equation of the system @6#. However, the complexity of the kinetic equations for inelastic systems has required the introduction of largely uncontrolled approximations when solving them. In particular, the standard Chapman-Enskog method @7#, in which an expansion in powers of the gradients of the hydrodynamic fields is carried out, cannot be directly applied to a granular fluid. Time evolution and space gradients are coupled not only among themselves, but also to dissipation and this has very significant consequences. For instance, it is known that freely evolving granular fluids do not reach a steady state, showing in addition the spontaneous formation of dense clusters @8,9#. Also, when submitted to homogeneous boundary conditions, a granular fluid can reach a steady nonuniform state, as it is the case in the situation to be considered here. Gradients are controlled by dissipation in collisions and not only by the boundary and initial conditions. In this sense, the limit of small gradients of the hydrodynamic fields is closely related to the limit of quasielastic collisions and this must be taken into account when looking for equations describing the hydrodynamic fields. To allow for a detailed and exact analysis of granular flows, it is useful to consider model kinetic equations that preserve the critical features of the original system. Here we will use a model for low-density inelastic gases @10#, obtained as an approximation of the Boltzmann equation. It has previously been applied to steady shear driven states and shown to lead to results in good agreement with simulations of the Boltzmann equation over a wide range of values of the parameter characterizing dissipation in collisions @11#. We will consider here the steady state reached by a system between two infinite parallel boundaries that are kept at constant temperature. For molecular fluids such a state is trivial since it is the Maxwellian equilibrium one. Nevertheless, for granular systems space gradients are developed in the system as a consequence of dissipation in collisions and the steady state is highly inhomogeneous. Closely related states have been studied recently. Grossman, Zhou, and BenNaim @12#have considered a two-dimensional system of hard disks in a box where one wall was kept at a fixed temperature and the other three were reflecting. The corresponding one-dimensional case had been considered previously @13#. Grossman et al. used transport equations, derived by means of heuristic arguments, which were expected to describe both highand low-density regions. Their analysis is restricted to first order in the gradients ~Navier-Stokes order! and to the quasielastic limit, leading to uniform pressure. Here we will show that the pressure is not uniform, except in the asymptotic limit of small dissipation. Steady states of inelastic systems in the presence of an energy source also have been investigated experimentally by considering a system of spherical particles rolling on a smooth rectangular surface @14#. One of the sidewalls was PHYSICAL REVIEW E FEBRUARY 1998VOLUME 57, NUMBER 2 57 1063-651X/98/57~2!/2019~11!/$15.00 2019 © 1998 The American Physical Society displaced periodically, thus supplying energy to the system. Hydrodynamic profiles and also clustering away from the energy source, similar to those obtained in Ref. @12#, were observed. Nevertheless, the effective coefficient of restitution was too far from unity to expect a quantitative agreement. An important phenomenon we have found is the presence of normal stress differences in our system. This is a well known effect for both molecular and granular sheared flows that shows up to Burnett order ~expansion of the fluxes up to second order in the gradients of the hydrodynamic fields!, but, to the best of our knowledge, has not been reported before in granular systems without velocity flow. In addition, normal stress differences are significant, in the sense that they are already clearly observed for quite small values of the dissipation parameter. This appears to be another unexpected phenomenon peculiar to rapid granular flows and associated with rheological effects and not very small dissipation. Let us point out that the components of the pressure tensor were not considered in any of the studies of related steady states mentioned above. In order to identify the origin of the anisotropy of the diagonal terms of the pressure tensor and eventually the nonuniformity of the pressure, it is necessary to go beyond the Navier-Stokes approximation. As mentioned above, this is a very hard task for a general situation because of the complexity introduced by the energy sink term @10#. Nevertheless, for the state considered here, we can exploit the fact that the gradients are determined by the dissipation parameter and the system approaches the equilibrium state in the elastic limit. Therefore, an expansion in the dissipation parameter is not only convenient because of practical reasons, but also the only consistent way of carrying out a perturbation expansion of the solution of the kinetic equation for the steady state. A similar approach has been previously used by Sela, Goldhirsch, and Noskowicz for a sheared two-dimensional granular gas @15#. Nevertheless, while in the steady sheared state considered by the above authors the only present gradient, the shear rate, is constant throughout the system, in the steady state we will deal with all the hydrodynamic fields, which are, in principle, highly nonuniform. Consequently, we have to include in our description space derivatives of order higher than one. We have also investigated the same state by means of the direct Monte Carlo simulation method @16#, which has been developed to obtain numerical solutions to the Boltzmann equation and can be also applied to the case of inelastic collisions @11,17#. The numerical results agree very well with the predictions of the model kinetic equation in the limit of very small inelasticity, while the agreement is only qualitative for larger values of the inelasticity. This is due not only to the simplification inherent to the model kinetic equation but also, and mainly, to the asymptotic character of the series expansion carried out. The indication following from our calculations is that such an expansion is divergent, as is the case of the usual Chapman-Enskog expansion @18#. The structure of the paper is as follows. In Sec. II the kinetic model is briefly reviewed and the steady state to be studied is introduced. In addition, the Navier-Stokes approximation is discussed. We believe that this is important in order to establish that it is appropriate only in the limit of very small dissipation and to motivate the series expansions carried out in Sec. IV. Direct Monte Carlo simulation results are presented in Sec. III. It is shown that normal stress differences appear for quite small values of the inelasticity and that, for slightly larger values, the pressure is not uniform as well. The series expansion solution to the model kinetic equation developed in Sec. IV correctly describes the anisotropy effect, but fails to give a correct prediction for the inhomogeneity of the pressure. The reasons for this conflict are discussed. A comparison of theory and simulation is also carried out for the one-particle distribution function. The agreement is excellent at very low dissipation and for thermal velocities. Finally, Sec. V provides a short summary and conclusions. II. THE KINETIC MODEL AND THE NAVIER-STOKES APPROXIMATION The kinetic model we will use has been introduced recently as a simplification of the Boltzmann equation. Since the details of its motivation and derivation are described elsewhere @6,10#, we give only the results here. The model kinetic equation for the one-particle distribution function f(r,v,t) of a low-density granular gas of smooth hard disks (d52) or spheres (d53) of diameter s and mass m, whose collisions are characterized by a constant coefficient of restitution a ,is ] ] tf1v•¹f52 n ~f2fl!21 nkBT~12 a 2! w ~V! v @f u f#fl, ~1! where V~r,t!5v2u~r,t!,~2! w ~V!5mV2 dkBT21, ~3! and flis the local equilibrium distribution fl~r,v,t!5n S m 2 p kBT D d/2 exp S 2mV2 2kBT D .~4! Moreover, n is an effective frequency given by n 5Cn s d21 S p kBT m D 1/2 ,~5! with Cbeing a dimensionless constant to be fixed later on and v @f u f#a source term describing the rate of dissipation in collisions. It is a functional of the distribution function defined by v @g,h#5m p ~d21!/2 s d21 8G S d13 2 D E dv1dv2 u v12v2 u 3g~v1!h~v2! 5 v @h u g#.~6! In the above expressions, nis the local number density, uis the local flow velocity, and Tis the local temperature. They are defined in terms of f(r,v,t) in the usual way, n~r,t!5 E dvf~r,v,t!,~7! 2020 57J. J. BREY AND D. CUBERO n~r,t!u~r,t!5 E dvvf ~ r,v,t ! ,~8! d 2n ~ r,t ! k B T ~ r,t ! 5 E dv1 2mV2~r,t!f~r,v,t!.~9! From Eq. ~1!the balance equations for mass, momentum, and energy are obtained by taking the appropriate velocity moments, ] n ] t1¹•~nu!50, ~10! ] u ] t1u•¹u1~nm!21¹•P50, ~11! d 2nkB ] T ] t1d 2nkBu•¹T52P:~¹u!2¹•q2~12 a 2! v . ~12! The pressure tensor Pand the heat flux qare given by P~r,t!5 E dvmVVf~v,r,t!,~13! q~r,t!5 E dvmV2 2Vf~v,r,t!.~14! Let us point out that the balance equations are preserved by the model. Not only do they have the same form as obtained from the Boltzmann equation, but also the fluxes P and qand the source term v are given by the same functionals of the distribution function. We want to investigate the properties of a steady state with no macroscopic velocity field. In addition, we consider a geometry of the system such that there are gradients in only one direction, chosen as the x axis. It is convenient to introduce a new scale sdefined by s~x!5 E 0 xdx8 n ~x8!.~15! In terms of it, particularization of Eq. ~1!for the steady state reads vx ] ] sf~s!52f~s!1fl ~s!2~12 a 2! v ~s! n~s!kBT~s! n ~s! w ~s!~v!fl ~s!. ~16! The superscript sindicates that the quantity refers to the steady state. Then the local equilibrium distribution fl (s)is particularized for u50,n5n(s)(s), and T5T(s)(s). Furthermore, when applied to this steady state, Eq. ~10!is verified identically, while Eqs. ~11!and ~12!reduce to ] ] sPxi ~s!50~17! for all iand ] ] sqx ~s!52~12 a 2! v ~s! n ~s!.~18! In the remainder of this section we are going to consider the Navier-Stokes approximation, in which the pressure tensor and heat flux are given by @6,10,19# Pij5p d ij2 h S ¹iuj1¹jui22 d d ij¹•u D ,~19! qi52l¹iT2 m ¹in,~20! where p5nkBTis the hydrostatic pressure, h the shear viscosity, lthe thermal conductivity, and m a coefficient that has no analog for elastic fluids and vanishes in the limit a 51. The values of these transport coefficients are h 52 h 0~22 g !21,l5l0~122 g !21, m 52l0T n g ~122 g !21~223 g !21,~21! with g 5~12 a 2!2 p ~d22!/2 CdG S d 2 D .~22! Here h 05p/ n and l05(d12)kBp/2 n mare the shear viscosity and thermal conductivity in the elastic limit, respectively. Thus, in the Navier-Stokes approximation we have from Eq. ~19! Pij ~s!5p~s! d ij ~23! and Eq. ~17!implies that the pressure is uniform in the steady state. Equation ~20!leads to qx ~s!52a~ g !p~s! ] ] sT~s!~s!,~24! where a~ g !5~d12!~225 g !kB 2m~122 g !~223 g !.~25! When the above expression for the heat flux is substituted into Eq. ~18!, an equation for the temperature profile is obtained, ] 2 ] s2T~s!~s!5 g d 2a~ g !.~26! In the case of elastic collisions, the right-hand side of this equation vanishes and one recovers the linear in stemperature profile characteristic of molecular gases. Nevertheless, let us note the qualitative difference between the linear and the quadratic terms in the expression of the temperature profile. While the coefficient of the former is determined by the boundary conditions, the coefficient of the latter is a given 57 2021STEADY STATE OF A FLUIDIZED GRANULAR MEDIUM . . . function of the coefficient of restitution a and does not depend on the nature of the boundaries of the system. As a consequence, the steady states of a granular fluid with vanishing flow velocity are intrinsically inhomogeneous. This is closely related to the fact that the homogeneous state of an isolated granular fluid, the homogeneous cooling state, is characterized by a monotonically decreasing in time temperature @1,9,17#. An important consequence of Eq. ~26!is that gradients and dissipation in collisions, measured by the restitution coefficient, are not independent. Therefore, retaining contributions up to a given order in one of them implies a limitation also in the range of values of the other. In order to fully specify the model, we have to fix the value of the constant Cappearing in the expression of the effective collision frequency n , Eq. ~5!. Since the heat flux is the only macroscopic flux in our system, a sensible choice appears to be to require that the model gives the same value for the Navier-Stokes thermal conductivity in the elastic limit a 51 as the Boltzmann equation. This leads to C.1 for d52 and C.32/15 for d53@7#. Grossman et al. @12#have studied a system of nearly elastic hard disks in a steady state similar to the one we are dealing with here. By using heuristic arguments, they obtain an equation for the density profile, which in the low-density limit is equivalent in our notation to T~s!1/2 ] ] x S T~s!1/2 ] ] xT~s! D 58~12 a 2! s 2p~s!2 lkB 2.~27! Here lis an unknown parameter whose value is determined by the authors by fitting molecular-dynamics simulation data obtained for a *0.99. They get l.6.03. It is easily verified that Eq. ~26!can also be written in the form of Eq. ~27!with l516~225 g ! C p ~122 g !~223 g !,~28! which for 0.99< a <1 and C51 leads to values of lin the interval 5.09&l&5.30, which are close to the value found in Ref. @12#. Nevertheless, in the present calculation, which is not restricted to the quasielastic limit, ldepends quite strongly on the value of the restitution coefficient a . III. DIRECT MONTE CARLO SIMULATION To test the theoretical predictions obtained in the preceding section, we have carried out direct Monte Carlo simulations of the Boltzmann equation for an inelastic gas. This simulation method has been applied extensively in the case of molecular gases @16#and can be easily extended to granular flows @11,17#. We refer the reader to the literature for the details of the method. We considered a dilute system of Nsmooth inelastic hard disks between two infinite plates perpendicular to the xaxis and separated by a distance L. The two plates were treated as diffusive thermal walls at a temperature TW. This is implemented in the simulation by assigning to the particles, after collision with the walls, a velocity drawn from the MaxwellBoltzmann distribution at temperature TW@20,21#. Periodic boundary conditions were applied in the ydirection, perpendicular to the thermal walls. The system was split into M 5MxMysquare cells of the same size. We took advantage of the symmetry of the system and the number of cells in the y direction, My, was smaller than in the xdirection, Mx. The initial state was taken to be homogeneous, i.e., all cells had the same number of particles Nc5N/M, and the velocities were distributed according to a Maxwellian with temperature TW. Reduced units defined by m51, kBTW51/2, and l 51, where l5@2 A 2n(0) s #21is the initial mean free path, were used. The time step over which it is assumed in the simulation that free motion and collisions are not correlated was taken to be Dt50.2. In all the simulations the values of the parameters were L530, Mx5120, My510, and Nc 530. After an initial transient time, the system reached a steady state in which the averages over different trajectories of the properties of the system became time independent. The results we will report in the following correspond to time averages over a number of trajectories once the system was in the steady state. Let us first consider the limit of a very close to unity, namely, a >0.99. In Figs. 1 and 2 we present the results obtained for the pressure tensor and the temperature for a 50.99, as a function of the scaled coordinate s. The FIG. 1. Profiles of the diagonal components of the pressure tensor Pii in the steady state for a 50.99. The crosses correspond to Pxx and the circles to Pyy . The solid line is the hydrodynamic pressure p. Distance is measured on the scale sdefined in the text. The lowest temperature in the system is Tmin.0.57TW, while the highest density is nmax.1.15n(0). FIG. 2. Slope of the temperature profile as a function of the position for the same state as in Fig. 1. Quantities are measured in the reduced units defined in the main text. 2022 57 J. J. BREY AND D. CUBERO results have been averaged over 100 trajectories. The origin for x~and for s) has been taken at the same distance from both walls. Therefore, the system is symmetric around s 50. It is seen that, outside the boundary layers, the pressure is uniform and the temperature profile is accurately described by a parabola, as predicted by Eq. ~26!. The solid line in Fig. 2 is the numerical fit used to determine the value of ] 2T(s)/ ] s2in the bulk. Nevertheless, Pxx (s)is clearly larger than Pyy (s), contrary to the prediction of the Navier-Stokes approximation ~23!. Although the boundary layers can be studied in detail by applying the same methods as developed for molecular gases, they will not be discussed here. As a quantitative test, we have plotted in Fig. 3 the simulation values of ] 2T(s)/ ] s2and also the theoretical expression, i.e., the right-hand side of Eq. ~26!, as a function of g ( a ), for a >0.99. This will be referred to as the quasielastic region in the following. The simulation values have been obtained by fitting the numerical data for ] T(s)/ ] sto a straight line as indicated in Fig. 2. It is seen that the agreement is excellent. Nevertheless, let us note that over the range of values of g shown, the g dependence of the factor a( g ) appearing in Eq. ~26!is negligible and the curves obtained with a( g ) and a(0) are indistinghishable over the scale of the figure. In other words, the right-hand side of Eq. ~26!can be accurately approximated by a linear function of g in the quasielastic region. What happens when the value of a is decreased below the quasielastic region? The anisotropy of the diagonal terms of the pressure tensor increases as expected, but, in addition, the hydrodynamic pressure becomes nonuniform in the bulk of the system. The nonuniformity comes from the Pyy (s)component, while the Pxx (s)component remains homogeneous, as required by the exact balance of momentum ~17!. As an example, we present in Fig. 4 the components of the pressure tensor as functions of the scaled position sfor a 50.95. The slope of the temperature for the same situation is shown in Fig. 5. There is a small but perceptible curvature, indicating that the Navier-Stokes approximation is no longer valid for the temperature either. Nevertheless, let us mention that if the slope of the temperature in the bulk is approximated by a straight line, the resulting numerical value for the second derivative of the temperature is close to the prediction of Eq. ~26!. This is easily understood since the temperature gradient is still very small in the bulk of the system. We have investigated the behavior of the system up to a 50.8 and in all cases a steady state similar to that for a 50.95 is reached. The main effect of decreasing the value of a is to increase the anisotropy of the pressure tensor and also the influence of higher-order terms in the temperature profile. In the simulations we have also computed the one-particle distribution function of the system. The results will be discussed in Sec. IV. In the above figures we have presented the hydrodynamic profiles as functions of the scaled variable s. This is a convenient representation in order to carry out a comparison with the predictions of our model kinetic equation. Of course, it is possible to relate sto the spatial coordinate xby means of the definition of the former, Eq. ~15!. This can be done using the hydrodynamic profiles obtained from the simulation and evaluating the integral numerically. Alternatively, we can use the theoretical profiles provided by Eqs. ~23!and ~26!. In this latter case we only need the values of the temperature at the origin T(0) and the pressure in the bulk p. In Fig. 6 we have plotted sas a function of xfor the system with a 50.99. It is seen that the predictions of the theory match very well the numerical data over all the size the system. To avoid misunderstandings, let us note that the variable sis measured from the center of the system and therefore the influence of the boundary layers is reduced. Due the complexity of the hydrodynamic profiles in the steady state, we have not been able to construct an exact FIG. 3. Second derivative of the temperature with respect to the scaled position sas a function of the dissipation parameter g 51 2 a 2. The crosses are values obtained from simulation by fitting the bulk data as in Fig. 2 and the solid line is the prediction of the Navier-Stokes approximation ~26!. Quantities are measured in the reduced units defined in the text. FIG. 4. Same as Fig. 1 for a 50.95. Now it is Tmin.0.14TWand nmax.1.18n(0). FIG. 5. Same as Fig. 2 for a 50.95. 57 2023STEADY STATE OF A FLUIDIZED GRANULAR MEDIUM . . . solution to the model kinetic equation describing the state for arbitrary values of the restitution coefficient. Therefore, in order to try to understand the origin of the rheological effects observed in the simulation, it seems worthwhile to look for a perturbative solution. Given that in the steady state we are considering the gradients of the hydrodynamic fields are induced by the inelasticity in collisions, an expansion in the latter seems the more appropriate one. This will be the subject of the next section. IV. EXPANSION IN THE DISSIPATION PARAMETER The distribution function of the steady state is given by the solution of Eq. ~16!, which can be we written in the equivalent form vx ] ] sf52f1fl2 e 2 v dnkB n ] ] Tfl,~29! where e [12 a 2is the degree of inelasticity. In addition, in order to simplify the notation we have omitted the superscript sand here and in the following ] / ] T( ] / ] n) is understood to be taken at constant density ~temperature!. We are going to use a modified Chapman-Enskog expansion method to get a solution of Eq. ~29!in the form of a series expansion around e 50. The physical reason is that we are interested in a system whose boundaries are kept at constant temperature. In such a situation, we know that, neglecting boundary effects, the only steady state for an elastic gas corresponds to ~uniform!equilibrium given by the MaxwellBoltzmann distribution. Therefore, the gradients present in the steady state of a granular fluid must be functions of e vanishing in the limit e →0. Then we formally expand f5f01 e 1/2f11 e f21 e 3/2f31••• ~30! and ] ] s5 e 1/2 ] 11 e ] 21 e 3/2 ] 31•••.~31! In particular, for the temperature and the density we consider ] T ] s5 e 1/2 ] 1T1 e ] 2T1 e 3/2 ] 3T1•••,~32! ] n ] s5 e 1/2 ] 1n1 e ] 2n1 e 3/2 ] 3n1•••.~33! The expansion in powers of e 1/2 is motivated by the results obtained in the Navier-Stokes approximation discussed in Sec. II, which are expected to be valid in the asymptotic limit e →0. There it was found @see Eq. ~26!# that ] 2T/ ] s2 ; e , i.e., ] T/ ] s; e 1/2. Moreover, the pressure was uniform, so that the leading contribution to ] n/ ] sis also of the same order e 1/2. Then, from Eq. ~29!it follows that a balance to order e 1/2 is possible only if fhas contributions of that order. In any case, the verification of the consistency conditions to be discussed below will indicate whether the expansions are correct, at least up to the order considered in the calculations. Of course, a different question is the convergence of the expansion. We will return to this point at the end of this section. When the expansion in Eq. ~30!is introduced into the expression of v , Eq. ~6!, one gets v @f u f#5 v 01 e 1/2 v 11 e v 21 e 3/2 v 31•••,~34! with v 05 v @f0 u f0#, v 152 v @f0 u f1#, v 25 v @f1 u f1#12 v @f0 u f2#,.... ~35! Substitution of Eqs. ~30!,~31!, and ~34!into Eq. ~29!and equating coefficients of the same power of e 1/2 leads to the equations f05fl,~36! f152vx ] 1fl,~37! and fr52vx( q51 r ] qfr2q22 v r22 dnkB n ] ] Tfl~38! for r>2. From Eq. ~36!it follows that to zeroth order the distribution function is the local equilibrium one and the corresponding contributions to the pressure tensor and heat flux defined by Eqs. ~13!and ~14!are P0,ij5p d ij,q050.~39! Since the local equilibrium distribution reproduces, by definition, the exact values of the hydrodynamic fields, we have the consistency conditions E dvfr5 E dvvfr 5 E dv v 2fr 50~40! for r>1. When applied to f1, given by Eq. ~37!, these conditions require n ] 1T1T ] 1n50, ~41! i.e., ] 1p50. Then one gets f1~s,v!52vx ] 1T S ] ] T21 T D fl~42! and P1,ij50, q1,i5 d ix ~d12!nkB 2T 2m ] 1T.~43! FIG. 6. Scaled variable sas a function of the spatial coordinate x~both in reduced units!for a 50.99. The symbols are simulation results and the solid line is the prediction of the theory. 2024 57 J. J. BREY AND D. CUBERO To obtain f2we need v 0defined in Eq. ~35!. A simple calculation gives v 05d 2 g ˜ p n ,~44! where we have introduced g ˜ [ g / e , with g given by Eq. ~22!. Then, from Eq. ~38!with r52 we get f2~s,v!52vx ] 2T ] ] Tfl2vx~ ] 2n!fl n1vx 2~ ] 1 2T! S ] ] T21 T D fl 1vx 2~ ] 1T!2 S ] 2 ] T222 T ] ] T12 T2 D fl2 g ˜ T ] ] Tfl. ~45! The consistency conditions to this order read ] 2p50, ] 1 2T5 g ˜ md ~d12!kB.~46! In fact, the above equation for the temperature gradient can also be obtained by substituting the expression of the heat flux to first order ~43!into the energy conservation law ~18!and by restricting ourselves to first order in e . A similar test can be carried out to each perturbation order to verify the self-consistency of the calculations. The results for the pressure tensor and heat flux contributions are P2,ij5 d ijp S kB m ] 1 2T2 g ˜ D 1 d ij d ix 2kBp m ] 1 2T 52 d ijp2 g ˜ ~d12!~12 d ixd!,~47! q2,i52 d ix ~d12!nkB 2T 2m ] 2T.~48! Upon writing Eq. ~47!we have used Eqs. ~46!. Therefore, anisotropy of the diagonal terms of the pressure tensor appears to order e . Note that it is associated with a contribution of second order in the temperature gradient to the pressure tensor. Hence a description at the level of the Navier-Stokes approximation cannot predict the normal stress differences, in agreement with the results obtained in Sec. II. Combination of Eqs. ~39!,~43!,~47!, and ~48!provides expressions for the pressure tensor and the heat flux, valid up to first order in e , Pij5 d ijp F 12~12 d ixd!2 g d12 G ,~49! qi52 d ix ~d12!nkB 2T 2m ] T ] s.~50! Therefore, the normal stress ratio in this approximation is given by Pxx Pyy 511 g 2d d12.~51! When these expressions for the fluxes are used in the conservation laws ~17!and ~18!, the following equations determining the hydrodynamic profiles are obtained: ] p ] s50, ] 2T ] s25 g md ~d12!kB.~52! Since the theory still predicts uniform pressure at this order, it is clear that the range of applicability of the above results is restricted to values of e lying in the quasielastic region as characterized in Sec. III, i.e., a *0.99. The equation for the temperature profile coincides with the linearization in g of the one derived in the Navier-Stokes approximation ~26!. Therefore, we already know from the discussion in Sec. III that it accurately fits the simulation results in the limit of quasielasticity. In Fig. 7 the ratio Pxx /Pyy is plotted as a function of g . The crosses correspond to the direct Monte Carlo simulation results and have been obtained by space averaging in the bulk. The values of the simulation parameters are the same as in the previous figures. The continuous line is the theoretical prediction given by Eq. ~51!.It is seen that the theory describes qualitatively well the asymmetry of the diagonal components of the pressure tensor, but there is a clear quantitative discrepancy. This is not surprising since g is proportional to the inverse of the constant C appearing in the expression of the effective collision frequency, Eq. ~5!. We have fixed its value by requiring the model to reproduce the Boltzmann heat conductivity in the absence of dissipation and then taking C.1. If we would have chosen Cto give the same value as the Boltzmann equation for the shear viscosity, it should be C.2 and we would have a much better agreement for the components of the pressure tensor. Nevertheless, it is clear from Eq. ~52! that the agreement for the temperature profile should be worse in this case. This is a well known limitation of single relaxation models of the Boltzmann equation. They cannot reproduce simultaneously the correct shear viscosity and the correct heat conductivity. Let us now consider the distribution function. Figure 8 depicts the marginal velocity distribution for vy, the component of the velocity parallel to the thermal walls. The coefFIG. 7. Ratio of normal stresses Pxx /Pyy as a function of the dissipation parameter g 512 a 2. The solid line corresponds to the present theory and the crosses are results from the Monte Carlo simulation. The dotted line is a linear fit of the data. 57 2025STEADY STATE OF A FLUIDIZED GRANULAR MEDIUM . . . ficient of restitution is a 50.99 and the results shown correspond to the layer located at s524.59 (x5212.62). This position has been chosen such that it is as far as possible from the center of the system, but still inside the bulk region where the pressure is constant ~see Fig. 1!. Similar results are obtained at other positions in the bulk. The solid line corresponds to the integration with respect vxof the solution of the model up to order e ,f(2)5fl1 e 1/2f11 e f2. It is seen that the agreement is excellent in the velocity range shown. Of course, the distributions is symmetric around vy50, as required by the own symmetry of the system. On the other hand, the marginal distribution for the velocity vxperpendicular to the thermal walls is asymmetric ~see Fig. 9!. This reflects that particles coming from the nearest wall have more energy that those coming from the most distant one. The asymmetry of the distribution function can be more clearly appreciated by plotting the ratio between the marginal distribution and the corresponding local equilibrium distribution fl,x(vx), as shown in Fig. 10. There it is also seen that the agreement between theory and simulation becomes much worse when large velocities are considered. This is to be expected since in our model kinetic equation the evolution of the one-particle distribution function is governed by only the first five velocity moments. Similar results are found in other cases, although the asymmetry increases as the distance from the center of the system increases and also as the coefficient of restitution decreases. Given the accuracy of our model to describe the behavior of the system in the quasielastic limit, it is tempting to try to extend the above calculations to the case of stronger dissipation. In particular, it should be interesting to explain the origin of the inhomogeneity of the pressure. Therefore, we have evaluated the thirdand fourth-order contributions in the e 1/2 expansion of the distribution function. The derivation is straightforward but lengthy and tedious, and we only quote the results here. Some details of the calculations are given in the Appendix. The final expressions for the pressure tensor and the heat flux, valid up to order e 2, are Pij5 d ijp H 12~12 d ixd! F 2 g d1213 g 2~24d21d12! 2~d12!3 G 2~12 d ixd!3 g ~1514d!kB 16d~d12!mT S ] T ] s D 2 J ,~53! qi52 d ix nkB 2T m F d12 214 g ~2d218d21! d12 G ] T ] s.~54! Equations for the pressure and the temperature follow by substituting the above expressions into Eqs. ~17!and ~18!.Of course, terms of order higher than g 5/2 have to be consistently neglected. The result is ] p ] s52p2b~ g !~d21!kB dm ] ] s F 1 T S ] T ] s D 2 G ,~55! ] 2T ] s25c~ g !1b~ g !1 T S ] T ] s D 2 ,~56! where b~ g !5 g 3~1514d! 32~d12!,~57! c~ g !5 F g d d122 g 264d31253d2238d 4~d12!3 G m kB.~58! FIG. 8. Marginal distributions for the velocity component vy perpendicular to the gradients. The symbols are simulation data and the solid line results of the present theory. The coefficient of restitution is a 50.99 and we have considered the gas layer located at x5212.62. Quantities are measured in the reduced units defined in the text. FIG. 9. Marginal velocity distribution in the direction of the gradients. The system and location of the layer considered are the same as in Fig. 8. FIG. 10. Same as Fig. 9, but the marginal distribution is scaled with the Maxwellian defined by the local density and temperature. 2026 57 J. J. BREY AND D. CUBERO A comparison of the terms of order g with those of order g 2, for instance, in the expression of c( g ), shows that the numerical coefficients of the latter are much larger than those of the former. This clearly indicates that the expansion we have carried out, although probably asymptotic, is divergent. In fact, if one compares the predictions of Eqs. ~55!and ~56! with the simulation results for a 50.99, worse agreement than by using the e approximation is obtained. Similar behavior has been found in the usual Chapman-Enskog expansion applied to a molecular fluid under uniform shear flow @18#. Equation ~55!deserves several comments. A simple estimation of the terms appearing on the right-hand side leads to ] p/ ] s}s, with a negative constant of proportionality. Therefore, it predicts a pressure profile with a curvature that has opposite sign to the one observed in the simulations ~see Fig. 4!. Although the origin of this strong discrepancy is not clear, we think that it is due to the divergent character of the e 1/2 expansion and not to the inaccuracy of the model kinetic equation. This is a point that deserves future work. V. CONCLUSIONS The objective of this work has been to study an inhomogeneous steady state of a granular fluid with homogeneous boundary conditions. The existence of such a state is a consequence of the inelastic character of collisions and has no analog in molecular fluids. By means of direct Monte Carlo simulations, it has been shown that the system exhibits anisotropy of the diagonal elements of the pressure tensor. This rather unexpected effect for a system without velocity flow appears even in the low dissipation limit, in which the pressure is uniform in the bulk, i.e., outside the boundary layers. The observation of the normal stress differences and the identification of its origin is one of the central physical results presented here. We have also studied the velocity distribution function, which is non-Gaussian and asymmetric, presenting an overpopulation of particles coming from the nearest wall with large velocities. In the quasielastic region, roughly defined by a coefficient of restitution larger than 0.99, the simulation results are accurately described, both qualitatively and quantitatively, by a simple model kinetic equation proposed recently @6,10#. This refers not only to the hydrodynamic fluxes and fields, but also to the more detailed information provided by the velocity distribution, at least in the region of thermal velocities. Nevertheless, this requires going beyond the Navier-Stokes approximation. In fact, the peculiarity of the state under consideration makes more appropriate an expansion in powers of the inelasticity parameter than the usual Chapman-Enskog expansion in the gradients of the hydrodynamic fields. This is a consequence of the strong coupling between gradients and dissipation or, more precisely, of the fact that gradients are induced by the inelastic character of collisions. In the limit of zero inelasticity the equilibrium solution of the Boltzmann equation is recovered. This renders a perturbative approach in powers of e 1/2 possible. However, it must be pointed out that the usefulness of this kind of expansion can be limited since our results strongly suggest that they are asymptotic but divergent. This does not mean that the model kinetic equation is not relevant for the description of low-density granular fluids with strong inelasticity, but that we must look for solutions that are valid in the nonperturbative regime. In fact, this procedure has proved to be very fruitful for molecular gases @23#. Finally, let us stress that although our study has been restricted to a low density granular fluid, there is no reason to expect that the reported effects are negligible for highdensity granular fluids. In fact, the presence of highly inhomogeneous spatial distributions already has been noted elsewhere @12#and we believe that the asymmetry of the normal components of the pressure tensor is also present and can be analyzed both theoretically, by means of the extension of the present model to the revised Enskog kinetic theory @6#, and by using molecular-dynamics simulation. ACKNOWLEDGMENTS The authors want to thank M.J. Ruiz Montero for stimulating discussions and a critical reading of the manuscript. This research was partially supported by Grant No. PB960534 from the Direccio ´n General de Investigacio ´n Cientı ´fica yTe ´ cnica ~Spain!. APPENDIX: THIRDAND FOURTH-ORDER CONTRIBUTIONS The expression for f3contains v 1. Using Eq. ~42!in the expression of v 1given by Eq. ~35!we have v 152 m s d21 p ~d21!/2 4G S d13 2 D ] 1T E dv1 E dv2 u v12v2 u 3 3fl~v1!v2x S ] ] T21 T D fl~v2!50~A1! since the integrand is odd with respect to the change v1,v2 →2v1,2v2. This is an expected result since v is a scalar that cannot couple linearly to f1, which is a vector. Then, from Eq. ~38!we get f3~s,v!52vx ] 3fl2vx ] 2f12vx ] 1f2 52vx ] 3fl1vx 2~ ] 2 ] 11 ] 1 ] 2!fl 2vx 3 ] 1 3fl1 g ˜ vx ] 1 S T ] ] Tfl D .~A2! Although the above expression can be written in a more explicit form, we have found it more suitable for calculations. Let us consider the consistency conditions ~40!. For the first one we have E dvf3~s,v!5 ] 3 E dvvxfl1~ ] 2 ] 11 ] 1 ] 2! E dvvx 2fl 2 ] 1 3 E dvvx 3fl1 g ˜ ] 1 S T ] ] T E dvfl D 5~ ] 2 ] 11 ] 1 ] 2!nkBT m,~A3! which vanishes identically since we have previously obtained ] 1p5 ] 2p50@see Eqs. ~41!and ~46!#. Similarly, we get 57 2027STEADY STATE OF A FLUIDIZED GRANULAR MEDIUM . . .