High energy tail in the velocity distribution of a granular gas

High energy tail in the velocity distribution of a granular gas J. Javier Brey, D. Cubero, and M. J. Ruiz-Montero Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain ~Received 23 June 1998! The Enskog-Boltzmann equation for a homogeneous freely evolving system of smooth hard disks colliding inelastically is solved by means of the direct simulation Monte Carlo method. The distribution function shows an exponential high velocity tail, while it is Gaussian for small velocities. The numerical results are compared with recent predictions of approximate analytical theories and quite good agreement is found. @S1063-651X~99!01601-3# PACS number~s!: 81.05.Rm, 05.20.Dd, 51.10.1y I. INTRODUCTION The simplest realistic model for rapid granular flows is a system of smooth hard disks or spheres colliding inelastically @1#. Energy loss in collisions is described by means of a constant coefficient of restitution 0, a <1, which measures the decrease in magnitude of the normal component of the relative velocity of two colliding particles. It is trivial to formally modify the standard kinetic equations for molecular fluids, e.g., the Boltzmann and Enskog equations, adapting them to this dissipative dynamics. Of course, a different and nontrivial question is whether and to what extent inelasticity affects the range of validity of the kinetic equations. In this work it is assumed that the ~inelastic!Enskog-Boltzmann equation is an accurate tool to describe the time evolution of a granular gas with an arbitrary value of the coefficient of restitution. As a consequence of dissipation in collisions, there is no homogeneous steady state. The simplest possible state of a granular fluid is homogeneous but with a monotonically decreasing temperature @2#. This homogeneous cooling state ~HCS!plays an essential role in any perturbative approach to the description of rapid granular flows @3–5#. The velocity distribution function corresponding to the HCS is not known exactly, although a Gaussian or Maxwell distribution ~with a time dependent temperature!provides a good approximation for the region of thermal velocities, i.e., velocities of the order of v0(t)5@2kBT(t)/m#1/2 @6#. Deviations from Gaussian behavior in the HCS have been considered recently through different quantities. Goldshtein and Shapiro @7#computed the fourth velocity moment of the distribution function associated with the HCS from the inelastic Boltzmann equation and found a small correction to the Gaussian values. Although it is not an exact calculation since it is based in the first Sonine approximation, the method leads to results that are in quite good agreement with numerical solutions of the Boltzmann equation to be discussed later on. Let us point out that the expression given in Ref. @7#contains a numerical error @8,6#. More direct and qualitatively relevant evidence of non-Gaussianity has been given by Esipov and Po ¨schel @9#. They studied the asymptotic behavior, in the large velocity region, of the solution to the Boltzmann equation for the HCS and obtained that the decay of the distribution is not Gaussian but exponential. The aim of this paper is to present some numerical results on the high energy tail of the solution of the Enskog- Boltzmann equation for the HCS and compare them with the theoretical predictions. II. LARGE ENERGY TAILS The nonlinear Enskog equation for the one-particle distribution function f(v,t) of a homogeneous system of inelastic smooth hard spheres (d53) or disks (d52) is @7,10,11# ] tf~v1,t!5ge~ s ! s d21 E dv2 E d s ˆQ~g• s ˆ!~g• s ˆ! 3@ a 22f~v1 8,t!f~v2 8,t!2f~v1,t!f~v2,t!#.~1! Here s is the diameter of the particles, d s ˆis the element of the solid angle defined by the unit vector s ˆalong the line of centers of the colliding particles at contact, Qis the Heaviside step function, and g5v12v2. The primes on the velocities denote precollisional velocities yielding v1and v2as postcollisional ones. They are given by v1,2 85v1,271 2 a 21~11 a !~ s ˆ•g! s ˆ.~2! The only difference between Eq. ~1!and the Boltzmann equation for inelastic dilute gases is in the presence of the equilibrium pair correlation function at contact ge.Itisa function of the density that takes into account the increase of the collision frequency due to the finite size of the particles. By properly scaling the time, the solutions of both the Boltzmann and the Enskog equations for the HCS are the same. Goldshtein and Shapiro @7#noticed that the normal solution to Eq. ~1!describing the HCS has the form fH~v,t!5nv0 2d~t! f S v v0~t! D .~3! Substitution of this expression into the Boltzmann equation provides a differential equation for the function f (c), which becomes closed if the evolution equation for the temperature is used. The equation has not been solved exactly up to now and only approximate solutions are known. Here we are interested in the high energy tail of the velocity distribution, which is determined by the asymptotic behavior of f (c) for PHYSICAL REVIEW E JANUARY 1999VOLUME 59, NUMBER 1 PRE 59 1063-651X/99/59~1!/1256~3!/$15.00 1256 ©1999 The American Physical Society large c. Esipov and Po ¨schel @9#realized that for large velocities the gain term in the Boltzmann equation is negligible compared to the loss term. Thus f ~c!;Aexp~2Ac!,~4! is easily obtained for c→`, where Ais an undetermined integration constant and Ais another constant introduced when solving the differential equation by separation of variables. The behavior in Eq. ~4!represents an enhanced population for large energies compared to the usual Gaussian law. van Noije and Ernst @8#have obtained an approximated value for Aby solving the Enskog-Boltzmann equation in the first Sonine approximation. The result is A21~ a !5 G S d11 2 D ~12 a 2! S 113 16a2 D A 2G S d 2 D d ,~5! where a2is a quantity related with the fourth moment of f and therefore with the fourth velocity moment of the one particle distribution fH. In the first Sonine approximation and neglecting contributions to the velocity moments that are nonlinear in a2it is found that @7,8# a2516~12 a !~122 a 2! 9124d18 a d241 a 130~12 a ! a 2.~6! The above expression for a2has been compared with results obtained by means of the direct simulation Monte Carlo method and excellent agreement has been found @6#. This comparison is extended in the following to Eqs. ~4!and ~5!. The direct simulation Monte Carlo method was introduced as an efficient way of obtaining numerical solutions of the nonlinear Boltzmann equation @12#. The only modifications needed to apply it to the case of inelastic collisions are the expressions of the postcollision velocities. We refer to the literature for details of the method and only specific information about our simulations will be given here. The first point to be stressed is that we have solved the homogeneous Boltzmann equation, as given in Eq. ~1!. In the technical language of the direct simulation method, we have considered a single cell since positions of the particles do not appear in the equation and therefore there is no need to compute and store them. For the same reason, no boundary conditions are needed. Of course, this implies that the possibility of spontaneous formation of spatial inhomogeneities is eliminated in the simulation. In fact, it is the capacity of the method for forcing the system to stay in the HCS, even for small values of the coefficient of restitution, that renders it feasible to identify and measure the high velocity tails. In our simulations we have considered a two-dimensional system of 106particles. In addition, the results have been averaged over a number of different trajectories, going from 700 for a 50.1 up to 2500 for a 50.7 and 0.9. First, let us consider the case a 50.1. In Fig. 1 we present the simulation results for the scaled distribution function f as a function of the scaled velocity c5v/v0(t). Here and in the following the velocity chas been discretized in intervals Dc50.05. It is observed that ln f (c) shows a linear profile for c*2, indicating an exponential velocity tail, as predicted by Eq. ~4!. Also plotted for comparison is the Gaussian f M 5 p 21exp(2c2). The discrepancy between both functions is clear in the exponential region of f . Similar behaviors have been observed for a 50.2, 0.4, and 0.6. The main difference is that as the value of a increases the velocity range for which the exponential behavior shows up moves towards larger velocities. Thus the relative population corresponding to the asymptotic velocity tail decreases and the simulation noise increases, rendering it very hard to detect the exponential decay law. The highest value of a for which we have been able to obtain a relevant indication of exponential behavior is a 50.7. The results are shown in Fig. 2, where again the Gaussian has been included for reference. The quantitative deviation of f from the Gaussian in the velocity range considered is much smaller than in Fig. 1 and in fact it is quite difficult to identify a linear region in the plot of ln f (c). In order to get some additional information, we also present the FIG. 1. Scaled velocity distribution f as a function of the dimensionless velocity c5v/v0(t), where v0(t) is the thermal velocity at time t. The coefficient of restitution is a 50.1. The solid line is from the Monte Carlo simulation and the dashed line is the Gaussian f M5 p 21exp(2c2). FIG. 2. Same as in Fig. 1 but for a 50.7. Besides the numerical slope dln f /dc ~crosses!and dln f M/dc ~dashed line!are shown in the inset. The horizontal line is the theoretical prediction from Eqs. ~5!and ~6!. PRE 59 1257BRIEF REPORTS values of the slope of ln f and the corresponding Gaussian values. Although the error of the slope becomes very large for c*4 due to the very small population for those velocities for which f is of the order of 1028, a tendency towards a constant value of the slope is clearly identified. Moreover, this value seems to be very close to the theoretical prediction of Eqs. ~5!and ~6!, also indicated in the figure. The fact that the beginning of the exponential region shifts towards higher velocities as the system becomes more elastic is easy to understand since we know that the equilibrium velocity distribution is exactly Gaussian for all velocities in the elastic limit. More quantitatively, the asymptotic analysis carried out in Ref. @9#suggests that the exponential behavior requires c@(12 a 2)21. For a .0.7 we could not reach such a region, although in all cases the results were compatible with an exponential tail, in the sense that a tendency towards a constant slope of the logarithm of the distribution as the velocity increases was observed. In Fig. 3 we compare the simulation values for the exponent A( a ) introduced in Eq. ~4!with the theoretical results given by Eqs. ~5!and ~6!. The former have been obtained by fitting the numerical data for ln f in the velocity region where a linear decay was found. In order to improve the statistics, each of the displayed points is the result of averaging over the value of the slope at ten different times, so that each trajectory has provided ten raw values of A. The arrow in the value corresponding to a 50.7 indicates that it is only a lower bound, as discussed above. The agreement between theory and simulation is quite remarkable. III. DISCUSSION We have verified the exponential high energy tail of the velocity distribution function of a granular gas described by the Enskog-Boltzmann equation in the homogeneous cooling state. Also the accuracy of Eqs. ~5!and ~6!for the constant A has been verified. Although the exponential region in the velocity distribution has been clearly identified only for a &0.7, the conclusions are expected to apply for all 0, a <1 since no singular values of a , aside from 0 and 1, appear in the Boltzmann equation. In this context, let us note that for a &0.7 the HCS is known to be highly unstable against spatial perturbations @13#. For this reason, it is a very difficult task to obtain the high energy tail from molecular dynamics simulation techniques. The system becomes inhomogeneous and departs from the HCS very soon so that there is no way of reaching a significant level of statistics. Although the same comment is in principle true for the complete ~inhomogeneous!Boltzmann equation @6#, we have taken advantage of the possibility of simulating and solving numerically the homogeneous Boltzmann equation ~1!, from which the HCS is a stable solution for an arbitrary value of the coefficient of restitution. A practical result following from the numerical resolution of the Boltzmann equation is that the transition from the Gaussian behavior to the exponential one takes place over quite a narrow interval of thermal velocities c. In fact, an accurate prediction of the velocity for which the exponential tail is observed is obtained by equating the slopes of the Gaussian and the exponential distributions. ACKNOWLEDGMENTS We are very grateful to Mathieu Ernst and T. P. C. van Noije for providing us with results prior to publication. This research was partially supported by Grant No. PB96-0534 from the Direccio ´n General de Investigacio ´n Cientı ´fica y Te ´cnica ~Spain!. @1#C. S. Campbell, Annu. Rev. Fluid Mech. 22,57~1990!. @2#P. K. Haff, J. Fluid Mech. 134, 401 ~1983!. @3#C. K. W. Lun, S. B. Savage, D. J. Jeffrey, and N. Chepurniy, J. Fluid Mech. 140, 223 ~1984!. @4#N. Sela and I. Goldhirsch, J. Fluid Mech. 361,41~1998!. @5#J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys. Rev. E58, 4638 ~1998!. @6#J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 54, 3664 ~1996!. @7#A. Goldshtein and M. Shapiro, J. Fluid Mech. 282,75~1995!. @8#T. P. C. van Noije and M. H. Ernst, Granular Matter 1,57 ~1998!. @9#S. E. Esipov and T. Po ¨schel, J. Stat. Phys. 86, 1385 ~1997!. @10#J. J. Brey, J. W. Dufty, and A. Santos, J. Stat. Phys. 87, 1051 ~1997!. @11#J. T. Jenkins and M. W. Richman, Arch. Ration. Mech. Anal. 87, 355 ~1985!; Phys. Fluids 28, 345 ~1986!. @12#G. Bird, Molecular Gas Dynamics ~Oxford University Press, New York, 1976!;Molecular Gas Dynamics and the Direct Simulation of Gas Flows ~Clarendon, Oxford, 1994!. @13#I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 ~1993!; S. McNamara and W. R. Young, Phys. Rev. E 53, 5089 ~1996!. FIG. 3. Slope of the logarithm of the high energy tail of the velocity distribution in the homogeneous cooling state. The crosses are from the Monte Carlo simulation of the Boltzmann equation and the solid line is the theoretical prediction given by Eqs. ~5!and ~6!. 1258 PRE 59 BRIEF REPORTS