# Shear state of freely evolving granular gases

###### Abstract

Hydrodynamic equations are used to identify the final state reached by a freely evolving granular gas above but close to its shear instability. The theory predicts the formation of a two bands shear state with a steady density profile. There is a modulation between temperature and density profiles as a consequence of the energy balance, the density fluctuations remaining small, without producing clustering. Moreover, the time dependence of the velocity field can be scaled out with the squared root of the average temperature of the system. The latter follows the Haff law, but with an effective cooling rate that is smaller than that of the free homogeneous state. The theoretical predictions are compared with numerical results for inelastic hard disks obtained by using the direct Monte Carlo simulation method, and a good agreement is obtained for low inelasticity.

###### pacs:

45.70.Mg, 45.70.Qj, 47.20.Ky## I Introduction

Granular gases are ensembles of macroscopic particles whose dynamics is controlled by inelastic binary collisions, separated by ballistic motion Ca90 . The forces are usually short range and repulsive and, therefore, in the simplest form granular gases are modeled as a collection of smooth inelastic hard spheres or disks. In recent years, the study of granular gases has attracted a lot of attention ByP04 ; Du00 , among other reasons because they serve as a sensitive proving ground for kinetic theory and non-equilibrium statistical mechanics.

A peculiarity of granular gases is their tendency to spontaneously develop pattern forming instabilities when evolving freely. This includes the so-called shearing and clustering instabilities GyZ93 ; McyY94 , in which the system develops patterns both in the flow field and the density field. Analysis of the hydrodynamic Navier-Stokes equations for granular gases indicates the existence of several possible scenarios explaining the development of these instabilities and their mutual influence. An illuminating discussion of them is given in ref. PAFyM08 . For instance, the existence of non-stationary channel flows in freely cooling gases leading to an attempted finite-time density blowup, which is arrested by heat diffusion, has been shown recently ELyM05 ; MFyV08 .

In this paper, a shear flow state of freely evolving dilute granular gases will be described. Its origin is tied in with the shear mode instability and the analysis to be presented here is restricted to systems whose size is slightly larger than the critical size for this instability. The latter happens to be smaller than the critical size for the instability of one of the longitudinal modes associated with the density BDKyS98 ; WByE02 . Consequently, the spatial inhomogeneities of the hydrodynamic fields are expected to be small. Although the flow is non-stationary, all the time dependence can be expressed through the spatially averaged temperature and, therefore, it can be described in terms of dimensionless time-independent quantities. More precisely, it is an inhomogeneous two band shear state. It must be stressed that this state is exhibited by a system with, for instance, periodic or elastic boundary conditions, without being driven by the boundaries as it happens with the so-called simple or uniform shear flow Lu84 ; JyR88 ; SGyN96 . Non-stationary states in which all the time dependence occurs through the temperature, like the homogenous cooling state Ha83 , as well as non-uniform states where all the spatial dependence also occurs through the temperature are peculiar of granular gases BCMyR01 .

The existence of this free shear state has been suggested previously in the framework of a time-dependent Ginzburg-Landau model for granular gases WByE02 , and also as the result of a nonlinear analysis of the shearing instability SMyM00 . The theory presented in this paper differs from those studies both in the method and the results, where rather strong discrepancies occur as it will be discussed along the paper. The interest here is focussed on the complete identification of the hydrodynamic fields characterizing the free shear state and the comparison of the theoretical predictions with the results of numerical simulations of the particles composing the system. Of course, this provides a very demanding test of the validity of the macroscopic description provided by the hydrodynamic Navier-Stokes equations for granular gases.

The remainder of the paper is organized as follows. In Sec. II, the nonlinear Navier-Stokes hydrodynamic equations for a dilute granular gas are shortly reviewed, and particularized for the free shear state. This is macroscopically characterized by presenting gradients only in one direction and by a velocity field perpendicular to the gradients. Moreover, it is assumed that the density inhomogeneities are small. By introducing appropriate dimensionless length and time scales, a closed equation for the velocity flow is derived in Sec. III. The solution of this equation shows that the amplitude of the velocity field decays monotonically in time. Moreover, by analyzing the temperature equation, an explicit expression for the amplitude of the steady density fluctuation is obtained. Also, it is seen that the ratio between the average temperature of the system and the square of the amplitude of the velocity field is time independent.

In Sec. IV, the simulation method used is indicated. It is restricted to low density gases, whose one-particle distribution function obeys the Boltzamnn equation, but it must be stressed that it is an -particle simulation method, therefore providing information on all the space and time correlations. It is observed that the system actually tends to a steady situation having the assumed properties of the free shear state. Then, the measured hydrodynamic quantities are compared with the analytical theoretical expressions derived from the Navier-Stokes equations. A good agreement is obtained for small inelasticity and system sizes close to the critical value for the shear instability. Finally, Sec. V summarizes the obtained results and presents a brief discussion.

## Ii Hydrodynamic equations

The balance equations for the number density , the flow velocity , and the granular temperature of a system composed by smooth inelastic hard spheres () or disks () of mass , diameter , and coefficient of normal restitution , have the form

(1) |

(2) |

(3) |

For a dilute gas and to Navier-Stokes order, the pressure tensor , the heat flux , and the cooling rate are given by the constitutive relations BDKyS98 ; ByC01 ; Go03

(4) |

(5) |

(6) |

Here is the hydrostatic pressure, is the shear viscosity, the (thermal) heat conductivity, and a transport coefficient peculiar of granular fluids that is referred to as the diffusive heat conductivity. The term denotes the contribution to the cooling rate of zeroth order in the gradients of the hydrodynamic fields, while denotes the second order in the gradients contributions. The latter will be neglected in the following, as it is done in most of the calculations, since it is expected to give very small corrections to the hydrodynamic equations as compared with the similar terms coming from the pressure tensor and the heat flux BDKyS98 . The transport coefficients appearing in Eqs. (4) and (5) can be expressed as

(7) |

(8) |

(9) |

where and are the low density (Boltzmann) values of the shear viscosity and the thermal heat conductivity, respectively. Their explicit expressions as well as those of the dimensionless functions , , and are given in Appendix A. Finally, the expression of the zeroth order cooling rate reads

(10) |

The function is also given in Appendix A.

The hydrodynamic Navier-Stokes equations for the inelastic gas are obtained by substituting Eqs. (4)-(6) into Eqs. (1)-(3). They admit a simple solution defined by a constant and uniform density , a vanishing velocity field , and a time dependent temperature obeying the equation

(11) |

This is the so-called homogeneous cooling state (HCS) Ha83 . This state is known to be unstable with respect to long wavelength spatial perturbations. More precisely, linear stability analysis of the hydrodynamic equations shows that, for large enough systems, the transversal component of the velocity relative to the square root of the temperature grows in time leading to the instability of the system GyZ93 ; McyY94 ; BDKyS98 . This shear instability has been confirmed by molecular dynamics simulations GyZ93 ; McyY94 and also by Monte Carlo simulations of the effective dynamics associated to the Boltzmann equation BRyC96 . Of course, after a short time interval, the theoretical predictions based on the linearized hydrodynamic equations are not valid any longer since neglected nonlinear effects become very important. During the initial stages of the development of the instability, density inhomogeneities occur in the system. For dilute systems whose linear size is larger (and comparable) to the critical system size for the shear instability, the spatial inhomogeneities seem to be dominated by a nonlinear coupling of the density with the transversal velocity field BRyC99 . Although this could indicate that the system is going into a clustering regime, in which the particles tend to group together forming very high density regions coexisting with very dilute regions, simulation results indicate that the density actually tends to a quite smooth profile BRyC99 .

Because of continuous cooling due to the inelasticity of collisions, stationary states are not possible in freely evolving granular gases (for instance, with periodic boundary conditions). Nevertheless, it is still possible that the “final state” reached by the system when the HCS is unstable, exhibits some scaling properties so that it can be simply identified at least at the hydrodynamic level of description. Here, attention will be focused on systems just above the shear mode instability threshold. Well above this threshold, the physical scenario might be quite different ELyM05 ; MFyV08 ; PAFyM08 .

Consider hydrodynamic flows verifying the following conditions: i) there are gradients in only one direction taken as the axis, and ii) the hydrodynamic velocity field has the form . Other specifications of the state of the system will be made when appropriate. Use of the above conditions into Eqs. (1)-(3) shows that the density profile must be stationary, , and the component of the pressure tensor must be uniform. In the framework of the Navier-Stokes approximation, the latter implies that the pressure is also uniform, . Moreover, the equations for the velocity and temperature fields become

(12) |

(13) |

To solve the above equations, the boundary conditions must be specified. A system of particles enclosed in a cubic () or square () box of side will be considered. Periodic boundary conditions will be assumed at all the boundaries in order to avoid undesirable wall effects. It will be convenient for the purposes here to introduce a function by

(14) |

with . The conservation of the number of particles and the periodic boundary conditions imply that

(15) |

Moreover, it will be assumed that , i.e. the spatial density inhomogeneities are supposed to be small. It will be shown in the following that this actually reduces the range of applicability of the theory to systems near (and above) the threshold of the shear instability. Then, the temperature profile in the state being considered is given by

(16) |

where

(17) |

is the spatial average temperature of the system at time . It must be realized that assuming and, consistently, retaining terms up to first order in it when solving Eqs. (12) and (13), is not the same as linearizing around the HCS, since nothing is in principle assumed about the amplitude of the temperature or the velocity field . This will be discussed in detail in the next section.

## Iii The free shear state

The hydrodynamic equations (12) and (13) can be simplified by using dimensionless length and time scales defined by

(18) |

and

(19) |

respectively. Here

(20) |

is a characteristic frequency, proportional to the Boltzmann collision frequency , namely . Note that the pre-factor in Eq. (18) does not depend on time and the length scaling is made with a characteristic length which is proportional to the (time-independent) average mean free path of the system. In terms of the new scales, Eq. (12) can be approximated by

(21) |

where terms of order have been neglected. The solution of the above equation is a superposition of monocromatic waves of the form

(22) |

with

(23) |

and

(24) |

where and are arbitrary constants. The possible values of are restricted by the periodic boundary conditions which imply

(25) |

being a positive integer and the value of for , i.e.

(26) |

The form of Eq. (24) indicates that in the limit of large time , the dynamics of the system is governed by the fundamental mode corresponding to the lowest possible value of , i.e. . Therefore, for large enough times,

(27) |

with and . In the same approximation under consideration, , the evolution equation for the temperature (13) leads to

(28) |

Using Eq. (23) this becomes

(29) |

By requiring the sum of the position dependent terms to cancel it is obtained that

(30) |

where

(31) |

and

(32) |

Consistency requires that be actually independent of and, moreover, that it be defined positive. The former of these conditions together with Eq. (21) leads to

(33) |

This equation will be used later on to determine the effective cooling rate of the average temperature of the system in the free shear state. On the other hand, equating to zero the part of Eq. (III) that is position independent yields, after employing Eq. (33),

(34) |

As mentioned above this quantity must be positive. As a consequence, the existence of the mathematical solution of the hydrodynamic Navier-Stokes equations under consideration and hence of the free shear state is only posible if

(35) |

or, equivalently, , with the critical size given by

(36) |

This coincides with the condition determining the instability region of the shear mode found in the linear stability analysis of the hydrodynamic equations around the HCS BDKyS98 .

Substitution of Eq. (34) into Eq. (30) provides the explicit form of the steady density profile in the free shear state near the threshold of the instability,

(37) |

(38) |

Therefore, the condition formally implies that . How restrictive this condition actually is will be seen in the next section.

A particularly simple expression for the amplitude of the velocity field is obtained by scaling with the average temperature. Equations (32) and (34) yield

(39) |

Therefore, the scaled macroscopic velocity field is time-independent and all its dependence on the average density and inelasticity occurs through the critical length . A similar expression can be obtained for the ratio between the effective cooling rate of the shear state and the cooling rate of the HCS, . The former is identified by writing Eq. (33) in the form

(40) |

with

(41) |

that leads to

(42) |

Of course, in the limit , and the law for the HCS given by Eq. (11) is recovered.

The average energy per particle is

(43) |

with the bar denoting spatial average. Using Eqs. (17), (22), (23), (32), and (34) it is found:

(44) |

For , the equilibrium relation is recovered as expected. This equation has been previously obtained by Wakou et al. WByE02 , in the context of a Landau-Ginzburg-type equation of motion derived from the hydrodynamic equations, under certain restrictions and in the limit of nearly elastic collisions (). At this point, it is worth mentioning that the assumed periodic boundary conditions do not play an essential role in the theory developed here, although they must be compatible with the symmetry of the shear state. For instance, it is easily realized that completely equivalent results are obtained if elastic walls were used at the boundaries of the system.

Soto et al., SMyM00 have carried out a nonlinear analysis of the hydrodynamic equations of a system of inelastic hard disks close to the instability threshold. When is slightly smaller than , vorticity modes with wavenumber exhibit a critical slowing down, so that all the other hydrodynamic modes can be considered as enslaved by them. Then, it is possible to derive closed equations for the amplitudes of the vorticity modes with by using the adiabatic elimination method. This approach can be expected to be formally consistent with the one presented here, but quantitative and qualitative discrepancies occur when comparing the results from both approaches. Although the authors of SMyM00 do not give almost any detail of their analysis of the nonlinear hydrodynamic equations, we have identified a twofold origin of the discrepancies. Firstly, a nonlinear term involving the vorticity seems to have been inconsistently neglected. Secondly, the linearization around the HCS and the adiabatic method used in SMyM00 is not equivalent to the approximation scheme presented here, in which the velocity field obeys the closed Eq. (21).

## Iv Direct Monte Carlo simulations

To test the theoretical predictions presented in the previous sections and, in particular, the existence itself of the free shear state, we have employed the direct simulation Monte Carlo (DSMC) method Bi94 ; Ga00 , to numerically generate the dynamics of a system of inelastic hard disks. The DSMC method is a many-particle algorithm designed to mimic the effective dynamics of the particles of a gas in the low density limit. Therefore, in this limit it is expected to lead to the same results as, for instance, molecular dynamics simulations, with the advantage of a much larger statistical accuracy. In all the simulations, the initial state was homogeneous and isotropic with a Gaussian velocity distribution, and a square box with periodic boundary conditions was used. The linear size of the system was always larger than, but close to, the critical one , predicted by Eq. (36), for the considered value of the restitution coefficient . Then, according to the theory, the system is expected to generate the non-linear free shear flow described in Sec. III, if it is stable.

As usual in DSMC simulations, the mass of the particles will be taken as the unit of mass, the average mean free path as the unit of length, and , where is the initial temperature, as the unit of energy. The effect of a collision between particles and is to instantaneously modify their velocities according to the rule

(45) |

(46) |

where is the relative velocity and is the unit vector pointing from the center of particle to the center of particle at contact. This corresponds to the scenario in which the constitutive relations (4)-(6) were derived BDKyS98 ; ByC01 .

One of the technical difficulties when numerically simulating a freely evolving granular fluid is the continuous cooling of the system, i.e. the decrease in magnitude of the typical velocity of the particles. As a consequence, the numerical inaccuracies become very large after some time interval. To avoid this difficulty, a procedure was introduced based on a change in the time scale being used. More specifically, a logarithmic time scale is introduced Lu01 ; BRyM04 . Then, the dynamics of systems in states where all the time dependence comes through the average temperature because of inelastic cooling, can be easily mapped into the dynamics around steady states. In this case, the logarithmic time scale is proportional to the cumulated number of collisions per particle.

### iv.1 Transient dynamics

In most of the simulations to be reported, a similar time sequence was observed. The system remains homogeneous with no macroscopic velocity field for an initial period of time, developing afterwards a state with two vortices, and finally a shear state characterized by two counterflows parallel to one of the sides of the system. In this state, there is a coupling between the velocity and density fields. The spatial inhomogeneities remain stationary and the state looks time independent in the logarithmic time scale used in the numerical simulations. The particular location of the transient vortices and their direction, and hence of the bands in the shear state, depends on the initial state or, in a statistical sense, on the fluctuation taking the system away from the HCS. An example of the described behavior is given in Fig. 1, where the scaled velocity field is plotted at four different times for a system with and (. The formation of the vortices, their distortion, and the formation of shear bands are clearly identified. The latter remain unchanged after their formation, and this was observed in all cases, at least as long as the system is close enough to the shear instability. The spontaneous symmetry breaking occurring in the final shear state was observed in both perpendicular directions and, to compare with the theoretical predictions, the -axis was always chosen perpendicular to the velocity field. More will be commented on this issue in the final section of the paper.

To get some additional insight about the time evolution of the system before reaching the free shear state, in Fig. 2 the ratio between the spatial average temperature and the HCS temperature at the same time as predicted by the Haff law, Eq. (11), is shown as a function of the average number of collisions per particle . The system is the same as in Fig. 1 and three different simulation trajectories are plotted. Although the time evolution of the temperature is not exactly the same in all cases, a quite similar trend is observed. There is a time interval in which the average temperature is well described by the Haff law () in spite of the system having already well developed vortices (see Fig. 1). Afterwards, for , the average temperature decays much slower than in the HCS, so that the ratio grows very fast, until it saturates at a given value. The region with the largest growth corresponds to the distortion of the velocity vortices. The constant value reached in the final shear state indicates that the average temperature still obeys a Haff-like law, but with a different cooling rate, in qualitative agreement with Eq. (40).

The evolution of the density is illustrated in Fig. 3, also for the system with and . Two Fourier components of the relative density field are plotted, and . The former corresponds to the second density mode that is the one predicted to survive by the theory, Eq. (37). In the simulations, it appears parallel to any of the sides of the system as discussed above and no distinction is made when reporting the simulation results. The other component, is the lowest mode along the diagonal of the system. Again, the several curves are different simulation trajectories. In all cases, the second transversal mode of the density grows from the noise level to a constant value as the shear instability develops. Moreover, the final steady value is the same for all trajectories. On the other hand, although the diagonal mode often also grows in the time interval in which the two vortices are losing their shape transforming into the shear state, it decays to the noise level once this state is reached.

### iv.2 Shear state results

Now the results obtained for the properties of the system once in the free shear state will be reported, and compared with the theoretical predictions obtained in this paper. Consider first the density profile. The simulation data are very well fitted by a cosine function as given in Eq. (37). This has been checked both by plotting directly the profiles and by computing their Fourier components. The amplitude of the measured perturbation is plotted in Fig. 4 as a function of for three values of the restitution coefficient, , and . The lines are the theoretical predictions given by Eq. (38), the highest one corresponding to the smallest value of and the lowest one to the greatest value of . It is seen that the agreement is quite good when is close to , the discrepancies increasing as the size of the system increases above its critical value. Also, the deviation is larger the smaller the coefficient of restitution. In any case, it must be kept in mind that the theory developed here is by construction restricted to states with .

In Fig. 5, the amplitude of the steady cosine profile for is presented for the same values of the parameters as considered in Fig. 4. According to the theory developed here, Eq. (16), this amplitude should be the same as for the density profile, i.e. . Nevertheless, the simulation data indicate that this is not the case, except for values of close to unity. The deviations of the numerical data from the theoretical predictions are in opposite directions for and , being larger for the former. This indicates that the pressure is not actually strictly uniform as predicted by the Navier-Stokes equations for dilute granular gases, but exhibits some oscillatory profile. On the other hand, the component of the pressure tensor was found to be uniform, as required by the balance equations (1)-(3). Nevertheless, for small enough, the agreement between theory and simulation can be qualified as satisfactory. Let us mention that if the elastic values of the transport coefficients () were used, the theoretical prediction for as a function of would become independent of contrary to what is observed in the simulations.

A quite strong prediction of the theory is provided by Eq. (39), where the amplitude of the velocity field scaled with the square root of the average granular temperature is expressed as a simple time and -independent function of the ratio . This latter property was seen to ve verified in the simulation within the statistical uncertainties. Moreover, the results displayed in Fig. 6 show that the simulation data are in good agreement with the theoretical expression, although again systematic deviations are observed as the value of the coefficient of restitution decreases and/or the size of the system as compared with its critical value increases.

The simulation results also indicate that the decay of the average temperature in the free shear state is governed by a Haff-like law, Eq. (39). In Fig. 7, the ratio between the cooling rate of the free shear state, , and the one of the HCS, , is plotted, for the same systems being considered along this section. The theoretical prediction is provided by Eq. (42), implying that the free shear state cools slower than the associated HCS, i.e. with the same density and initial temperature. Again a quite satisfactory agreement is found between theory and simulations, with the discrepancies exhibiting the same trends as in all the previous comparisons.

## V Summary and discussion

It has been shown that the hydrodynamic Navier-Stokes equations for granular gases predict the existence of a shear state for freely evolving systems, whose size is slightly larger than the critical size for the shear mode instability. Explicit analytical expressions for the hydrodynamic fields have been derived. Although it is a non-stationary state due to cooling, all the time dependence of the hydrodynamic fields occurs through the average temperature and, therefore, it can be eliminated by introducing appropriate dimensionless quantities. These theoretical predictions have been found to be in good agreement with the numerical results obtained by the DSMC method, for values of the restitution coefficient close to unity. When the system is more inelastic, significant deviations from the theory are observed. This was expected, since the free shear state is characterized, as many other non-equilibrium states of granular gases, by a strong coupling between gradients and inelasticity. In the present case, this is easily identified through the dependence on of . As a consequence, a first order gradient expansion like the one leading to the Navier-Stokes equations also implies a limitation in the range of values of for which the theory applies. Moreover, it is probably true that the free shear state is inherently non-Newtonian, as it is the case for the steady uniform shear state of a granular gas SGyD04 . A possible indication of this is the inhomogeneity of the pressure observed in the simulations as implied by the difference between the measured values of and .

In ref. PAFyM08 several possible hydrodynamic scenarios are described for the behavior of a freely evolving granular fluid inside its instability region. In that classification, the situation considered in this paper is called scenario 4. Here the final state reached by the system has been investigated and a quantitative test for the complete scenario has been provided. The free shear state seems stable, in the sense that no deviations from it have been observed in the simulations. In this context, it is worth mentioning that, in many cases, the two-vortex state was clearly identified for a quite large period of time before the system moved to the shear state. This may indicate that it is a metastable state with a large escape time. This issue clearly deserves more attention. On the other hand, when becomes much larger than (typically, ), other hydrodynamic modes become relevant and the free shear state is not expected to occur in the system. This has been confirmed by the simulation results.

It is also worth comparing in some detail the approach followed here with the nonlinear stability analysis carried out in ref. SMyM00 . As already mentioned, the results reported there suggest that a term nonlinear in the velocity has been omitted although it is of the same order as others that are kept. This term dramatically modifies the results of the stability analysis DMGByR08 . Moreover, in SMyM00 an adiabatic elimination method is used, in which the time derivative of the hydrodynamic fields other than the transverse velocity are set equal to zero. In this way, these fields can be expressed in terms of , and when the expressions are inserted into the equation for , a closed equation is obtained for the latter. The approximation followed here is different. To lowest order, the transversal velocity field obeys a closed equation by itself, Eq.(21), once the function has been scaled out by the change of variables. This equation differs from the one obtained by the adiabatic elimination method. It is not fully clear to us presently the physical origin of this strong discrepancy.

One relevant question that always arises when using smooth inelastic hard particles to model granular gases, is to what extent the results would be modified if more realistic models, including for instance rotational degrees of freedom GNyB05 and/or velocity dependent restitution coefficients ByP04 , were considered. Although the hydrodynamic Navier-Stokes equations including these effects are known in some limiting cases, their analysis is far beyond our present reach. Nevertheless, it can be expected that an extension of the energy balance appearing in the simple case discussed here will hold when more dissipation mechanisms are included. Then it is our conjecture that a similar shear state but including the new rotational effects will also show up.

## Vi Acknowledgements

This research was supported by the Ministerio de Educación y Ciencia (Spain) through Grant No. FIS2008-01339 (partially financed by FEDER funds).

## Appendix A Navier-Stokes transport coefficients

In this Appendix, the explicit expressions of the transport coefficients and the cooling rate introduced in Eqs. (4)-(10) are given. The elastic shear viscosity and heat conductivity are

(47) |

(48) |

respectively. As usual in the context of granular fluids, the Boltzmann constant has been set equal to unity. The factors accounting for the dependence on the restitution coefficient are given by

(49) |

(50) |

(51) |

(52) |

In the above expressions,

(53) |

(54) |

(55) |

It is easily verified that and tend to unity when goes to one, while and vanish in this limit.

## References

- (1) C.S. Campbell, Annu. Rev. Fluid Mech. 22, 57 (1990).
- (2) N.V. Brilliantov and T. Pöschel, Kinetic theory of granular gases (Oxford University Press, Oxford, 2004).
- (3) J.W. Dufty, J. Phys.: Condens. Matter 12, A47 (2000).
- (4) I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993); I. Goldhirsch, M.L. Tan, and G. Zanetti, J. Sci. Comp. 8, 1 (1993).
- (5) S. McNamara and W.R. Young, Phys. Rev. E 50, R28 (1994); ibid. 53, 5089 (1996).
- (6) A. Puglisi, M. Assaf, I. Fouxon, and B. Meerson, Phys. Rev. E 77, 021305 (2008).
- (7) E. Efrati, E. Livne, and B. Meerson, Phys. Rev. Lett. 94, 088001 (2005).
- (8) B. Meerson, I. Fouxon, and A. Vilenkin, Phys. Rev. E 77, 021307 (2008).
- (9) J.J. Brey, J.W. Dufty, C.S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998).
- (10) J. Wakou, R. Brito, and M.H. Ernst, J. Stat. Phys. 107, 3 (2002).
- (11) C.K.K. Lun, S.B. Savage, D.J. Jeffrey, and N. Chepurniy, J. Fluid Mech. 140, 223 (1984)
- (12) J.T. Jenkins and M.W. Richman, J. Fluid Mech. 192, 313 (1988).
- (13) N. Sela, I. Goldhirsch, and S.H. Noskowicz, Phys. Fluids 8, 2337 (1996).
- (14) P.K. Haff, J. Fluid Mech. 134, 401 (1983).
- (15) J.J. Brey, D. Cubero, F. Moreno, and M.J. Ruiz-Montero, Europhys. Lett. 53, 432 (2001).
- (16) R. Soto, M. Mareschal, and M. Malek Mansour, Phys. Rev. E 62, 3836 (2000).
- (17) J.J. Brey and D. Cubero, in Granular Gases, edited by T. Pöschel and S. Luding (Springer-Verlag, Berlin, 2001).
- (18) I. Goldhirsch, Annu. Rev. Fluid Mech. 35, 267 (2003).
- (19) J.J. Brey, M.J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 54, 3664 (1996).
- (20) J.J. Brey, M.J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 60, 3150 (1999).
- (21) G.A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon Press, Oxford, 1994).
- (22) A.L. García, Numerical Methods for Physics (Prentice Hall, Englewood Cliffs, NJ, 2000).
- (23) J.F. Lutsko, Phys. Rev. E 63, 061211 (2001).
- (24) J.J. Brey, M.J. Ruiz-Montero, and F. Moreno, Phys. Rev. E 69, 051303 (2004).
- (25) A. Santos, V. Garzó, and J.W. Dufty, Phys. Rev. E 69, 061303 (2004).
- (26) A. Domínguez, P. Maynar, M.I. García de Soria, J.J. Brey, and M.J. Ruiz-Montero, unpublished.
- (27) I. Goldhirsch, S.H. Noskowicz, and O. Bar-Lev, J. Phys.: Condens. Matter 17, S2591 (2005).