Steady uniform shear flow in a low density granular gas

Steady uniform shear flow in a low density granular gas J. J. Brey, M. J. Ruiz-Montero, and F. Moreno Fı ´sica Teo ´rica, Facultad de Fı ´sica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain ~Received 8 October 1996! The steady uniform shear flow of a low density granular sytem is studied by means of a kinetic model equation and also by direct Monte Carlo simulation of the Boltzmann equation. The former can be exactly solved for arbitrary shear rate and dissipation. Explicit expressions for the one-particle distribution function and the pressure tensor are obtained. Comparison of the results with those of previous theories is presented in the appropriate limits. The simulation shows that the model reproduces fairly well the values of the stresses and, in particular, the phenomenon of normal stress differences. The agreement is also very good for the velocity distribution function in the thermal velocity region, although significant discrepancies appear for large velocities. @S1063-651X~97!12803-3# PACS number~s!: 05.20.Dd, 47.50.1d, 47.20.2k I. INTRODUCTION The development of kinetic model equations describing the dynamics of particles, which collide inelastically, seems to be an important step towards the understanding of the complex phenomena taking place in rapid granular flows @1–3#. Nevertheless, the nonconservation of energy in collisions introduces quite drastic changes in the physics of the problem, and the generalization of the existing theories for ordinary fluids is far from being an easy task. This refers to both the own derivation or proposal of a kinetic equation for the distribution function of the system and also to solving the equation to obtain consistent hydrodynamic-like equations describing the macroscopic evolution. Although a hydrodynamic description, in terms of the density, the macroscopic flow velocity, and the granular temperature is suggested by analogy with normal fluids and also by experiments and computer simulations, a derivation from a more fundamental description of the system is needed. In the last years, several kinetic equations for the oneparticle distribution function of rapid granular flows have been proposed, extending the heuristic arguments leading to the Boltzmann and Enskog equations for ordinary gases @4–7#. In addition, a more basic approach to the problem, starting from the Liouville equation, has been presented @8#. Nevertheless, even if it is assumed that these equations are relevant to the description of rapid granular flows, a fundamental difficulty arises. The standard procedure of obtaining solutions to kinetic equations is by means of the ChapmanEnskog method @9#, in which an expansion in the gradients of the hydrodynamic fields is carried out. For molecular fluids, the reference state about which gradients are considered is the Maxwellian local equilibrium distribution. The corresponding reference distribution for inelastic particles is complex and has not been determined to date. As a consequence, the limit of small dissipation ~i.e., the coefficient of restitution a close to 1!is usually considered. Nevertheless, many of the phenomena peculiar to rapid granular fluids are clearly associated to nonlinear and rheological effects and to not very small dissipation. This includes properties such as significant normal stress differences under shear @2,10,11#, spontaneous formation of dense clusters surrounded by regions of low density @12–15#, and inelastic collapse @13#. The reasons mentioned above indicate that it is instructive to consider simple model kinetic equations of inelastic gases which allow controlled and detailed analysis. Kinetic models have proven to be very useful in the study of far from equilibrium states of dilute molecular gases @16#. For many physical situations, exact solutions to the models have been derived, and comparison with numerical solutions of the exact equations obtained by computer simulation shows a fairly good agreement. Very recently @8#, two kinetic models for systems of inelastic hard spheres have been proposed. Both models are formulated as approximations of the revised Enskog equation, but with different quantitative accuracy. Here we will consider the limit of low densities and large length scales as compared with the diameter of the particles. In this limit, the two models lead to the same equation, and become a kinetic model of the Boltzmann equation for inelastic hard particles. Most of the simplest far from equilibrium physical situations correspond to steady states. In this paper, we study an unbounded, steady, and uniform shear flow ~USF!. While in molecular fluids such a state is not possible due to viscous heating, in granular fluids this effect can be balanced by dissipation in collisions. The simplicity of the model allows us to obtain the exact solution describing the steady USF and also all the physically relevant quantities. This is done without introducing any expansion in the gradients or in the inelasticity of collisions. The results for the pressure tensor show anisotropy, i.e., normal stress differences. It must be noted that our approach is quite different from others, in which a specific form is assumed for the oneparticle distribution function describing a given state. For instance, in the theory by Jenkins and Richman @17#a generalized Maxwellian is conjectured to model the distribution function of the steady USF. The velocity-independent coefficients are identified through the balance equation for the second order fluctuating velocity correlation tensor. This theory, which in principle also keeps all orders in the shear rate and dissipation, leads to normal stress differences which are in good agreement with the results of molecular dynamics simulations. On the other hand, the model we use has been formulated for arbitrary conditions, and no specific apPHYSICAL REVIEW E MARCH 1997VOLUME 55, NUMBER 3 55 1063-651X/97/55~3!/2846~11!/$10.00 2846 © 1997 The American Physical Society proximation is introduced to apply it to the steady USF. To obtain the one-particle distribution function one has to solve the kinetic model equation under the appropriate conditions. Another interesting approach to the same problem @18# has been carried out by expanding the Boltzmann equation to Burnett order. The perturbative expansion method is specific for the steady USF, and it is based on the observation that in this state the shear rate scales as A e , where e 512 a 2. The theory results in a prediction for the normal stress differences which are very close to those of Jenkins and Richman, when the latter are expanded to the same order of approximation. The direct Monte Carlo simulation method @19#was developped to obtain numerical solutions to the Boltzmann equation, and has been successfully applied to a great variety of situations. Since it is formulated for an arbitrary scattering law, there is no difficulty in using it to study granular flows @20#. Here we present results obtained for the steady USF, and compare them with the predictions of the model kinetic equation. For the components of the pressure tensor the agreement turns out to be very good over a wide range of values of the coefficient of restitution. Regarding the oneparticle distribution function, the agreement is only semiquantitative. This is not surprising, since the detailed Boltzmann collision operator is replaced in the model by a much simpler effective term, which is determined by the second moments of the one-particle distribution. The structure of the paper is as follows. In Sec. II the model is briefly discussed, and the basic equations are presented. In addition, a further simplification to the formulation in Ref. @8#in introduced. It consists of an evaluation of the dissipation source term in the local equilibrium approximation. We believe this is consistent with the spirit of the model, and makes the calculations much simpler. The model kinetic equation is particularized for the USF in Sec. III, and the steady pressure tensor is computed in Sec. IV. The results are compared to the theories of Jenkins and Richman @17# and of Sela, Goldhirsch, and Noskowitz @18#. Section V deals with the one-particle distribution function, which is also compared with the above theories in the low dissipation limit. In Sec. VI we discuss the Monte Carlo simulation results and, finally, Sec. VII provides a short summary and conclusions. II. BASIC EQUATIONS OF THE MODEL We consider a gas of identical smooth disks (d52) or spheres (d53) of diameter s and mass m, whose collisions are characterized by a constant coefficient of normal restitution a in the interval (0,1#. The Boltzmann equation for the one-particle distribution function, f(r,v,t), is @7,8# S ] ] t1v1•“1 D f~r1,v1,t!5JB@r1,v1,t u f~t!#,~1! where JBis the ~inelastic!Boltzmann collision operator JB@r1,v1,t u f~t!#5 s d21 E dv2 E d s ˆQ~g• s ˆ!~g• s ˆ! 3@ a 22f~r1,v1 8,t!f~r1,v2 8,t! 2f~r1,v1,t!f~r1,v2,t!#.~2! In the above expression, s ˆis the unit vector pointing from the center of particle 2 to the center of particle 1 at contact, Qis the Heaviside step function, g5v12v2, and v1 85v1211 a 2 a ~g• s ˆ! s ˆ,~3! v2 85v2111 a 2 a ~g• s ˆ! s ˆ~4! are precollisional velocities leading after collision to velocities v1and v2. The local particle number density n, flow velocity u, and temperature Tare defined in the usual way, n~r,t!5 E dvf~r,v,t!,~5! n~r,t!u~r,t!5 E dvvf~r,v,t!,~6! d 2n~r,t!kBT~r,t!5 E dv1 2m@v2u~r,t!#2f~r,v,t!.~7! From Eq. ~1!the following evolution equations for these fields are easily obtained: ] n ] t1“•~nu!50, ~8! ] u ] t1u•“u1~nm!21“•P50, ~9! d 2nkB ] T ] t1d 2nkBu•“T52~¹u!:P2“•q2~12 a 2! v , ~10! where P~r,t!5 E dvm~v2u!~v2u!f~r,v,t!~11! is the pressure tensor, q~r,t!5 E dvm 2~v2u!2f~r,v,t!~12! is the heat flux, and v ~r,t!5m s d21 8 E dv1 E dv2 E d s ˆQ~g• s ˆ! 3~g• s ˆ!3f~r,v1,t!f~r,v2,t!~13! 55 2847STEADY UNIFORM SHEAR FLOW IN A LOW DENSITY . . . is a source term describing the rate of dissipation in collisions. Solving Eq. ~1!for a given situation is a formidable task. Even in the case of ordinary gases, very little is known about solutions of the Boltzmann equation for far from equilibrium states, and the presence of dissipation complicates in a nontrivial way the structure of the Boltzmann equation. The simplest physical situation one can think of for a granular fluid is the so-called homogeneous cooling state ~HCS!, for which the system is uniform and the time dependence occurs entirely through the temperature. Although it has been shown @7,21#that Eq. ~1!admits a solution describing such state, the exact form of the distribution function is not known. For these reasons, it is worth looking for model kinetic equations similar to those which have proven to be so fruitful for ordinary fluids @16#. The basic idea is to replace the Boltzmann collision operator with a simpler form, while preserving the essential properties of the exact equation. Very recently @8#, a kinetic model along the lines sketched above has been proposed. The Boltzmann collision operator is separated into two parts: one in the velocity subspace spanned by 1, v, and v2, and the other one in the subspace orthogonal to it. The first contribution is retained exactly, while the second one is approximated by a single relaxation time term. This guarantees that the balance equations ~8!– ~10!are preserved by the model and that the fluxes Pand q, and also the source term w, are the same functionals of the distribution function fas in the Boltzmann equation, i.e., they are given by Eqs. ~11!–~13!. Since the details are described elsewhere @8,22#, we only give the result here. The Boltzmann collision operator is replaced by JB→2 n ~f2fl!21 nkBTfl c ~V!~12 a 2! v ,~14! where we have introduced the peculiar velocity V5v2u, c ~v!5mv2 dkBT21, ~15! and flis the local equilibrium distribution, fl~r,v,t!5n~r,t! F m 2 p kBT~r,t! G d/2 exp F 2mV2~r,t! 2kBT~r,t! G . ~16! Also, n @n(r,t),T(r,t)#is an effective local collision frequency, specified in more detail below. The model defined by Eq. ~14!, although much simpler than the original Boltzmann equation, is still too complicated to allow exact calculations in most of the applications. This is due to the nonlinear functional dependence of v on the distribution function. Then, we go a step further in the simplification, and approximate v by its local equilibrium value, i.e., v ~r,t!→ v l~r,t!5m s d21 8 E dv1 E dv2 E d s ˆQ~g• s ˆ! 3~g• s ˆ!3fl~r,v1,t!fl~r,v2,t! 5~d21! S p m D 1/2 n2 s d21~kBT!3/2.~17! In summary, our model kinetic equation is given by S ] ] t1v•“ D f52 n ~f2fl!21 nkBTfl c ~V!~12 a 2! v l. ~18! Trivially, this equation also leads to the balance equations ~8!–~10!, with the only difference that in the last of them the source term v is replaced by v l. In the elastic collision limit, the kinetic model reduces to the well-known Bhatnagar-Gross-Krook ~BGK!equation @23#. Application of Eq. ~18!to the HCS is very simple. The solution is given by a Maxwellian with a time dependent temperature. The time evolution of this latter quantity is given by Eq. ~10!, which for the HCS has the solution T~t!5T~0! S 11t t0 D 22 ,~19! with t0 215~12 a 2!d21 d s d21n S p kBT~0! m D 1/2 .~20! As we have already mentioned, the exact solution of the Boltzmann equation for the HCS is not known and, in particular, it is easy to check that it is not a Maxwellian ~except in the elastic limit, a 51). Nevertheless, both theoretical studies @7,21#and Monte Carlo simulations @20#have shown that deviations from Maxwellian are quantitatively small for thermal velocities and Eq. ~19!provides a very accurate approximation for the time evolution of the temperature. III. UNIFORM SHEAR FLOW We want to study the uniform shear flow ~USF!, which is characterized by a constant linear velocity profile u~r!5a•r,~21! where aij5a d ix d jy ,abeing the constant shear rate. For simplicity we consider an unbounded system. In addition, the heat flux vanishes and the density n(0) and temperature T(0) are uniform. For this flow, Eqs. ~8!,~9!, and ~10!imply that the density is constant in time, the pressure tensor is uniform, and the temperature obeys the energy balance d 2n~0!kB ] T~0! ] t52aPxy ~0!2~12 a 2! v l ~0!.~22! .This equation clearly shows the quite different behavior of elastic and inelastic fluids under shear. While for molecular fluids the USF is always a time-dependent state whose temperature increases monotonically in time due to viscous 2848 55J. J. BREY, M. J. RUIZ-MONTERO, AND F. MORENO heating, in granular media the variation of the temperature arises from the balance of two opposite effects: viscosity and dissipation in collisions. As a consequence, a steady state is possible when both effects cancel each other. Our aim in the following will be to analyze such a state. It is worth mentioning that sheared granular media posses multiple steady states, depending on the initial and boundary conditions @24#. Here we restrict ourselves to the USF as described above, assuming that the system has reached such a state, and without paying any attention to its stability. A detailed analysis of these questions will be published elsewhere. Let us introduce position and velocities coordinates with respect to a frame moving with the flow velocity u.We define the Lagrangian coordinates by R5L~t!•r,V5v2a•r,~23! where Lij~t!5 d ij2aijt.~24! In terms of these, Eq. ~18!becomes ] f ] t1LijVj ] f ] Ri 2aijVj ] f ] Vi 52 n ~f2fl!21 nkBTfl c ~12 a 2! v l. ~25! Since the USF is macroscopically homogeneous in the Lagrangian frame, we expect the associated distribution function, f(0), to have the same property, i.e., to be independent of R. Therefore, Eq. ~25!reduces to ] f~0! ] t2aijVj ] f~0! ] Vi 52 n ~0!~f~0!2fl ~0!!21 n~0!kBT~0!fl ~0! c ~0! 3 ~12 a 2! v l ~0!.~26! This equation implies that ] ^ Vi & ~0! ] t52a d ix ^ Vy & ~0!,~27! where the angular brackets denote averaging with respect to f. Therefore, if the initial state verifies ^ V & (0)50, this property is kept in time by Eq. ~26!. In the following we will limit our consideration to initial conditions of this kind. This is the reason why, for velocities in the Lagrangian frame, we used the same symbol as for the peculiar velocities introduced in Sec. II. Then the local equilibrium distribution appearing in Eq. ~26!has the form fl ~0!~V,t!5n~0! S m 2 p kBT~0!~t! D d/2 exp S 2mV2 2kBT~0!~t! D , ~28! with all the possible time dependence occurring through the temperature. IV. PRESSURE TENSOR By multiplying Eq. ~26!by ViVjand integrating over V,it is straightforward to obtain the evolution equations for the components of the pressure tensor in the USF, ] ] tPij ~0!1aikPjk ~0!1ajkPik ~0!52 n ~0!~Pij ~0!2 d ijp~0!! 22 d d ij~12 a 2! v l ~0!,~29! where p5(iPii /d5nkBTis the hydrostatic pressure. For a 51 we recover the equations for an ordinary gas under USF in the BGK approximation, which have been extensively studied @25–27#. From Eq. ~29!, in particular, one obtains ] ] tp~0!52 2 daPxy ~0!22 d~12 a 2! v l ~0!,~30! ] ] tPxy ~0!52aPyy ~0!2 n ~0!Pxy ~0!,~31! ] ] tPyy ~0!52 n ~0!~Pyy ~0!2p~0!!22 d~12 a 2! v l ~0!.~32! This closed system of equations has the trivial solution Pxy (0)5Pyy (0)5p(0)50. A second steady solution, denoted by an asterisk, is given by Pxy *52 a n * F p*22 d~12 a 2! v l * n * G ,~33! Pyy *52 n * aPxy *,~34! a2 n * F p*22 d~12 a 2! v l * n * G 2~12 a 2! v l *50. ~35! For given values of the shear rate a, the number density n, and the restitution coefficient a , Eq. ~35!determines the value of the pressure ~or temperature!at which the steady state is possible. It is convenient at this point to specify the form of the effective collision frequency n . Dimensional analysis requires that n 5cn s d21 S p kBT m D 1/2 ,~36! where cis a dimensionless constant. We are going to fix this by requiring that the model gives the correct ~Boltzmann! value for the Navier-Stokes shear viscosity in the elastic limit a 51. This leads to c516/5b, with b.1.016, for d53~spheres!, and c52/b, with b.1.022, for d52~disks! @28#. When Eq. ~36!is used in Eqs. ~33!–~35!they reduce to P ˜ xy52 H ~12 a 2!~d21! c F 122~d21! dc ~12 a 2! G J 1/2 , ~37! 55 2849STEADY UNIFORM SHEAR FLOW IN A LOW DENSITY . . . P ˜ yy52 P ˜ xy a ˜ 5122~d21! dc ~12 a 2!,~38! a ˜ 25d~d21!~12 a 2! dc22~d21!~12 a 2!.~39! We have introduced the dimensionless pressure tensor P ˜ ij and shear rate a ˜ in the steady state as P ˜ ij5Pij * p*,a ˜ 5a n *.~40! Equation ~39!shows that in the steady state the ratio of the shear rate to the collision frequency depends only on the restitution coefficient. The simple dependence of the presure tensor on the shear rate found here contrasts with the complex one found in the time-dependent uniform shear flow of ordinary gases @25–27#. Let us suppose we carry out a series of experiences, all of them in the same system, i.e., with given values of a and n, but with different values of the shear rate, measuring the shear viscosity in the steady state. Taking into account that Eq. ~39!implies that T*}a2,it follows from Eq. ~37!that the generalized shear viscosity h *is linear in the shear rate, h *~a![2Pxy * a}a.~41! Of course, the proportionality constant in this expression depends on the coefficient of restitution. The other components of the pressure tensor are easily computed by means of Eq. ~29!, particularized for the steady state. One obtains P ˜ xx5112~d21!2~12 a 2! dc ,~42! and, for a three-dimensional system, P ˜ zz5P ˜ yy ,P ˜ xz5P ˜ yz50. ~43! Let us note that the reduced pressure tensor in the steady state is expressed as a function of a alone, being independent, in particular, of the applied shear rate. Then, for instance, for the normal stress ratios we have Pxx * Pyy *5dc12~d21!2~12 a 2! dc22~d21!~12 a 2!,~44! Pyy * Pzz *51. ~45! Therefore, the model predicts anisotropy of the diagonal terms of the pressure tensor in the shear plane, but not in the plane perpendicular to the flow. The existence of normal stress differences in shear flows of normal fluids is very well known. It is a viscometric effect which appears beyond the Navier-Stokes approximation @29#. For a low density USF they have been calculated using the BGK model @27#. There, differences between Pxx and Pyy but not between Pyy and Pzz were found. As we have already mentioned, an approximated theory for the steady USF of smooth inelastic disks has been formulated by Jenkins and Richman @17#. Their results for the second velocity moments in the low density limit read, in our notation, a ˜ ~JR!25~11 a !~12 a 2!~723 a !2 c2~917 a !,~46! P ˜ xy ~JR!52 4 2529 a @~917 a !~12 a !#1/2,~47! P ˜ yy ~JR!5917 a 2529 a ,~48! P ˜ xx ~JR!541225 a 2529 a .~49! Although these expressions may appear rather different from the corresponding results obtained from our model kinetic equation, Figs. 1 and 2 show that the discrepancies are quite small for low dissipation, and tend to vanish when a approaches unity. Let us also mention that recent molecular dynamics simulations @24#have found that the values of the steady temperature for a dilute gas of inelastic disks can be closely fitted by an empirical relation, which in our units reads a ˜ 22516Ac2 p S 1 12 a 21B D ,~50! with A.0.080 and B.20.54. The above expression has the same functional dependence on the restitution coefficient as Eq. ~39!. In addition, from this latter equation one identifies A.0.100 and B.20.511, which are surprisingly close to the simulation values, taking into account that no excluded volume effects are considered in our model, and also that the flow observed in Ref. @24#is highly nonuniform and it exhibits time-dependent microstructures. Sela, Goldhirsch, and Noskowitz @18#carried out a study of a two-dimensional granular USF to Burnett order using the Boltzmann equation. They arrived at expressions for the steady temperature and the pressure tensor of the forms a ˜ 25C~12 a 2!1D~12 a 2!3/21O~12 a 2!2,~51! P ˜ xy5E~12 a 2!1/21O~12 a 2!3/2,~52! P ˜ xx511F~12 a 2!1O~12 a 2!2,~53! P ˜ yy511G~12 a 2!1O~12 a 2!2.~54! Here C,D,E,F, and Gare dimensionless constants. In Table I we compare the values for these constants given in Ref. @18#with those obtained by expanding the expressions derived in this section. Also given are the results obtained from the expansion of the Jenkins and Richman theory. The values from the three theories are very close and, in particu2850 55J. J. BREY, M. J. RUIZ-MONTERO, AND F. MORENO lar, the agreement between our results and the Burnett expansion of the Boltzmann equation is excellent. V. VELOCITY DISTRIBUTION FUNCTION The formal solution of Eq. ~26!can be written as f~0!~V,t!5e2setaijVj ] ] Vif~0!~V,0! 1 E 0 tdt8e2~s2s8!e~t2t8!aijVj~ ] / ] Vi! F n ~0!~t8 ! 21 n~0!kBT~0!~t8 ! c ~0!~V,t8 !~12 a 2! v l ~0!~t8 ! G 3fl ~0!~V,t8 !,~55! where s[ E 0 tdt8 n ~0!~t8!~56! is a measure of the average number of collisions per particle between 0 and t,s85s(t8), and exp@taijVj( ] / ] Vi)‡is a shift operator in velocity space, etaijVj~ ] / ] Vi!h~V!5h@L ~2t!•V#.~57! Of course, the time-dependent temperature in Eq. ~55! must be consistently calculated, i.e., it is given by the solution of Eqs. ~30!–~32!. Let us now consider the long time limit of f(0)(t). Assuming that the steady state discussed in Sec. IV is reached, we know that the temperature tends to the constant value given by Eq. ~39!, and sdiverges. The conclusion is that f(0)(V,t) approaches the steady form FIG. 1. Reduced shear rate a ˜ as a function of the coefficient of restitution a for a twodimensional steady USF. The solid line corresponds to the present theory, and the dashed line to that of Jenkins and Richman. FIG. 2. The pressure tensor for a twodimensional steady USF. The solid lines correspond to the present theory, and the dashed ones to that of Jenkins and Richman. 55 2851STEADY UNIFORM SHEAR FLOW IN A LOW DENSITY . . . f*~V!5 E 0 `dte2 n *tetaijVj~ ] / ] Vi! 3 F n *21 n~0!kBT* c *~V!~12 a 2! v l * G fl *~V!, ~58! with c *~V!5mV2 dkBT*21. ~59! It can be checked that this expression reproduces the values of the steady pressure tensor found in Sec. IV, and also that it gives a vanishing heat flux. Now we introduce dimensionless quantities by V ˜ 5 S 2kBT* m D 21/2 V,~60! f ˜ 51 n~0! S 2kBT* m D d/2 f*.~61! Then Eq. ~58!reduces to f ˜ ~V ˜ !5 p 2d/2 E 0 `dse2sesa ˜ V ˜ y~ ] / ] V ˜ x! 3 F 12d21 c S 2 dV ˜ 221 D ~12 a 2! G e2V ˜ 2,~62! where use has been made of Eqs. ~17!and ~36!. In the case of a system of hard spheres, the marginal distribution for the zcomponent of the velocity is easily obtained, f ˜ z~V ˜ z![ E 2` `dV ˜ x E 2` `dV ˜ yf ˜ ~V ˜ !5 p 21/2 3 F 122~12 a 2! 3c~2V ˜ z 221! G e2V ˜ z 2.~63! This expression shows a limitation of our model kinetic equation, namely, that the distribution function takes unphysical negative values for large enough velocities. Nevertheless, the velocity values for which Eq. ~63!is negative are given by Vz 2. F 3c 2~12 a 2!11 G kBT* m,~64! and, therefore, the distribution function is positive in all the range of thermal velocities, even in the limit a →0. It is then possible that Eq. ~62!provides an accurate approximation of the solution of the Boltzmann equation in all the relevant region of the distribution function. This point will be analyzed in detail in Sec. VI. The marginal velocity distribution in the shear plane, also for d53, is f ˜ xy~V ˜ x,V ˜ y![ E 2` `dV ˜ zf ˜ ~V ˜ ! 5 p 21 E 0 `dse2se2[~V ˜ x1sa ˜ V ˜ y ! 2 1V ˜ y 2 ] 3 H 12 4 ~ 12 a 2 ! 3c@~V ˜ x1sa ˜ V ˜ y!21V ˜ y 2# J . ~65! In order to compare with previous theories, we have carried out a perturbative expansion of f ˜ for d52 in powers of 12 a 2. The result is f ˜ ~V ˜ !5 p 21e2V ˜ 2 F 12V ˜ 2sin2 u A c~12 a 2!1/2 1~122V ˜ 21V ˜ 2cos2 u 11 2V ˜ 421 2V ˜ 4cos4 u ! 312 a 2 c1O~12 a 2!3/2 G ,~66! Here, we have introduced a polar representation for the velocity, V ˜ (V ˜ , u ), where u is measured with respect to the streaming direction x. Jenkins and Richman @17#assumed that the steady distribution function has the form of a generalized Gaussian, f ˜ ~JR!51 p A u K u exp~2V ˜ •K21•V ˜ !,~67! where Kis a matrix which accounts for the anisotropy existing in the USF. When the low density limit of the above distribution is considered and the result expanded to order (12 a 2) one obtains f ˜ ~JR!~V ˜ !5 p 21e2V ˜ 2 F 12V ˜ 2sin2 u A 2~12 a 2!1/2 1~122V ˜ 212V ˜ 2cos2 u 11 2V ˜ 421 2V ˜ 4cos4 u ! 312 a 2 41O~12 a 2!3/2 G .~68! If we neglect the difference between c(;1.96) and 2, both expressions agree to order (12 a 2)1/2. The terms of order (12 a 2) have the same dependence on the velocity modulus TABLE I. Comparison of the values of the coefficients in expansions given by Eqs. ~51!–~54!obtained with different theories and in the present work. Model Sela, Goldhirsch, and Noskowitz Jenkins and Richman C0.511 0.511 0.522 D00 0 E20.714 20.714 20.707 F0.511 0.522 0.5 G0.511 0.522 0.5 2852 55 J. J. BREY, M. J. RUIZ-MONTERO, AND F. MORENO and also on u , but all the constant prefactors differ. In spite of this difference, let us notice that both expressions ~66!and ~68!make the same contribution to the diagonal part of the pressure tensor when the constant cis again approximated by 2. This can be verified by a direct calculation ~see also Table I!.The distribution function for the USF of a twodimensional granular gas has also been calculated at the same order by Sela, Goldhirsch, and Noskowitz in Ref. @18#. There, a much more complicated dependence on both polar components of the velocity is found. By comparing with Jenkins and Richman results, the authors conclude that the constants appearing as prefactors for the second harmonics in u in Eq. ~68!can be considered as rough averages of the corresponding functions obtained by them. We refer the reader to their paper for details. VI. DIRECT MONTE CARLO SIMULATION The direct simulation Monte Carlo method @19#has proved to be a very useful tool to obtain numerical solutions of the Boltzmann equation for molecular fluids. Very recently, it has also been applied to study the HCS of a low density granular flow @20#. Since the details of the method have been extensively discussed in Ref. @19#, they will not be given here. We saw in Sec. III that the distribution function for the USF becomes homogeneous in the Lagrangian frame. The consideration of homogeneous states allows a great simplification of the simulation. Therefore, we have carried out the simulation in the Lagrangian frame and restricted ourselves to solutions of the Boltzmann equation which stay homogeneous in that frame. In other words, we numerically solved the equation @compare with Eq. ~26!# ] f~0! ] t2aijVj ] f~0! ] Vi 5J B@V u f~0!#.~69! Of course, this implies that the possibility of spontaneous formation of spatial inhomogeneities is eliminated in the FIG. 3. Time evolution of the reduced shear rate a ¯ for a 50.8, and different values of the shear rate a. In all cases, the initial distribution was homogeneous with a Maxwellian velocity distribution in the Lagrangian frame. Time is measured in units of c„2 A pn (0)…21. FIG. 4. Steady reduced shear rate a ˜ as a function of the coefficient of restitution a for a dilute system of hard disks. The solid line corresponds to the kinetic model, and the symbols are results from the Monte Carlo simulation. 55 2853STEADY UNIFORM SHEAR FLOW IN A LOW DENSITY . . . simulation. Our aim here is to study the properties of the homogeneous steady state. Its stability will be analyzed elsewhere. For homogeneous systems there is no need to split the system into cells and, consequently, the spatial coordinates of the particles do not play any role in the simulation. In addition, no boundary conditions must be introduced. In our simulation we have considered a three-dimensional system. The number of particles is N51000, and the results have been averaged over 500 different trajectories. The time interval Dtover which it is assumed that free motion, including the effect of the inertial force, and collisions, are uncoupled has been taken d t50.025 n 0 21, where n 05(2 A p /c) n (t), with n (t) being the instantaneous value of the collision frequency given by Eq. ~36!. We considered a nonconstant time step in order to guarantee that it always remains much smaller than the average time between collisions, in spite of the change of the temperature @30#. The initial velocity distribution in all cases is a Maxwellian. In Fig. 3 we present the time evolution of a ¯ [a/ n (t) for a 50.8 and four different values of the shear rate, namely, a50.1, 0.25, 0.5, and 1. Time is measured in units of c„2 A pn (0)…21, where n (0) is the initial collision frequency. After an initial transient period, all curves converge to the same steady value, as predicted by Eq. ~39!. The same qualitative behavior has been found for all the reduced components of the pressure tensor Pij(t)/p(t). Therefore, in the following we will concentrate on the dependence of the steady values of the reduced quantities on the restitution coefficient a , once we have checked they do not depend on the shear rate. The results obtained for the steady reduced shear rate a ˜ for different values of a are shown in Fig. 4. The statistical errors are smaller than the symbols used to represent the data. Also plotted is the prediction of our model kinetic equation, i.e. Eq. ~39!with d53. It is seen that the agreement is remarkable at low dissipation, although the discrepancy increases as the restitution coefficient decreases. We FIG. 5. The same as in Fig. 4 for the reduced pressure tensor P ˜ ij. FIG. 6. Marginal velocity distribution function in the direction perpendicular to the shear plane. The symbols are simulation data, and the solid line corresponds to the model kinetic equation. The restitution coefficient is a 50.8. 2854 55 J. J. BREY, M. J. RUIZ-MONTERO, AND F. MORENO